Integral of exponential- Incorrect result

Hello Mathworks community,

I have to integrate a function in Matlab and I follow 2 different ways to compare the results: sumation and integral. This a reduced example:

xsum = 0; 
dx = 10;
xmin = 0; xmax = 40000E3;
%Way 1----------------------------------------------------
for x = xmin: dx :xmax
    xsum  = xsum + ((((10^-5)*x)^10)* exp(-x/1000))*dx;
    %Or the same as: xsum  = xsum + ((((10.^-5).*x).^10).* exp(-x/1000))*dx;        
end
xsum
%Way 2----------------------------------------------------   
fun = @(x) ((((10.^-5).*x).^10).* exp(-x/1000));  
xint = integral(fun,xmin,xmax)

And the results:

xsum = 3.6288e-11

xint = 1.3259e-11

does anyone know where the error in the integration is? (I am using Matlab R2019a)

Thank you in advance!!

[Johan]

Seeing that the results converges for smaller xmax I would guess the integral function ignores very small values or something like it ?

xsum = 0; 
dx = 1;
xmin = 0; xmax = 400E3;
%Way 1----------------------------------------------------
for x = xmin: dx :xmax
    xsum  = xsum + ((((10^-5)*x)^10)* exp(-x/1000))*dx;
    %Or the same as: xsum  = xsum + ((((10.^-5).*x).^10).* exp(-x/1000))*dx;        
end
xsum

xsum = 3.6288e-11

%Way 2----------------------------------------------------   
fun = @(x) ((((10.^-5).*x).^10).* exp(-x/1000));  
xint = integral(fun,xmin,xmax)

xint = 3.6288e-11

[Monse F]

Hi Johan,

Thank you for your answer.

I also noticed that both results are the same when xmax is much more smaller, but not sure why. This function is just a part of a much longer function with 3 exponentials. I noticed that if I remove the exponentials, both result (from sumation and integral) are the same. My suspicion is that the integral of the exponential is causing problems. Any idea about this?

Answer 1

[Jan]

Reduce the tolerances:

format long
fun  = @(x) (((1e-5 * x) .^10) .* exp(-x / 1000));
xmin = 0;
xmax = 40000E3;
xint1 = integral(fun, xmin, xmax, 'AbsTol', 1e-12, 'RelTol', 1e-12)

xint1 = 
     3.628800302951016e-11

xint2 = integral(fun, xmin, xmax, 'AbsTol', 1e-14, 'RelTol', 1e-14)

xint2 = 
     3.628800000063254e-11

xint3 = sum(fun(xmin:0.1:xmax)) * 0.1

xint3 = 
     3.628800000000002e-11

[John D'Errico]

As Jan says, you can reduce the tolerances. In some cases, that will solve the problem, as long as double precision is capable of handing the arithmetic without floating point issues. At some point, if your powers were sufficiently high or the exponential sufficiently nasty or you are going too far out, double precision will fail you. Then a symbolic integral, or vpaintegral will help to resolve the issues.

[Monse F]

Thank you, Jan! Now it works.

I just had to adapt the values of 'AbsTol' and 'RelTol' for the original function I was working on.