# Please explain what's going on with this Geodesic Equation of a Sphere

I'm using the Christoffel Symbols found on this link to generate a set of three coupled differential equations as solutions to the Geodesic Equation of a Sphere, I have: $$\frac{d^2r}{d\lambda^2}-r\frac{d\theta}{d\lambda}\frac{d\theta}{d\lambda}-r\space Sin(\theta)^2\frac{d\phi}{d\lambda}\frac{d\phi}{d\lambda}=0$$ $$\frac{d^2\theta}{d\lambda^2}+\frac{2}{r}\frac{dr}{d\lambda}\frac{d\theta}{d\lambda}-Cos(\theta)Sin(\theta)\frac{d\phi}{d\lambda}\frac{d\phi}{d\lambda}=0$$ $$\frac{d^2\phi}{d\lambda^2}+\frac{2}{r}\frac{d\theta}{d\lambda}\frac{dr}{d\lambda}+2\space Cot(\theta)\frac{d\phi}{d\lambda}\frac{d\theta}{d\lambda}=0$$ So far, so good. My eventual goal here is to find the geodesic on a sphere that expands with time (that is, a preliminary step to a full FLRW model). Now I want to program this formula in Mathematica to see what happens on the surface of a sphere. Here's the code:

X[r_, theta_, phi_] := r*Cos[phi]*Sin[theta];
Y[r_, theta_, phi_] := r*Sin[phi]*Sin[theta];
Z[r_, theta_, phi_] := r*Cos[theta];
R[r_, theta_, phi_] := {X[r, theta, phi], Y[r, theta, phi], Z[r, theta, phi]};

equations =
{
r[0] == 1.,
theta[0] == 1.,
phi[0] == 1.,
Derivative[1][r][0] == 0.,
Derivative[1][theta][0] == 0.5,
Derivative[1][phi][0] == 1.,
Derivative[2][r][lambda] - r[lambda]*D[theta[lambda], lambda]*D[theta[lambda], lambda] - r[lambda]*Sin[theta[lambda]]^2*D[phi[lambda], lambda]*D[phi[lambda], lambda] == 0,
Derivative[2][theta][lambda] + (2/r[lambda])*D[r[lambda], lambda]*D[theta[lambda], lambda] - Cos[theta[lambda]]*Sin[theta[lambda]]*D[phi[lambda], lambda]*D[phi[lambda], lambda] == 0,
Derivative[2][phi][lambda] + (2/r[lambda])*D[phi[lambda], lambda]*D[r[lambda], lambda] + 2*Cot[theta[lambda]]*D[phi[lambda], lambda]*D[theta[lambda], lambda] == 0
};
S[lambda_] = {r[lambda], theta[lambda], phi[lambda]} /. First[NDSolve[equations, {r, theta, phi}, {lambda, 0, 5}]];
surface = ParametricPlot3D[{R[1., theta, phi]}, {theta, 0, Pi}, {phi, 0, 2*Pi}];
Show[{surface, ParametricPlot3D[R[S[lambda][[1]], S[lambda][[2]], S[lambda][[3]]], {lambda, 0, 5},
PlotStyle -> Directive[Red, Thick]]}, ImageSize -> Medium]


Now, here's where I need some help. If I plot this, I get:

Which appears to be a tangent vector. I was expecting the $$r$$ value to remain constant around the sphere since I set the initial distance to 1. and the velocity to 0. If I change the second derivative formula to:

  Derivative[2][r][lambda] == 0,


Then I get the expected Geodesic: Can someone explain to me what is going on? Did I make a mistake? Or should I not expect the $$r$$ coordinate to follow the Geodesic because it's normal to the surface?

There is an obvious typo in the derivation of the equations of motion on a sphere. First, we write down the coordinates and the metric tensor for a sphere with a given radius r0=1 immersed in 4D:

