How to numerically evaluate this singular integral in integral2?

Hello,

I am trying to numerically evaluate a 2D integral with an integrand singular along an elliptical curve, please see the simple code below

a = 1;
b = 2;
c = 3;

k = @(x,y) 1./(sqrt(x.^2/a^2+y.^2/b^2)-c);

h = integral2(@(x,y) k(x,y) ,0,10,-10,10);

This yields warnings that "the result fails the global error test" which is evidently because the integrand is singular when

I am curious as to how this integral can be restated to make a numerical integration tractable?

Thanks!

[Torsten]

Can you prove that the integrals exists at all ? Because the integrand is singular on a complete ellipse if c > 0.

[Ryan]

It may not in-fact, I was trying to evaluate a more difficult integral and thought perhaps if I could evaluate this simpler one I could apply whatever method I used to the more complicated one.

[Walter Roberson]

format long g
a = 1;
b = 2;
c = 3;
k = @(x,y) 1./(sqrt(x.^2/a^2+y.^2/b^2)-c);
N = 20000;
M = 10;
Areas = zeros(1,M);

for K = 1 : M
    x = sort(rand(1,N) * (10-0) + 0) .';
    y = sort(rand(1,N) * (10-(-10)) + (-10));

    f = k(x, y);
    Areas(K) = trapz(y, trapz(x, f));
end

Areas

Areas = 1x10

129.808911190055 -1.11824424236561 88.0158960832912 145.159907352496 -11.9500169387748 21.6013629977695 67.9199992329145 39.8938824225843 155.213969423048 125.799087585423

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Random sampling is showing a wide range of values. It is unlikely that the integral converges.

[Matt J]

Let us make the change of variables, .

Then the integral becomes,

Now, converting to polar coordinates,

which is obviously a non-convergent integral if and only if for any . Since the region of integration is a rectangle of width 10/b and height 10/a, that will happen if c<norm([10/a,10/b]), which it is in this case.

[Ryan]

Thank you Walter and Matt, I was hoping there was a way to render it numerically integratable but as demonstrated it appears that it isn't convergent to begin with.

Answer 1

[Matt J]

As demonstrated in the comment above, the integral is theoretically non-convergent unless c>=norm([10/a,10/b]). Here is a further, numerical test:

a = 1;
b = 2;
c0 = norm([10/a,10/b]); %critical threshold

c=c0*1.0001; %slightly greater than threshold

k = @(x,y) 1./(sqrt(x.^2/a^2+y.^2/b^2)-c);
h = integral2(@(x,y) k(x,y) ,0,10,-10,10)

h = -60.2647

c=c0*0.9999;%slightly less than threshold

k = @(x,y) 1./(sqrt(x.^2/a^2+y.^2/b^2)-c);
h = integral2(@(x,y) k(x,y) ,0,10,-10,10)

Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.

h = -60.4654

[Ryan]

Thanks for that, this is an interesting insight into the behaviour of the integral