2
$\begingroup$

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*)
Rs = 1; (*Schwarzschild radius*)
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"]
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    May 21, 2023 at 21:31
  • $\begingroup$ 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. $\endgroup$ May 21, 2023 at 21:37
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    May 22, 2023 at 1:02
  • 1
    $\begingroup$ 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. $\endgroup$
    – user87932
    May 22, 2023 at 2:12
  • 1
    $\begingroup$ @kevinTahN. Did you seen this post mathematica.stackexchange.com/questions/189961/… ? $\endgroup$ May 22, 2023 at 2:12

1 Answer 1

3
$\begingroup$

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

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

enter image description here

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*)
Rs = 1; (*Schwarzschild radius*)
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}]
$\endgroup$

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.