2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ xyz is not defined. $\endgroup$ – Anton Antonov Aug 16 '18 at 16:23
  • $\begingroup$ In the formula, integration is carried out with respect to x', whereas in the code with respect to x. How really? $\endgroup$ – Alex Trounev Aug 16 '18 at 16:29
  • $\begingroup$ 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 $\endgroup$ – Jan SE Aug 16 '18 at 17:14
  • $\begingroup$ In the Mathematica code the integration is also carried out with respect to x'. x' is indicated as xs in the code $\endgroup$ – Jan SE Aug 16 '18 at 17:16
  • $\begingroup$ 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. $\endgroup$ – Bill Watts Aug 16 '18 at 18:35
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.