The Wikipedia article about geodesics talks about the equivalence of obtaining the geodesic by either minimizing the length functional $L$, or by minimizing the energy functional $L^2/2$, cf. the Phys.SE post Geodesic Equation from variation: Is the squared lagrangian equivalent?
In my courses so far, I haven't heard about the latter technique, so I was curious to see whether both approaches yield the same geodesic. I chose the geodesics on the unit sphere because it has the well-known solution of great circles.
Minimizing the length functional:
Given the spherical symmetry of the problem, start with the line element in spherical coordinates:
\begin{equation} ds^2=dr^2+r^2d\theta^2+r^2\sin^2\theta\,d\varphi^2 \end{equation}
For $r=1$ and factoring out $d\theta$, we get the arc length as the integral over $ds$: \begin{equation} l=\int ds=\int\limits_{\theta_1}^{\theta_2}\underbrace{\sqrt{1+\sin^2\theta\,\varphi'(\theta)^2}}_{\displaystyle L}\,d\theta \end{equation}
Then, plug in the Lagrangian $L$ in the Euler-Lagrange equation \begin{equation} \frac{\partial L}{\partial \varphi}-\frac{d}{d\theta}\frac{\partial L}{\partial \varphi'}=0 \end{equation}
until we arrive at $\displaystyle\varphi(\theta)=\int\frac{C_1}{\sin\theta\sqrt{\sin^2\theta-C_1^2}}\,d\theta$ . This result can be simplified using eq. (16) to yield
\begin{equation} \varphi(\theta)=-\arcsin\left(\frac{C_1\cot(\theta)}{\sqrt{1-C_1^2}}\right)+C_2 \end{equation}
which, as far as I understand, is the necessary condition for the curve to be a geodesic on the sphere. Insert into the parametrization of the unit sphere: \begin{equation} \theta\mapsto\begin{pmatrix} \sin\theta\cos\varphi(\theta) \\ \sin\theta\sin\varphi(\theta) \\ \cos\theta \end{pmatrix} \end{equation}
Mathematica produces the following plot:
Minimizing the energy functional:
Perform the same calculation as above using the squared Lagrangian, multiplied by $1/2$: \begin{equation} L^*=\frac{L^2}{2}=\frac{1}{2}\left(1+\sin^2\theta\,\varphi'(\theta)^2\right) \end{equation}
The E-L equation yields \begin{equation} \varphi(\theta)=\int \frac{C_1}{\sin^2\theta}\,d\theta=-C_1\cot(\theta)+C_2 \end{equation}
This time, the resulting curve does not look like a great circle whenever $C_1\neq 0$. What did I miss?
Mathematica code for the first plot:
\[CurlyPhi][\[Theta]_, C1_, C2_] := -ArcSin[(C1 Cot[\[Theta]])/Sqrt[1 - C1^2]] + C2
Manipulate[Show[ParametricPlot3D[sphere[{\[Theta], \[CurlyPhi]}], {\[Theta], 0, \[Pi]}, {\[CurlyPhi], 0, 2 \[Pi]}, PlotStyle -> {White, Opacity[0.6]}, ImageSize -> 500, AxesLabel -> {x, y, z}], ParametricPlot3D[{Sin[\[Theta]] Cos[\[CurlyPhi][\[Theta], C1, C2]], Sin[\[Theta]] Sin[\[CurlyPhi][\[Theta], C1, C2]], Cos[\[Theta]]}, {\[Theta], 0, t}, PlotStyle -> Red, PlotPoints -> 100, MaxRecursion -> 5]], {{C1, 0}, -1, 1, Appearance -> "Labeled"}, {{C2, 0}, -2 \[Pi], 2 \[Pi], Appearance -> "Labeled"}, {{t, \[Pi]/2}, 0, 2 \[Pi], Appearance -> "Labeled"}]