Order of integration in double integral

I have an analytical expression for function , which is fairly complicated and uses many smaller functions, and I want to numerically calculate its double integral for different values of the constant a:

The function has one singularity from one of its terms but it is on the boundary of the area of integration. I have included image of that function for , x from 0 to 5000 and y from 0 to 8. Note that the large blue section is not equal to 0 but very close to 0.

In the code I used parfor loop to get values of the integral for a going from 0 to 8 and I have included image of that calculation. The blue line represents one order of integration (), and the red the other one (). However, the result is significantly different for different order of integration, which should not be the case.

f = @(x,y) ...
    
I1 = integral2(f, 0,inf, 0,a, 'method','iterated','AbsTol',absTol,'RelTol',relTol);              %dxdy
I2 = integral2(@(y,x) f(x,y), 0,a, 0,inf, 'method','iterated','AbsTol',absTol,'RelTol',relTol);  %dydx

I got no error/warning messages in either case. I’ve done similar calculation for different functions and I’ve got the correct results in all of them, so this may be somewhat of peculiar function for MATLAB. These calculations take a lot of computing time, even hours on clusters, and I’m sorry I couldn’t include the whole code.

What could be the cause of this erroneous calculation? How can I determine if any of this results is correct? Should I somehow divide the area of integration into smaller pieces? If so, what could be the best way of doing that?

Answer 1

[David Goodmanson]

Hi A.K.

If you are convinced that the problem arises from the square root singularity, one approach is a change of variable

y = a-z^2   dy = -2z dz

Then the integral

Int{0,a} Int{0,inf} f(x,y) / sqrt(a^2-y^2) dx dy   

becomes

Int{0,sqrt(a)} Int{0,inf} 2 f(x,a-z^2) / sqrt(2a-z^2) dx dz   

and the singularity goes away. It's possible that the z^2 dependence could cause problems if f(x,a-z^2) is highly oscillatory in what used to be the y variable, but since the domain of z integration is limited, hopefully that will not occur. Otherwise it would be necessary to split the z integration into two parts.

Incidentally, your initial expression for the integral has dx and dy reversed from the common notation.

[A. K.]

Hi David,

Thank you very much for your suggestion. Unfortunately, it did not work. The results are still the same as before, as it can be seen from the image bellow. Only in this calculation I used much larger tolerances so I could get just a quick overview of the curves, so there is some 'noise'.

[David Goodmanson]

Hi A.K.

Given that result, one could speculate that the singularity per se is not the source of the problem. It's interesting that things agree right up until a = 1 and diverge after that. I'm guessing that making the variable change y = aw, dy = a dw and integrating from 0 to 1 (a similar thing could be done with z = sqrt(a)*w ) would not make a difference either.

[A. K.]

Hi David,

You were right, nothing significantly different happened with that variable change. Do you maybe know of some instances where integral2 is known to give the incorrect result? I would like to compare my example to those and hopefully get some insight into what is happening in this case.

[David Goodmanson]

HI AK

All I have been able find is instances where integration gives a number and reverse order integration gives a result like NaN-cant-make-the-integration-interval-any-smaller. I don't have cases where both integration and reverse order integration give a number.

Not that I am suggesting this is involved, but the results as a function of 'a' are reminiscent of the fact that sqrt(1+a^2) has a convergent Taylor series for a<1 and a different, convergent Laurent series for a>1.

[A. K.]

Hi David,

Thank you very much for your trouble.