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

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: mathematica.stackexchange.com/questions/173253/…. Try user Henrik Schumacher answer. You must only to modify the code. Jul 15, 2018 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... Jul 15, 2018 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. Jul 15, 2018 at 19:55
• @MariuszIwaniuk it sure may help. Thanks! Jul 15, 2018 at 19:57
• 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). Jul 15, 2018 at 20:19

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

• How much time did it take your system to compute int1? Jul 16, 2018 at 11:09
• It takes 16 seconds. Jul 16, 2018 at 13:17
• 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? Jul 16, 2018 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 . Jul 16, 2018 at 13:33
• I see. I'm on 11.3 and 11.2. It could be the issue here. Algorithms might be redefined. Jul 16, 2018 at 13:42