How to increase performance of this code for plotting a contour plot?

Question

I have an equation which I need to triple integrate over a unit cube. The equation is

pot = NIntegrate[1/Sqrt[(x - h)^2 + (y - k)^2 + (z - l)^2], {h,-1,1}, {k,-1, 
1},{l,-1,1}];

As soon as I enter Shift+Enter it immediately processes the command. But now what I want is to plot its ContourPlot for different ${z}$ values (I chose $z=0.5$). So I give the command

ContourPlot[pot /. {z -> 0.5}, {x, -2, 2}, {y, -2, 2}]

But this piece of code takes just forever to process. I just keep on waiting and waiting but processing never ends (it takes really really long time). I am not sure that how is this such a computationally heavy task. For $z$ other than $0$ it takes longer time.

Is there something that I am doing wrong? I don't think this is a drawback of the device I am using. Is there a way to improve the performance of this code I am using?

P.S. It's been more than 10 minutes but the code for $z=0.5$ has not processed.

For your reference, I am attaching the contour plot for $z=0$.

This is the output for $z=0.5$ from the code above (it took about 10 minutes)

Maybe this helps: https://mathematica.stackexchange.com/questions/173253/mathematica-taking-forever-to-compute/173259#173259. Try user Henrik Schumacher answer. You must only to modify the code. — Mariusz Iwaniuk, Jul 15 '18 at 19:21
The problem with these integrals is that they are singular and three-dimensional. That makes it quite expensive. Unfortunately, one also cannot exploit symmetry (e.g. by polar coordinates) in an obvious way... — Henrik Schumacher, Jul 15 '18 at 19:46
But even Vector3D plots are taking such a long time that I could not output them. I did the same with other 3 dimensional Integrals but they took at most a few seconds. This one's not. — diffusiondiver11, Jul 15 '18 at 19:55
pot has no definition of the variables within it, so of course it will never be evaluated for z = 0.5 (or indeed any other value). — David G. Stork, Jul 15 '18 at 20:19
It did was evaluated for both z=0 and z=0.5 . But it took really long time to to that. 10-15 minutes as I remember. The picture is the output for the same code I gave above in description for z=0. You can run the code if you want, it will give the same output. — diffusiondiver11, Jul 15 '18 at 20:22
The reason you got the same picture "for different values of $z$" is that pot is independent of z, so of course you're merely plotting the same function. — David G. Stork, Jul 15 '18 at 20:25
No. I got different outputs for both values. The code really works but is very inefficient. I have added the other output to the description of the problem. pot is not independent of $z$. Look the SquareRoot term has $(z-l)^2$ term. Infact pot is a function of $x$, $y$ and $z$. $pot(x,y,z)$ — diffusiondiver11, Jul 15 '18 at 20:28

score 6 · Accepted Answer · edited Jul 16 '18 at 15:21

One helpfull rule to get fast integration is, to do analytical integration as much as you can.

int1 = Integrate[1/Sqrt[(x - h)^2 + (y - k)^2 + (z - l)^2], {l, -1, 1}, 
         Assumptions -> -1 <= h <= 1 && -1 <= k <= 1 && x \[Element] Reals &&
                        y \[Element] Reals && z \[Element] Reals]

(* -Log[-z + Sqrt[h^2 + k^2 - 2 h x + x^2 - 2 k y + y^2 + z^2]] - 
    Log[z + Sqrt[h^2 + k^2 - 2 h x + x^2 - 2 k y + y^2 + z^2]] + 
    Log[1 - z + Sqrt[1 + h^2 + k^2 - 2 h x + x^2 - 2 k y + y^2 - 2 z + z^2]] + 
    Log[1 + z + Sqrt[1 + h^2 + k^2 - 2 h x + x^2 - 2 k y + y^2 + 2 z + z^2]] *)

