Caluclating an integral over a rectangle, with a singularity point.

Hi! I am trying to evaluate the following integral:

where bx,by are some set values.

So I wrote the following code:

funcl_nn = @(x,z) (1/(4*pi*epsilon_zero))*(1./sqrt(x.^2+z.^2));
nn = quad2d(funcl_nn,-dist/2,dist/2,-bz/2,bz/2,'AbsTol',1e-8,'Singular',true,MaxFunEvals=100000);

and also tried to use integral2 and tiled method, but it always returns: Warning: Reached the maximum number of function evaluations (100000). The result fails the global error test.

What's weird is that if i change the values of by,bx to be larger, then the error does not occur.

I assume it is caused by the point (0,0).

I would try to switch to polar coordinates if the domain wasn't rectangular... Can anyone tell me what can I do such a situation? Thanks in advance!

[Alan Stevens]

For what limits does it not work?

[Matt J]

The function is not integrable, so it isn't clear why you expect a result other than inf from the get-go.

[Alan Stevens]

@Matt Yes, it's integrable. Maple gives the following (ignoring the 4pi.epsilon):

[Matt J]

Ah, well. Problem solved.

[Dyuman Joshi]

As @Alan Stevens asked, for what limits does it not work?

integral2 works fine for smaller and larger numbers

eps0 = 8.8541878128*1e-12;
f = @(x,z) (1/(4*pi*eps0))*(1./sqrt(x.^2+z.^2));
%rand returns a value between 0 and 1
a = rand;
b = rand;
integral2(f, -a/2, a/2, -b/2, b/2)

ans = 1.3207e+10

a = 1e5*rand;
b = 1e5*rand;
integral2(f, -a/2, a/2, -b/2, b/2)

ans = 8.0193e+14

[Matt J]

Well, there's no real reason to be pursuing integral2 anymore, now that we have an explicit closed-form formula for the interval

[John D'Errico]

Easier than that even.

Answer 1

[Ayush]

I understand that you are unable to perform numerical integration using quad2d and integral2 function. According to MathWorks documentation, they perform best when the singularities are on the integration boundary. So, as the function to be evaluated is symmetric in all quadrants you can change the limits such that point of singularity lies on the boundary of the integration and multiply the results accordingly. Follow the documentation links here:

1. Numerically evaluate double integral - MATLAB integral2 - MathWorks India
2. Numerically evaluate double integral — tiled method - MATLAB quad2d - MathWorks India

Answer 2

[John D'Errico]

As is often the case, this is FAR simpler than you think. And, yeah, sometimes things are far more difficult. They never seem to lie right in the middle though. ;-)

Is your function symmetrical, in the sense that you can split it into 4 pieces, each of which is identically the same? YES!

So break the problem into 4 pieces, integating from 0 to bx/2 in x, and 0 to by/2 in y. Then multiply by 4.

syms x y real positive
syms bx by real positive
syms eps0

eps0 is irrelevant, since it is on the outside, just a scalar multiplier. So pull it out of the integral.

Ival = 4/(4*sym(pi)*eps0)*int(int(1/sqrt(x^2 + y^2),x,[0,bx/2]),y,[0,by/2])

Ival = 

And you should see the problem is now trivial to solve. The complicated form in the comments above probably reduces to this.

As well, even integral2 will now have no problems, IF you had numerical values. For example, I'll try it out, on a simple example. Remember that integral2 is happier if the singularity is in a corner.

subs(Ival,[bx,by],[2 3])

ans = 

vpa(ans)

ans = 

and now use integral 2. Will it get it right? Just assume eps0 is 1, since it is irrelevant anyway.

format long g
integral2(@(X,Y) 1./sqrt(X.^2 + Y.^2),0,2/2,0,3/2)/pi

ans = 
         0.678789758361919

So if I just disregard the value of eps0, integral2 agrees, out to about 8 digits. We probably can't do much better in double precision anyway.

Finally, I would point out that literally whenever a singularity resides inside an interval, you are probably best off finding a way to reduce the problem to one where the singularity lives at a boundary. Also, even if this problem had not had symmetric limits of integration, you ould still have reduced this to four separate integrals, all of the form I generated.

Oh, one other point. While you say that converting to polar coordinates would be difficult, that too is far easier than you think. Not worth doing the work here, since I've already given a complete solution.