# Gravitational potential of a cube

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?

• source of the mysterious equation here – andre314 Mar 3 '18 at 8:39
• Similar question was already asked here mathematica.stackexchange.com/questions/160755/… – yarchik Mar 3 '18 at 12:37
• the interesting aspect of OP's question is : how to retranscript a formula found in a book in Mathematica code, in a way that respects the formulation of the book ? The idea is to do a minimal effort retranscription too. – andre314 Mar 3 '18 at 23:01
• "every little xi must vanish after the last step" = mysterious text – andre314 Mar 3 '18 at 23:06
• What I mean is, that the x_i are only some variables for the integration, but the Potential, after all calculations are made, should only contian: G, rho, X_i and D_i. – M. K. Mar 4 '18 at 18:12

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},

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: 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 x x, r -> Sqrt[x^2 + x^2 + x^2]};

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

potentialV[X1_,X2_,X3_]=Fold[
(
(#1 /. #2 /. sign -> 1)
-(#1 /. #2 /. sign -> -1)
)&,
sumOverI,
{
x-> sign D1-X1,
x-> sign D2-X2,
x-> 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

• Thank you so much. Is there some other term to discribe the outer summation then, since its not a normal summation? I want to understand whats going on here, but I need to know what I am looking for. – M. K. Mar 4 '18 at 18:18
• @MichelK. It's more like a composition of functions than a summation. A big difference is that a composition leads to a exponential increase of the number of terms (composition x->f[x] - g[x] n times leads to 2^n terms. In our case f and g are replacement of some terms x_i). Hope I'm relatively clearer – andre314 Mar 4 '18 at 18:31
• After thinking for quite a while, I still dont get it, sorry. But thanks for your time. Where is the composition? It would help me, if someone has the time to explain the steps/the scheme this Equation wants me to do, if I would calculate this by hand. Also this equation produces 48=3*2^4 Terms. :( – M. K. Mar 8 '18 at 16:48

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, X, X] == -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

• Can you add into yours answer a gravitational potential plot? – Mariusz Iwaniuk Mar 3 '18 at 12:15
• @MariuszIwaniuk I don't think I can, because as I explain in the answer, the (implementation of the) expression for V contains free parameters x[j] (lower-case x's) – user42582 Mar 3 '18 at 13:12
• So this means this formula is wrong? This Potential doesnt make any sense, if there are free parameters besides X[i]. There must be some other approach, which satisfies my assertion, right? But anyways thank you for your contribution. – M. K. Mar 3 '18 at 15:45
• @MichelK. you're welcome; I think the notation is misleading or the function of the brackets is different from the the usual – user42582 Mar 3 '18 at 16:21