# Numerical solution of 3 dim integral with singularity

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...

• xyz is not defined. Commented Aug 16, 2018 at 16:23
• In the formula, integration is carried out with respect to x', whereas in the code with respect to x. How really? Commented Aug 16, 2018 at 16:29
• I edited my post and it should work now. You can play with the x and y values to see the feature which I mentioned in the post Commented Aug 16, 2018 at 17:14
• In the Mathematica code the integration is also carried out with respect to x'. x' is indicated as xs in the code Commented Aug 16, 2018 at 17:16
• I think the problem may be that at x = x', Sin[0]/0 converges while Cos[0]/0 does not. Try rationalizing all your numbers, increasing your WorkingPrecision and play with Sin[omega ...] instead of Exp[I*omega ...]. When I do that I get closer numbers for the two data sets than you. Commented Aug 16, 2018 at 18:35

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.