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

8 views (last 30 days)
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!

Accepted Answer

John D'Errico
John D'Errico on 10 Sep 2022
Edited: John D'Errico on 10 Sep 2022
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.
  4 Comments
Han F
Han F on 10 Sep 2022
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!

Sign in to comment.

More Answers (3)

Torsten
Torsten on 10 Sep 2022
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 Comment
Han F
Han F on 10 Sep 2022
Thank you for your answer. I think I should be careful about singularities when I use integral next time.

Sign in to comment.


Walter Roberson
Walter Roberson on 10 Sep 2022
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.
  1 Comment
Han F
Han F on 10 Sep 2022
Thank you for your detailed explanation. I have more understandings now.

Sign in to comment.


Bruno Luong
Bruno Luong on 10 Sep 2022
Edited: Bruno Luong on 10 Sep 2022
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')

Community Treasure Hunt

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

Start Hunting!