# Plotting / Animating a test planet around a star

I am trying to plot/animate the motion of a test planet around a star using Mathematica in the framework of general relativity. In fact, I want to see the perihelion shift. I am using as inspiration some of the code from Mathematica for theoretical physicists, although my code is a bit different.

I have managed to compute the Euler-Lagrange equations from the Schwarzschild Lagrangian using a package called variational methods.

I am not sure how to proceed to plot/animate the motion of a test planet around a star from here.

Below is the code that gets the Euler-Lagrange equations. How do I use this to plot and animate the orbit? So basically what lines of code do I add for this?

ds2 = c^2 *(1 - Rs/r) Dt[t]^2 - 1/(1 - Rs/r)*Dt[r]^2 -
r^2  *Dt[\[Phi]]^2

schwarzschildLagrangian =
c^2*(1 - Rs/r[s])*(D[t[s], s])^2 - (D[r[s], s])^2 /(1 - Rs/r[s]) -
1/2 * r[s]^2 *(D[\[Phi][s], s])^2

Needs["VariationalMethods"]

swEquations =
EulerEquations[schwarzschildLagrangian, {r[s], \[Phi][s], t[s]},
s];

swEquations // Expand // TableForm

angularMomentum =
Map[Integrate[#, s] &, {swEquations[[2, 1]]}][[1]] == 1/Sqrt[\
[Beta]]

energy =
MapAt[Integrate[#, s] &, {swEquations[[3, 1]]}, 1][[1]] == -
k^2 /(c^2 * Sqrt[\[Beta]])


Where r is the radial distance, t is the time, Rs is the Schwarzschild radius = $$\frac{2GM}{c^2}$$, and $$\phi$$ is the angle. (s) means they depend on s. Ok, I am trying NDSolve now

I have now added more code to NDSolve the equations of motion, and am now working on code to parametrically plot things.

(*Specify initial conditions*)
r0 = 1 + \[CurlyEpsilon]; (*Initial radial position with offset*)
\[Theta]0 = \[Pi]/4; (*Initial polar angle*)
\[Phi]0 = 0; (*Initial azimuthal angle*)
t0 = 0; (*Initial time*)
v0 = 0; (*Initial radial velocity*)
w0 = 1; (*Initial angular velocity*)
u0 = 1; (*Initial time velocity*)

(*Specify numerical parameters*)
c = 1; (*Speed of light*)
\[CurlyEpsilon] = 1*^-6; (*Small offset to avoid singularity*)
sRange = {0, 10};

NDSolve[{swEquations[[1]], swEquations[[2]], swEquations[[3]],
r[0] == r0, \[Phi][0] == \[Phi]0, t[0] == t0,
r'[0] == v0, \[Phi]'[0] == w0, t'[0] == u0}, {r[s], \[Phi][s],
t[s]}, {s, 1, 10}, Method -> "StiffnessSwitching"]

• Below is the code that gets the Euler-Lagrange equations. it helps if you explain to the reader what $r(s),Rs,t(s),s,\phi(s)$ are meant to be physically instead of just giving code. After all, to animate this, one needs to know what these represent. A picture/diagram will also help. Look at Plot3D and Manipulate and decide what sliders you need (i.e. what are your control variables). Ofcourse you need to first use NDSolve to solve the odes's. So need initial conditions as well. I would suggest you try to first solve the three equations of motion you show numerically first before animation. May 21, 2023 at 21:31
• r is the radial distance, t is the time, Rs is the Schwarzschild radius = 2GM/c^2, phi is the angle. (s) means they depend on s. Ok, I am trying NDSolve now. May 21, 2023 at 21:37
• I do not understand. Normally equations of motion depends on time. i.e. position is function of time. But here you have position depending on s and time also depends on s? So to move the planet, are you not supposed to advance time? Now you need to advance s. Where does does time come into play? What will the solution $t(s)$ be used for? Since to move the planet all what you need is $r(s),\phi(s)$ and changing s. May 22, 2023 at 1:02
• In general relativity, you're free to choose arbitrary coordinate frames, and time is regarded as one of the "coordinates'. So there's no single definition as there is in non-relativistic physics. I believe he's treating all four coordinates parametrically, using 's' as his parameter. For example, t might be coordinate time in the Schwartzchild frame, and 's' might represent proper time along a world line. So using s as "time", and solving for t[s] is reasonable for this case.
– user87932
May 22, 2023 at 2:12
• @kevinTahN. Did you seen this post mathematica.stackexchange.com/questions/189961/… ? May 22, 2023 at 2:12

Here is basic animation using only $$r(s),\phi(s)$$ in 2D.

Animate[frames[[n]], {n, 1, Length@pts, 1}]


## Code

ds2 = c^2*(1 - Rs/r) Dt[t]^2 - 1/(1 - Rs/r)*Dt[r]^2 -
r^2*Dt[ϕ]^2
schwarzschildLagrangian =
c^2*(1 - Rs/r[s])*(D[t[s], s])^2 - (D[r[s], s])^2/(1 - Rs/r[s]) -
1/2*r[s]^2*(D[ϕ[s], s])^2
Needs["VariationalMethods"]
swEquations =
EulerEquations[schwarzschildLagrangian, {r[s], ϕ[s], t[s]},
s];
(*Specify initial conditions*)
r0 = 1 + ε; (*Initial radial position with offset*)
θ0 = π/4; (*Initial polar angle*)
ϕ0 = 0; (*Initial azimuthal angle*)
t0 = 0; (*Initial time*)
v0 = 0; (*Initial radial velocity*)
w0 = 1; (*Initial angular velocity*)
u0 = 1; (*Initial time velocity*)

(*Specify numerical parameters*)
c = 1; (*Speed of light*)
ε = 1*^-6; (*Small offset to avoid singularity*)
maxS = 40;
sol = First@
NDSolve[{swEquations[[1]], swEquations[[2]], swEquations[[3]],
r[0] == r0, ϕ[0] == ϕ0, t[0] == t0,
r'[0] == v0, ϕ'[0] == w0, t'[0] == u0},
{r, ϕ, t},
{s, 1, maxS}, Method -> "StiffnessSwitching"];

pts = Table[{Evaluate[r[s0] /. sol]*Cos[Evaluate[ϕ[s0] /. sol]],
Evaluate[r[s0] /. sol]*Sin[Evaluate[ϕ[s0] /. sol]]}, {s0, 1,
maxS, .1}];
time = Table[Evaluate[t[s0] /. sol], {s0, 1, maxS, .1}];
frames =
Table[Grid[{{Row[{"time = ", time[[n]]}]}, {ListLinePlot[
pts[[1 ;; n]], PlotRange -> {{-3, 3}, {-3, 3}},
ImageSize -> 300]}}, Frame -> All, Spacings -> {1, 1}], {n, 1,
Length@pts}];
Animate[frames[[n]], {n, 1, Length@pts, 1}]