Solution (source code)

= Solution

For a parametrized <curve> on the surface, the <energy of a curve> is
$$
\mathcal E(\gamma)=\frac12\int_a^b|\dot\gamma(t)|^2\,dt=\frac12\int_a^b\left(E\dot u^2+2F\dot u\dot v+G\dot v^2\right)\,dt,
$$
where the <first fundamental form> coefficients are evaluated at $(u(t),v(t))$. A <fixed-endpoint variation of a curve> has coordinate variation fields $\eta,\xi$ vanishing at $a,b$. Differentiating under the integral and integrating their derivative terms by parts gives
$$
\delta\mathcal E=\int_a^b\!\left[\left(\frac{E_u\dot u^2+2F_u\dot u\dot v+G_u\dot v^2}{2}-\frac d{dt}(E\dot u+F\dot v)\right)\eta+\left(\frac{E_v\dot u^2+2F_v\dot u\dot v+G_v\dot v^2}{2}-\frac d{dt}(F\dot u+G\dot v)\right)\xi\right]dt.
$$
The fields can be chosen independently with compact support in $(a,b)$. Stationarity and the <fundamental lemma of the calculus of variations> therefore force both coefficients to vanish. These <Euler-Lagrange equations> are the <geodesic equations>:
$$
\boxed{\frac d{dt}(E\dot u+F\dot v)=\frac12(E_u\dot u^2+2F_u\dot u\dot v+G_u\dot v^2),\qquad\frac d{dt}(F\dot u+G\dot v)=\frac12(E_v\dot u^2+2F_v\dot u\dot v+G_v\dot v^2).}
$$