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

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!

Answer 1

[John D'Errico]

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.

[Torsten]

One even gets an analytical expression for the integral:

syms x
f = x^(-0.98)*(1-x)^(-0.98);
int(f,x,0,1)

ans = 

[John D'Errico]

Good point.

[Han F]

Yes, the answer of these gamma function is 99.9361

[Han F]

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!

Answer 2

[Torsten]

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.

[Han F]

Thank you for your answer. I think I should be careful about singularities when I use integral next time.

Answer 3

[Walter Roberson]

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.

[Han F]

Thank you for your detailed explanation. I have more understandings now.

Answer 4

[Bruno Luong]

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')