0
$\begingroup$

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:

enter image description here

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?

$\endgroup$
  • 2
    $\begingroup$ source of the mysterious equation here $\endgroup$ – andre314 Mar 3 '18 at 8:39
  • $\begingroup$ Similar question was already asked here mathematica.stackexchange.com/questions/160755/… $\endgroup$ – yarchik Mar 3 '18 at 12:37
  • $\begingroup$ 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. $\endgroup$ – andre314 Mar 3 '18 at 23:01
  • $\begingroup$ "every little xi must vanish after the last step" = mysterious text $\endgroup$ – andre314 Mar 3 '18 at 23:06
  • $\begingroup$ 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. $\endgroup$ – M. K. Mar 4 '18 at 18:12
2
$\begingroup$

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:

enter image description here

$\endgroup$
2
$\begingroup$

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

enter image description here

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]

enter image description here

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – M. K. Mar 4 '18 at 18:18
  • $\begingroup$ @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 $\endgroup$ – andre314 Mar 4 '18 at 18:31
  • $\begingroup$ 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. :( $\endgroup$ – M. K. Mar 8 '18 at 16:48
0
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Can you add into yours answer a gravitational potential plot? $\endgroup$ – Mariusz Iwaniuk Mar 3 '18 at 12:15
  • $\begingroup$ @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) $\endgroup$ – user42582 Mar 3 '18 at 13:12
  • $\begingroup$ 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. $\endgroup$ – M. K. Mar 3 '18 at 15:45
  • $\begingroup$ @MichelK. you're welcome; I think the notation is misleading or the function of the brackets is different from the the usual $\endgroup$ – user42582 Mar 3 '18 at 16:21

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.