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

Improve this question
$\endgroup$
8
  • 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
    Commented Jul 15, 2018 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
    Commented Jul 15, 2018 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$
    – diffusiondiver11
    Commented Jul 15, 2018 at 19:55
  • $\begingroup$ @MariuszIwaniuk it sure may help. Thanks! $\endgroup$
    – diffusiondiver11
    Commented Jul 15, 2018 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
    Commented Jul 15, 2018 at 20:19

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
5
  • $\begingroup$ How much time did it take your system to compute int1? $\endgroup$
    – diffusiondiver11
    Commented Jul 16, 2018 at 11:09
  • $\begingroup$ It takes 16 seconds. $\endgroup$
    – Akku14
    Commented Jul 16, 2018 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$
    – diffusiondiver11
    Commented Jul 16, 2018 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
    Commented Jul 16, 2018 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$
    – diffusiondiver11
    Commented Jul 16, 2018 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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