n = 4; r0=1;
coord = {t, r, \[Theta], \[Phi]};
metric ={{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, r0^2, 0}, {0, 0, 0,
r0^2 Sin[\[Theta]]^2}};


We calculate the Christoffel symbols for this case

inversemetric = Simplify[Inverse[metric]]

affine :=
affine = Simplify[
Table[(1/2)*
Sum[inversemetric[[i,
s]]*(D[metric[[s, j]], coord[[k]]] +
D[metric[[s, k]], coord[[j]]] -
D[metric[[j, k]], coord[[s]]]), {s, 1, n}], {i, 1, n}, {j,
1, n}, {k, 1, n}]];

listaffine :=
Table[If[UnsameQ[affine[[i, j, k]],
0], {ToString[\[CapitalGamma][i, j, k]], affine[[i, j, k]]}], {i,
1, n}, {j, 1, n}, {k, 1, j}];

TableForm[Partition[DeleteCases[Flatten[listaffine], Null], 2],
TableSpacing -> {2, 2}]


And so there are two non-zero components $$\Gamma[3,4,4]=-\cos\theta \sin\theta, \Gamma[4,4,3]=\cot \theta$$ Now we find the equations of motion $$\frac {d^2x^i}{ds^2}+\Gamma ^i_{kl}\frac{dx^k}{ds}\frac {dx^l}{ds}=0$$ Put

v = {t'[s], r'[s], \[Theta]'[s], \[Phi]'[s]};


And then calculate

Table[D[v[[i]], s] + ((affine.v).v)[[i]] == 0, {i, n}]


Finally, we have a system of equations

    eqs = {t''[s] == 0,
r''[s] ==
0, -Cos[\[Theta][s]] Sin[\[Theta][s]] \[Phi]'[s]^2 + \[Theta]''[
s] == 0,
2 Cot[\[Theta][s]] \[Theta]'[s] \[Phi]'[s] + \[Phi]''[s] == 0};

ic = {t[0] == 0, t'[0] == 1,
r[0] == 1., \[Theta][0] == 1., \[Phi][0] == 1.,
r'[0] == 0., \[Theta]'[0] == 0.5, \[Phi]'[0] == 1.};


Now we can display the trajectory on the sphere

X[r_, theta_, phi_] := r*Cos[phi]*Sin[theta];
Y[r_, theta_, phi_] := r*Sin[phi]*Sin[theta];
Z[r_, theta_, phi_] := r*Cos[theta];
R[r_, theta_, phi_] := {X[r, theta, phi], Y[r, theta, phi],
Z[r, theta, phi]};

S1 = NDSolveValue[
Flatten[{eqs, ic}], {r[s], \[Theta][s], \[Phi][s]}, {s, 0, 5}];
surface =
ParametricPlot3D[{R[1., theta, phi]}, {theta, 0, Pi}, {phi, 0,
2*Pi}];
Show[{surface,
ParametricPlot3D[R[S1[[1]], S1[[2]], S1[[3]]], {s, 0, 5},
PlotStyle -> Directive[Red, Thick]]}, ImageSize -> Medium]


• Could you direct me to my typo, please? Your solution is excellent, but it would be useful to me to see where my original solution broke. Dec 29, 2019 at 18:29
• You used Christoffel symbols for spherical coordinates in 3D, but not for a sphere. Dec 29, 2019 at 18:38
• Your equations describe a straight line in 3D (geodesic in empty space). Dec 29, 2019 at 21:19
• I believe your metric is incorrect. It should be: metric = {{1, 0, 0}, {0, r^2, 0}, {0, 0, r^2 Sin[\[Theta]]^2}}; Dec 29, 2019 at 21:29
• Alright, this is the concept I think I was missing. A geodesic on a 2D surface embedded on a 3D sphere would trace a great circle. But the geodesic in 3D space - without something like gravity to bend it - would be a straight line in both cartesian and polar coordinates. All I was doing was switching the coordinates to polar in my original post. Dec 29, 2019 at 22:29