I do the second integration with the rule based integrator (Rubi) by Albert Rich (see http://www.apmaths.uwo.ca/~arich/ ), because, in contrast to Mathematica, it gives an antiderivative without discontinuities.

rint2[x_, y_, z_, h_, k_] = Int[int1, k];

Take integration values at borders to get the definite integral.

rint2def[x_, y_, z_, h_] = 
     rint2[x, y, z, h, 1] - rint2[x, y, z, h, -1] // 
      Simplify[#, Assumptions -> -1 <= h <= 1 && -1 <= k <= 1 && 
         x \[Element] Reals && y \[Element] Reals && z \[Element] Reals] &

(*   -h ArcTan[((-1 + y) (-1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 - 2 z + z^2])] + 
x ArcTan[((-1 + y) (-1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 - 2 z + z^2])] + 
h ArcTan[((1 + y) (-1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2])] - 
x ArcTan[((1 + y) (-1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2])] - 
h ArcTan[((1 - y) (1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2])] + 
x ArcTan[((1 - y) (1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2])] - 
h ArcTan[((1 + y) (1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2])] + 
x ArcTan[((1 + y) (1 + z))/((h - x) Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2])] - (-1 + 
z) ArcTanh[(1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 - 2 z + z^2]] - 
y ArcTanh[Sqrt[2 + h^2 - 2 h x + x^2 - 2 y + y^2 - 2 z + z^2]/(
1 - z)] - 
ArcTanh[(-1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2]] + 
z ArcTanh[(-1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2]] + 
y ArcTanh[Sqrt[2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2]/(
1 - z)] + 
ArcTanh[(1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2]] + 
z ArcTanh[(1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2]] - 
y ArcTanh[Sqrt[2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2]/(
1 + z)] - 
ArcTanh[(-1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2]] - 
z ArcTanh[(-1 - y)/Sqrt[
2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2]] + 
y ArcTanh[Sqrt[2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2]/(
1 + z)] - Log[-z + Sqrt[1 + h^2 - 2 h x + x^2 - 2 y + y^2 + z^2]] -
Log[-z + Sqrt[1 + h^2 - 2 h x + x^2 + 2 y + y^2 + z^2]] + 
Log[1 - z + Sqrt[2 + h^2 - 2 h x + x^2 + 2 y + y^2 - 2 z + z^2]] + 
Log[((1 - z + Sqrt[
   2 + h^2 - 2 h x + x^2 - 2 y + y^2 - 2 z + z^2]) (1 + z + Sqrt[
   2 + h^2 - 2 h x + x^2 - 2 y + y^2 + 2 z + z^2]) (1 + z + Sqrt[
   2 + h^2 - 2 h x + x^2 + 2 y + y^2 + 2 z + z^2]))/((z + Sqrt[
   1 + h^2 - 2 h x + x^2 - 2 y + y^2 + z^2]) (z + Sqrt[
   1 + h^2 - 2 h x + x^2 + 2 y + y^2 + z^2]))]   *)

The last integration has to be done numericaly.

rint3[x_, y_, z_] := NIntegrate[rint2def[x, y, z, h], {h, -1, 1}]

ContourPlot now finishes within 21 seconds.

ContourPlot[rint3[x, y, 1/2], {x, -2, 2}, {y, -2, 2}, 
             ImageSize -> 400] // Timing

Strange that it takes very long for me. I had to break it into pieces that way it was fast. Like you did in rint2def. Any idea why could it be? — diffusiondiver11, Jul 16 '18 at 13:20
I am working with MMA Version 8.0. May be it works with other assumptions, which yield the same result int1 = Integrate[ 1/Sqrt[(x - h)^2 + (y - k)^2 + (z - l)^2], {l, -1, 1}, Assumptions -> -1 < h < 1 && -1 < k < 1 && x > 0 && y > 0 && z > 0] // Timing . — Akku14, Jul 16 '18 at 13:33
I see. I'm on 11.3 and 11.2. It could be the issue here. Algorithms might be redefined. — diffusiondiver11, Jul 16 '18 at 13:42

How to increase performance of this code for plotting a contour plot?

1 Answers1

Linked