2
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Maybe this helps: mathematica.stackexchange.com/questions/173253/…. Try user Henrik Schumacher answer. You must only to modify the code. $\endgroup$ – Mariusz Iwaniuk Jul 15 '18 at 19:21
  • $\begingroup$ 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... $\endgroup$ – Henrik Schumacher Jul 15 '18 at 19:46
  • $\begingroup$ 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. $\endgroup$ – シャシュワト Jul 15 '18 at 19:55
  • $\begingroup$ @MariuszIwaniuk it sure may help. Thanks! $\endgroup$ – シャシュワト Jul 15 '18 at 19:57
  • 1
    $\begingroup$ 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). $\endgroup$ – David G. Stork Jul 15 '18 at 20:19
6
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ How much time did it take your system to compute int1? $\endgroup$ – シャシュワト Jul 16 '18 at 11:09
  • $\begingroup$ It takes 16 seconds. $\endgroup$ – Akku14 Jul 16 '18 at 13:17
  • $\begingroup$ 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? $\endgroup$ – シャシュワト Jul 16 '18 at 13:20
  • $\begingroup$ 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 . $\endgroup$ – Akku14 Jul 16 '18 at 13:33
  • $\begingroup$ I see. I'm on 11.3 and 11.2. It could be the issue here. Algorithms might be redefined. $\endgroup$ – シャシュワト Jul 16 '18 at 13:42

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.