improper integral

Hi!!I need help in the calculation of the integral of the function below

    function out = fun(u,a,b)
          out = exp(-a.* exp(2.*u)).*exp(-(u + b.^2).^2./2.*b.^2)/...
              (sqrt(2.*pi.*b.^2));
     end

a=0; b=1; *quad('fun',-Inf,Inf,[],[],a,b);

Is there any idea? where is the mistake?

[Andrew Newell]

What version of MATLAB are you using? The syntax you are using for QUAD is outdated.

[Asatur Khurshudyan]

I guess you have to calculate the principal value (by Cauchy) of that integral, otherwise you`ll have Nan`s or Infinity`s. I also have such a problem with 3D plotting improper integrals with infinite upper limit. Have anyone did something like that?

[Roger Stafford]

For a = 0 this becomes just the integral of a normal distribution density function over its full range, and of course it must then be equal to 1. It is when a is greater than zero that the integral becomes more interesting. For negative a the integral is divergent.

Answer 1

[Andrew Newell]

I don't think quad ever allowed infinite limits. quadgk does. The answer you need depends on how old your MATLAB version is. I am going to express an answer in more recent syntax. If you try

a=0; b=1; 
fun = @(u) exp(-a.* exp(2.*u)).*exp(-(u + b.^2).^2./2.*b.^2)/...
  (sqrt(2.*pi.*b.^2));
quad(fun,-Inf,Inf)

then you get this message:

Warning: Infinite or Not-a-Number function value encountered.

and the answer is NaN. The likely reason is that there is underflow in these exponential expressions:

>> fun(1000)

ans =

   NaN

However, if you plot this function, you can see that this function approaches zero very rapidly. For example,

>> fun(0)

ans = 

    0

Therefore, you can do this calculation:

quadl(fun,-100,100)

ans = 

    1.0000

This is 1 to an accuracy of about 10^-8.

Note: the fun = ... line is an anonymous function.

Answer 2

[Matt Fig]

quadgk(@(x) fun(x,0,1),-inf,inf)  % 0 is a, 1 is b.

But for a = 0, b = 1, we encounter an infinite number or nan. Are these numbers correct? If so, you could try this:

quadgk(@(x) fun(x,0,1),-10,10) % Exponentially decreasing values.

Or look at QUADL.

quadl(@(x) fun(x,0,1),-200,200,1e-15)