Why ''int'' and ''integral'' return different answers

10 次查看（过去 30 天）
Han F
Han F 2022-9-10
编辑： Bruno Luong 2022-9-10
Hello, I am doing a simple integration:
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!

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2022-9-10
编辑：John D'Errico 2022-9-10
Ran in:
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
BETA Beta function. Y = BETA(Z,W) computes the beta function for corresponding elements of Z and W. The beta function is defined as beta(z,w) = integral from 0 to 1 of t.^(z-1) .* (1-t).^(w-1) dt. The arrays Z and W must be real and nonnegative. Both arrays must be the same size, or either can be scalar. Class support for inputs Z,W: float: double, single See also BETAINC, BETALN. Documentation for beta doc beta Other functions named beta codistributed/beta gpuArray/beta sym/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)
ans = 99.9361
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.
https://mathworld.wolfram.com/Jacobi-GaussQuadrature.html#:~:text=Numerical%20Integration%20Jacobi-Gauss%20Quadrature%20Jacobi-Gauss%20quadrature%2C%20also%20called,given%20by%20the%20roots%20of%20the%20Jacobi%20polynomials.
  4 个评论
显示 2更早的评论
Han F
Han F 2022-9-10
Yes, the answer of these gamma function is 99.9361
Han F
Han F 2022-9-10
Hi John. THank you for your answer. I think you catched the idea in my mind.
I found the discripency when I was trying to check if the integral of a beta distribution is 1 (of course is 1). If I calculate the equation by a beta function, the result is 99.9361, but when I do the integration by myself, the result is around 53, as the question described. That's way I posted the question here.
Anyway, thank you for your detailed answer, I learned a lot from you guys.
Thanks a lot!
Have a nice day!

请先登录，再进行评论。

更多回答（3 个）

Torsten
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.
  1 个评论
Han F
Han F 2022-9-10
Thank you for your answer. I think I should be careful about singularities when I use integral next time.

请先登录，再进行评论。

Walter Roberson
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
Bruno Luong 2022-9-10
编辑：Bruno Luong 2022-9-10
Ran in:
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')

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by