# 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 == 1.,
theta == 1.,
phi == 1.,
Derivative[r] == 0.,
Derivative[theta] == 0.5,
Derivative[phi] == 1.,
Derivative[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[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[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][], S[lambda][], S[lambda][]], {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[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, t' == 1,
r == 1., \[Theta] == 1., \[Phi] == 1.,
r' == 0., \[Theta]' == 0.5, \[Phi]' == 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[], S1[], S1[]], {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. – Quarkly Dec 29 '19 at 18:29
• You used Christoffel symbols for spherical coordinates in 3D, but not for a sphere. – Alex Trounev Dec 29 '19 at 18:38
• Your equations describe a straight line in 3D (geodesic in empty space). – Alex Trounev Dec 29 '19 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}}; – Quarkly Dec 29 '19 at 21:29
• You wrote down the metric of empty space in spherical coordinates. I described the metric of a sphere immersed in 4D. – Alex Trounev Dec 29 '19 at 21:36