# How to directly input calculated Geodesic differential equation into *NDSolve*

I am trying to find the trajectory of a test particle around a relativistic object. Initially, I was trying the Kerr metric as spacetime but the final simulations clearly were not correct as I saw the particle being repelled by the star/planet/black-hole which is not realistic. Hence there must be something wrong going on with the equations I have calculated and used. Since the equations were humongous I was not able to rectify the mistake (although it was not me but my computer that calculated the differential equations, so the equations would be wrongs seems not so plausible to me). I did modify a little bit of the equation I was calculated. I'll show how. (I'll use the Schwarzschild metric and not Kerr metric as the equations are very ugly)

First I defined the metric tensor for the schwarzchild geometry.

schwar = {{-(1 - (2 G M)/(r c^2))^-1, 0, 0, 0}, {0, -r^2, 0, 0}, {0,
0, -r^2 Sin[\[Theta]]^2, 0}, {0, 0, 0, (1 - (2 G M )/(c^2 r))}};


then I defined the inverse metric

inverseschwar = Inverse[schwar];


Now, I have to calculate the Christoffel symbols for the metric. I defined a function christ for this

christ[ a_, b_, c_] := christ[a, b, c] =   Sum[1/2 inverseschwar[[a, d]]
(D[schwar[[d, c]], coord[[b]]] + D[schwar[[d, b]],
coord[[c]]] - D[schwar[[b, c]], coord[[d]]]),
{d, 1, 4}]


Now I have to compute the geodesic differential equations from the geodesic equation

Manipulate[coord[[i]]'' + Sum[christs[i, h, k] coord[[h]]'
coord[[k]]', {h, 1, 4}, {k, 1, 4}], {i, {1, 2, 3, 4}}]


This above are the 4 differential equations for each coordinate $r$,$\theta$,$\phi$ and $t$. Notice that I have purposely put on the differentiation operators $''$ and $'$ . But these operators don't do differentiation in mathematica (Once you run this code you'll understand why. They come on on the variable making in into a new variable).

For example: output for the equation of $\theta$ coordinate which is

Here the $\theta$ and $\theta^{'}$ are being considered as two different coordinates and not as differentiation (w.r.t proper time $s$) of $\theta$. Moreover what I need is that all coordinates to be labeled as a function of $s$ so that integration process can be carried out. e.g. $\theta$ should be represented as $\theta[s]$. I cannot do this by just putting an $s$ in front of every variable before calculating these equations. I can only do it for derivative terms present in the equations but cannot for the terms that arise from the chirstoffel symbols as if suppose I put the argument $s$ before calculating the chirstoffel symbols then as the metric is being differentiated w.r.t coordinates $r$, $\theta$, $\phi$ and $t$ so every chirstoffel symbol will become zero as $s$ is independent of all coordinates.

So now I have to manually modify the equationsto make it look like this

After making it look like above now I can feed this and other $3$ equation to NDSolve function and get the output trajectory of the particle.

So, my question is, is there a way to directly have mathematica produce the equation that can directly be inserted into NDSolve function? I don't want to make the process manual at any step.

I am listing the 4 differential equations for your reference

(G M r'[s]^2)/(2 G M r[s] - c^2 r[s]^2) + (
G M (-2 G M + c^2 r[s]) t'[s]^2)/(
c^4 r[s]^3) + ((2 G M)/c^2 - r[s]) \[Theta]'[
s]^2 + ((2 G M - c^2 r[s]) Sin[\[Theta][s]]^2 \[Phi]'[s]^2)/c^2 +
r''[s]==0;


Theta Equation

(2 r'[s] \[Theta]'[s])/r[s] -
Cos[\[Theta][s]] Sin[\[Theta][s]] \[Phi]'[s]^2 + \[Theta]''[s]==0;


Phi Equation

(2 r'[s] \[Phi]'[s])/r[s] +
2 Cot[\[Theta][s]] \[Theta]'[s] \[Phi]'[s] + \[Phi]''[s]==0;


Time Equation

(2 G M r'[s] t'[s])/(r[s] (-2 G M + c^2 r[s])) + t''[s]==0;

• In mathematica you should write e.g. D[t[s],s] instead of t'[s] and D[t[s],s,s] instead of t''[s] – chris Jun 25 '18 at 10:08
• you also need to define your coordinates e.g. coord = {r, \[Theta], \[Phi], t} – chris Jun 25 '18 at 10:09
• Then Table[D[coord[[i]][s], s, s] + Sum[christ[i, h, k] D[coord[[h]][s], s] D[coord[[k]][s], s], {h, 1, 4}, {k, 1, 4}], {i, {1, 2, 3, 4}}] // FullSimplify // TableForm  this should work. – chris Jun 25 '18 at 10:12

I think the following works (check that it is correct please).

schwar = {{-(1 - (2 G M)/(r c^2))^-1, 0, 0, 0}, {0, -r^2, 0, 0}, {0,
0, -r^2 Sin[\[Theta]]^2, 0}, {0, 0, 0, (1 - (2 G M)/(c^2 r))}};
inverseschwar = Inverse[schwar];
christ[a_, b_, c_] :=
christ[a, b, c] =
Sum[1/2 inverseschwar[[a, d]] (D[schwar[[d, c]], coord[[b]]] +
D[schwar[[d, b]], coord[[c]]] -
D[schwar[[b, c]], coord[[d]]]), {d, 1, 4}]


Then

christ[1, 1, 1]

coord = {r, \[Theta], \[Phi], t}

Table[D[coord[[i]][s], s, s] +
Sum[christ[i, h, k] D[coord[[h]][s], s] D[coord[[k]][s], s], {h,
1, 4}, {k, 1, 4}] == 0, {i, {1, 2, 3, 4}}] //
FullSimplify // TableForm


Update

If you want the functions to appear explicitly in the solution

coord = {r, \[Theta], \[Phi], t};
rc = Thread[coord -> Map[#[s] &, coord]];
Clear[christ]
christ[a_, b_, c_] :=
Sum[1/2 inverseschwar[[a, d]] (D[schwar[[d, c]], coord[[b]]] +
D[schwar[[d, b]], coord[[c]]] -
D[schwar[[b, c]], coord[[d]]]), {d, 1, 4}] /. rc
christ[1, 1, 1]

Table[D[coord[[i]][s], s, s] +
Sum[christ[i, h, k] D[coord[[h]][s], s] D[coord[[k]][s], s], {h,
1, 4}, {k, 1, 4}] == 0, {i, {1, 2, 3, 4}}] // FullSimplify

• This is exactly what I faced. See the differential equations that is the final result has some terms without the $s$ argument. These are the terms from the Christoffel symbols. Let me rephrase my question. "Is there a way to print the Christoffel symbols so that each coordinate variable in it has a $s$ argument to it" Example: look at the coordinate $r$ in the denominator of all terms it stands alone without $s$ argument. – シャシュワト Jun 25 '18 at 10:25
• Does putting /. {r -> r[s], \[Theta] -> \[Theta][s], \[Phi] -> \[Phi][s], t -> t[s]} at the end of the definition for christ solve your issue? – KraZug Jun 25 '18 at 10:36
• YES! IT DID. Thanks!! Completely solved the problem I had. – シャシュワト Jun 25 '18 at 10:46
• I see that you did not read my answer on mathematica.stackexchange.com/questions/175933/… – Alex Trounev Jun 25 '18 at 17:32
• I'm sorry Chris, it was a message for Shashwat Sharan. – Alex Trounev Jun 26 '18 at 1:47