Numerical solution of 3 dim integral with singularity

Question

I want to solve the following integral numerically with Mathematica:

$\int_{0}^{L_x}\int_{0}^{L_y}\int_{0}^{L_z}d^3x'\frac{\sin^2(\frac{x'\pi}{L_x})\sin^2(\frac{y'\pi}{L_y})\sin^2(\frac{z'\pi}{L_z})}{|x-x'|}e^{i\omega |x-x'|}$

But I am not getting the result which I expect. Therefore I have the feeling that Mathematica cannot deal with the singularity of the integrand ($x\in [0,L_x]\times[0,L_y]\times[0,L_z]$). By the way: Does anyone know how to prove if the integral is finite?

I already tried different Methods/MinRecursion/MaxRecursion/MaxPoints in NIntegrate which more or less all yield the same result.

Here is my code:

Lx = 0.09;  
Ly = 0.03;  
Lz = 0.4;  
f = 10*10^9;  
c = 3*10^8;  
omega = 2*Pi*f/c;  
x = 0.045
y = 0.015
z = 0.15

B0 = Sin[Pi*xs/Lx]^2*Sin[Pi*ys/Ly]^2*Sin[Pi*zs/Lz]^2;  
Integrand = B0/((x-xs)^2+(y-ys)^2+(z-zs)^2)^(1/2)*Exp[I*omega*((x-xs)^2 +(y-ys)^2 +(z-zs)^2)^(1/2)];  


NIntegrate[Integrand, {xs, 0, Lx}, {ys, 0, Ly}, {zs, 0, Lz}];  

The thing which worries me a little is that when I use Lx=0.09, Ly=0.03, Lz=0.4 and evaluate the integral at (x=(0.045,0.015,0.15)m) I get -0.00508-i 0.000606. But when I evaluate it at (x=(0.044,0.014,0.15)) I get -0.00024+0.000198. So there seems to be a discontinuity at (x=(0.045,0.015,0.15)m). All other points ($x\in [0,L_x]\times[0,L_y]\times[0,L_z]$) yield a much smaller results than the point at (x=(0.045,0.015,0.15)m). When I plot the results in the xy plane I have to exclude this point such that one can see the structure...

All comments are welcomed.

Many thanks in advance.

Answer 1

You could try a different rule/method and indicate the singularity:

First integral:

On@Assert;
Block[{x = 0.045, y = 0.015, z = 0.15},
 Assert[0 <= x <= Lx && 0 <= y <= Ly && 0 <= z <= Lz]; (* warn if outside domain *)
 NIntegrate[integrand, {xs, 0, x, Lx}, {ys, 0, y, Ly}, {zs, 0, z, Lz},
   Method -> "CartesianRule"]
 ]
Off@Assert;

(*  -0.00024633 + 0.000202509 I  *)

2nd integral:

On@Assert;
Block[{x = 0.044, y = 0.014, z = 0.15},
 Assert[0 <= x <= Lx && 0 <= y <= Ly && 0 <= z <= Lz];
 NIntegrate[integrand, {xs, 0, x, Lx}, {ys, 0, y, Ly}, {zs, 0, z, Lz},
   Method -> "CartesianRule"]
 ]
Off@Assert;

(*  -0.000246293 + 0.000198098 I  *)

They seem reasonably close.

Note the original answer used the following, which emits a slow convergence warning but gives the same result to ten nonzero digits and is much faster:

Block[{x = 0.044, y = 0.014, z = 0.15},
 NIntegrate[integrand, {xs, 0, Lx}, {ys, 0, Ly}, {zs, 0, Lz}, 
  Method -> "CartesianRule"]
 ]

The NIntegrate::slwcon slow convergence warning does not in itself indicate an error. If the integral had failed to converge, it can help diagnose the problem. However, for the other parameter setting x = 0.045, y = 0.015, z = 0.15, the integral had a small but significant error and had only five digits of precision, which is not horrible.