Gravitational potential of a cube

Question

I want to write code that calculates the gravitational potential of an arbitrarily shaped celestial body. To understand the calculation, I started with an easy shape: a rectangle (or a cube). I realised this was allready hard enough, so I researched in the Internet and found a work of this problem. I found the equation:

Where v is

v = Subscript[x, 1]*Subscript[x, 2]*Subscript[x, 3],

and the $D_i$ are the lenghs of the edges of the cube.

I have 2 questions (but the first one will be answered too if the second one gets answered):

I try to understand how this sum looks like, but I dont quite get how this summation over the Integration Limits works. I guess every little $x_i$ must vanish after the last step.

So my actual question is:

How can I code these integration limit brackets and sum over them?

Answer 1

This is Not a answer of yours Question its only alternative solution.

G = 1;
ρ = 1;
a = 1;
b = 1;
V[X_, Y_] := -G*ρ*NIntegrate[1/Sqrt[(X - x)^2 + (Y - y)^2], {x, -a, a}, {y, -b, b}, 
Method -> "LocalAdaptive"] // Quiet

n = 1/15;
ListContourPlot[Partition[Flatten[Table[{x, y, V[x, y]}, {x, -2, 2, n}, {y, -2, 2, n}]], 3], 
Contours -> 40] 

Gravitational potential of a square:

Answer 2

I find that the notations used in the formula are misleading.

• the inner summation is a normal summation

• The outer "summation" describes in fact a list of nested (hence the use of the instruction Fold[] in the code below) replacement rules associated with summations.

Here is the transcription in Wolfram Language (= Mathematica) :

sumOverIBaseElt[i_] = v/x[i] Log[x[i]+r] - (x[i]^2)/2 ArcTan[v/(r*x[i]^2)] /. 
{v -> x[1] x[2] x[3], r -> Sqrt[x[1]^2 + x[2]^2 + x[3]^2]};

sumOverI=Sum[sumOverIBaseElt[i],{i,1,3}];

potentialV[X1_,X2_,X3_]=Fold[
(
   (#1 /. #2 /. sign -> 1)
   -(#1 /. #2 /. sign -> -1)
)&,
sumOverI,
{
x[1]-> sign D1-X1,
x[2]-> sign D2-X2,
x[3]-> sign D3-X3
}];

ContourPlot[
      Evaluate[potentialV[x,y,0] /. {D1->1,D2->1,D3->3}],
      {x,-2,2},{y,-2,2},
      PlotLegends->Automatic,
      Exclusions -> None,
      Epilog -> {Dashed,Line[{{-1,-1},{-1,1},{1,1},{1,-1},{-1,-1}}]}]

Calculate the gravitational field vectors.

StreamPlot[
Evaluate[{
    D[potentialV[x, y, 0] /. {D1->1,D2->1,D3->3}, x],
    D[potentialV[x, y, 0] /. {D1->1,D2->1,D3->3}, y]}],
{x, -2, 2},{y, -2, 2}, 
StreamPoints -> Fine]

EDIT version 2

also possible, and finally simpler :

potentialV[X1_,X2_,X3_]=
    (v/xi Log[xi+r] - (xi^2)/2 ArcTan[v/(r*xi^2)] /. 
    {v->x1 x2 x3,r->Sqrt[x1^2 + x2^2 + x3^2]}) /. 
    {{xi-> x1},{xi-> x2},{xi-> x3}} //
    (Plus @@ #)& //
    ((# /. x1->D1-X1)-(# /. x1->-D1-X1))& //
    ((# /. x2->D2-X2)-(# /. x2->-D2-X2))& //
    ((# /. x3->D3-X3)-(# /. x3->-D3-X3))&;  

-> same result

Answer 3

If my interpretation of how the square brackets are used (as is the case in definite integration) is correct, then the following code calculates $V(X_1,X_2,X_3)$.

(* localize symbols *)
Block[{x, r, Dd, X, G, ρ},
 (* expression inside the Sum inside the square brackets *)
 With[{v = Times @@ Array[x, 3]},
  f[i_] := (v/x[i]) Log[x[i] + r] - ((x[i]^2)/2) ArcTan[v/(x[i]^2 r)]
  ];
 (* the expression inside the square brackets *)
 With[{sum = Plus @@ Array[f, 3]},
  g[j_] := Plus @@ (
     {1, -1} (sum /. {
         {x[j] -> Dd[j] - X[j]},
         {x[j] -> -Dd[j] - X[j]}
         }
       )
     )
  ];
 (* actual definition of V *)
 V[args__] := With[{params = Array[X, 3]},
   With[{sum = Plus @@ Array[g, 3]},
    -G ρ sum /. Thread[params -> {args}]
    ]
   ]
 ]

Evaluating V[a, b, c] returns a lengthy expression (which I won't reproduce here). It seems, though, that the assertion in the question

"[...] I guess every little xi must vanish after the last step"

is not confirmed.

---

End notes

Verify that indeed V performs the calculations in the right order:

Evaluating

V[X[1], X[2], X[3]] == -G ρ Plus @@ (
  (
    (* perform the operation implied by the square brackets *)
    Plus @@ ({1, -1} (
      With[{v = Times @@ Array[x, 3]},
        (* expression inside square brackets *)
        Plus @@ Array[(v/x[#]) Log[x[#] + r] - ((x[#]^2)/2) ArcTan[v/(x[#]^2 r)] &, 3]
       ] /. #))
           (* construct the limits of the square brackets *)
   ) & /@  Array[{{x[#] -> Dd[#] - X[#]}, {x[#] -> -Dd[#] - X[#]}} &, 3]
 )

returns

True