Solution (source code)

= Solution
{c}

Independent endpoint-fixed variations of the two coordinates give the <Euler-Lagrange equations>
$$
\frac d{dt}F_{\dot u}-F_u=0,\qquad
\frac d{dt}F_{\dot v}-F_v=0.
$$
Differentiating the sphere parametrization gives the <Riemannian metric> $ds^2=a^2(d\theta^2+\sin^2\theta\,d\phi^2)$. Its energy functional is the stated $I$. Its stationary paths are <geodesics> with affine, constant-speed parametrization; conversely a geodesic can be so parametrized. This explains the use of energy rather than arc length.

For $F=a^2(\dot\theta^2+\sin^2\theta\dot\phi^2)/2$, the equations are
$$
\boxed{\ddot\theta-\sin\theta\cos\theta\,\dot\phi^2=0,\qquad
\frac d{dt}(\sin^2\theta\,\dot\phi)=0.}
$$
The second gives the conserved <generalized momentum> $\sin^2\theta\dot\phi=h$. Since $F$ has no explicit time dependence, its <Beltrami identity> gives constant $\dot\theta F_{\dot\theta}+\dot\phi F_{\dot\phi}-F=F$. Dividing by $a^2/2$ gives the second first integral,
$$
\boxed{\dot\theta^2+\sin^2\theta\dot\phi^2=E.}
$$
Thus $E$ is the squared speed divided by $a^2$, while $h$ records the rotational symmetry about the polar axis. The coordinate equations apply away from the polar coordinate singularities.