Search code examples
matlabprecisionintegral

Matlab double integral with heavily suppressed exponentials

I have recently been trying to calculate the double integral of the functionfun = @(v,x)(10^4)*0.648*(1+v*0.001).*( exp(-2.83./( 10^-8+(sqrt(1+2*v*0.001)).*(x.^2)) ) -1).*(exp(-(v.^2)*0.33)), in the range (-1000,1000) for v and (0,a) for x, where a is either a very large number or infinity. What I have found is that while in the case a = inf the value seems to be decently accurate (it reduces to a single integral which is less cumbersome to evaluate numerically), but if I integrate from 0 to a 10^9 and from 10^9 to infinity the integrals don't sum up to the correct value, with the latter one being zero. What I am really interested in is in the integral from 0 to 10^9, but these results make me wonder if I can trust it at all. In what I have done, I also had to use a large prefactor (10^200) in front of the function to "compensate" for the small numbers; otherwise the results were all nonsense. I have tried to use vpa, but with no success. What am I doing wrong?

Rob

Solution

  • Looks like your problem has to do with the different methods Matlab uses for different cases and the big numbers you are handling.

    We can see your function with ezsurf just to have an idea on how does it behave.

    So hint 1 is that the value is going to be a negative value, lets integrate over small limits to see an approximation on how much will it be.

    integral2(fun,-100,100,0,100)

%ans =
%  -5.9050e+04

    And assuming that the function tends to zero, we know the final value should be on the neighborhood.

    Now hint 2:

    integral2(fun,-1000,1000,0,100)

%ans =
%  -2.5613e-29

    This doesn't make much sense, by increasing the range of the limits the integral basically became zero. After checking the documentation of integral2

    'Method' — Integration method 'auto' (default) | 'tiled' | 'iterated'

    Integration method, specified as the comma-separated pair consisting of 'Method' and one of the methods described below.

    Integration Method Description

    'auto' For most cases, integral2 uses the 'tiled' method. It uses the 'iterated' method when any of the integration limits are infinite. This is the default method.

    'tiled' integral2 transforms the region of integration to a rectangular shape and subdivides it into smaller rectangular regions as needed. The integration limits must be finite.

    'iterated' integral2 calls integral to perform an iterated integral. The outer integral is evaluated over xmin ≤ x ≤ xmax. The inner integral is evaluated over ymin(x) ≤ y ≤ ymax(x). The integration limits can be infinite.

    Ok, so if we don't define a method it will use "tiled" if the limits are finite, and "interpolated" if they are infinte.

    Could it be that if the range is too big, the tiles created by the "tiled" method are too big to accurately calculate the integral? If that is the case then "iterated" should not have that problem, let's check

    integral2(fun,-1000,1000,0,100,'Method','iterated')

%ans =
%  -5.9050e+04

    Interesting, looks like we are into something. Let's try the original problem

    integral2(fun,-1000,1000,0,inf)

%ans =
%  -5.9616e+04

integral2(fun,-1000,1000,0,10^9,'Method','tiled')

%ans =
%  -2.1502e-33

integral2(fun,-1000,1000,0,10^9,'Method','iterated')

%ans =
%  -5.9616e+04

integral2(fun,-1000,1000,10^9,inf)

%ans =
%     0

    That looks better. So it looks like the 'tiled' method is the problem with your function because its characteristics and the size of the range of the limits. So as long as you use 'iterated' you should be ok.