Why ''int'' and ''integral'' return different answers
13 次查看(过去 30 天)
显示 更早的评论
Hello, I am doing a simple integration:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1121585/image.png)
For integral:
zzz = @(x) x.^(-0.98).*(1-x).^(-0.98)
integral(zzz,0,1)
ans =
52.3164
While for int:
syms x
zzz = x.^(-0.98).*(1-x).^(-0.98)
int(zzz,x,0,1)
ans =
99.9361
I know 99.9361 is the right numer, but why integral returns a different number? (I prefer to use 'integral')
A big thank you in advance!
0 个评论
采纳的回答
John D'Errico
2022-9-10
编辑:John D'Errico
2022-9-10
As others have said, the problem with integral is because of the singularities. Numerical methods don't like singularities.
zzz = @(x) x.^(-0.98).*(1-x).^(-0.98);
fplot(zzz,[0,1])
So it is a good idea to avoid integral on functions like that. Unless, of course, you see that the integral you are trying to evaluate is just the complete beta function. It helps to know your special functions!
help beta
As such, you should see the integral is given very nicely by beta. No explicit integration was necessary. Don't forget to add 1 to the exponents to give as arguments to the beta function.
beta(0.02,0.02)
Alternatively, I could probably do the same computation using a Gauss-Jacobi quadrature, but that seems to be a bit of overkill for this problem, since beta already exists in MATLAB, and does what you want very nicely.
更多回答(3 个)
Torsten
2022-9-10
Your function has singularities at x = 0 and x = 1.
Integral doesn't manage to estimate the area under the curve (which goes to Infinity at both places) correctly.
Walter Roberson
2022-9-10
You should pass in abstol and reltol options to integral()
integral() subdivides intervals until the contribution appears to be small. But there are functions for which a small contribution adds up a lot.
Consider for example sum of 1/x as x approaches infinity, or the similar integral. For any given threshold you can find x such that all further 1/x will be less than the threshold. Let the threshold be eps(1) and you have the situation where if you were looping adding 1/n that each 1/n contribution would individually be less than a 1 bit difference in the floating point sum, so it would be reasonable to stop adding them. But integral 1/x is log(x) and as x approaches infinity that would go to log(infinity) which is infinite. This shows that numeric methods can be infinitely wrong if they do not happen to do the right transform to match the function.
It is not really a bug in the numeric integration, it is a limitation: for any numeric integration algorithm you can create functions that the algorithm cannot integrate accurately.
Bruno Luong
2022-9-10
编辑:Bruno Luong
2022-9-10
I think the problem is more subttle than the simple presence of singularities. Some "mild" singularity can be handlle correctly by integration. The problem here is that for power that is close to -1, where the integration goes to infinity and it is not mild at all. It's more like of making a transition to diverging integration. No wonder why integrationn function cannot deal it well.
Normally for all w < 1, singularity exists but one observe that integration starts to have difficulty for w < 0.15, otherwise it converges and gives an accurate answer.
fclose=@(w)beta(w,w);
g=@(w,x)x.^(w-1).*(1-x).^(w-1);
fi=@(w)integral(@(x)g(w,x),0,1);
fi=@(w)arrayfun(fi,w);
warning('off','MATLAB:integral:MinStepSize')
close all
figure
ezplot(@(x)g(0.4,x),[0,1])
w = linspace(0.01,1);
plot(w,fclose(w))
hold on
plot(w,fi(w))
set(gca,'YScale','log')
legend('beta (exact)','integration')
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!