Can someone explain why I'm getting 0 as answer for the following integral?

Can someone explain why I'm getting 0 as answer for the following integral? I expect to get an answer of 1.38e-10 (which is what I get on my TI calculator). Why don't I get the same with Matlab?

y=@(x) exp(-2*3.62e9*x);

integral(y,0,Inf)

ans = 0

Answer 1

[Massimo Zanetti]

This is a case where quadrature formulas fail..

Indeed, by computing the (very simple) integral manually, one gets

integ_result = 1/(2*3.62e9) = 1.3812 e-10.

[Jordan Wiker]

So, what do I need to do so that I can use Matlab and be confident that I'm getting the actual correct answer when integrating?

[Massimo Zanetti]

I suggest using symbolic computation. See Steven Lord answer.

Answer 2

[Steven Lord]

The value of your integrand in double precision for all x values greater than about 1e-7 underflows to 0.

y([1e-8 1e-7 1e-6])
ans =
    3.60644663446566e-32     3.72216055251194e-315                         0

So when integral evaluates your integrand over the range [0, Inf] it's extremely likely that the value of the integrand is 0 at all the points (except x = 0) at which integral evaluates the integrand. When it tries to take a closer look at the value of the integrand near x = 0, it probably chooses points whose values may be nonzero, but as you can see from the example above it would need to look very, very closely in order for the value of the integrand to be different enough from 0 for it to "decide" that the function is not the zero function.

The problem is in your integrand. The exponential function exp(x) grows quickly, and when you start taking the exponential of a multiple of x it grows even more quickly.

ezplot(@exp)
figure
ezplot(@(x) exp(2*x))

The Y limit on the first plot is around 200, while the Y limit on the second is about 50000. [If you use exp(-x) or exp(-2*x) then the plot is simply reflected around x = 0.] You're computing the exponential of -3 billion, 620 million times x -- that curve drops to the ground like it was fired out of the Large Hadron Collider at a substantial fraction of the speed of light.

There are two potential solutions that I can think of for how to handle this. The first would be to rescale your problem -- if that coefficient is in units of millimeters per second (3620000000), for instance, rescale the problem so that the coefficient is instead in kilometers per second (3620.) The second would be to use Symbolic Math Toolbox.

syms x
y = exp(-2*3.62e9*x);
int(y, x, 0, Inf)