Shortest path between two points on a surface

Anthony Carapetis does a good job of describing geodesics as shortest paths. I want to address you equation of how the surface constraint is used to actually calculate the shape of geodesics. For this, we will use calculus of variations. In single variable calculus, we learn how to optimize a functional with respect to a variable. In multivariable calculus, we learn how to do the same with respect to several variables. But what about an infinite number of variables? Suppose we had a function $J$ which has the form $$J(f)=\int_a^b F(x,f(x),f'(x))\,dx.$$ Note that $J$ takes a differentiable function $f$ as its input, and gives a real number as its output. One of the questions we can attempt to answer with calculus of variations is "what function $f$ gives the minimum value of $J$?" (a function of the type of $J$ is called a functional). This is the kind of problem that you are asking. In your case, the problem is what is the minimum value of $$J(x,y,z)=\int_0^T\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}\,dt.$$ The only difference is that $J$ is a functional of three functions $x$, $y$, and $z$. To solve a problem like this in regular calculus, we would take a derivative or gradient, and set it equal to zero. We do the same thing here, only the derivative we take is a bit different. It is called the variational derivative. If the functional takes in only one function as input, the derivative takes the form: $$\frac{\partial F}{\partial f}-\frac{d}{dx}\frac{\partial F}{\partial f'}.$$ In this derivative we are taking partial derivatives with respect to the functions $f$ and $f'$. This is done just like $f$ and $f'$ are normal variables independent of each other. Notice also that we don't take the variational derivative of $J$ directly, but rather of $F$, the function inside the integral. Since out functional accepts three functions, we have to solve the simultaneous equations:

$$\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial x}-\frac{d}{dt}\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial x'}=0,$$

$$\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial y}-\frac{d}{dt}\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial y'}=0,$$

$$\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial z}-\frac{d}{dt}\frac{\partial \sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}{\partial z'}=0.$$

But wait! we haven't used the constraint that the solution to these equations has to lie in the surface. Indeed, if you solve the above equations (and you can, it's not too hard), you will get $x''(t)=y''(t)=z''(t)=0$, a straight line. To force the solution to lie in a surface, we have to do something similar to Lagrange's method of optimizing constrained multivariable optimization problems. Instead of applying the variational derivative to $F$, we apply it to $F-\lambda f$, for some number $\lambda$ (remember $f(x,y,z)=0$ represents the surface). So we need to solve

$$\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-\lambda f(x,y,z))}{\partial x}-\frac{d}{dt}\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-\lambda f(x,y,z))}{\partial x'}=0,$$

$$\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-f(x,y,z))}{\partial y}-\frac{d}{dt}\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-f(x,y,z))}{\partial y'}=0,$$

$$\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-\lambda f(x,y,z))}{\partial z}-\frac{d}{dt}\frac{\partial (\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}-\lambda f(x,y,z))}{\partial z'}=0.$$

Let's simplify this. Remember that we are treating $x$ and $x'$ as independent, so $\frac{\partial x}{\partial x'}=0$ and $\frac{\partial x'}{\partial x}=0$. The same for $y$ and $z$, and all mixed derivatives. Thus, the problem reduces to solving

$$-\lambda \frac{\partial f(x,y,z)}{\partial x}-\frac{d}{dt}\frac{x'(t)}{\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}}=0,$$

and the same for $y$ and $z$. This is really hard to solve. It becomes easier if we assume something about the parameterization. The parameterization is how fast we sweep through the curve. There are multiple parameterizations that give the same curve. There is a parameterization called the arc-length parameterization which is a constant speed. If we assume that the solution is arc-length parameterized, the derivative arc length with respect to time is constant, so there is a constant $C$ so that $\sqrt{x'(t)^2+y'(t)^2+z'(t)^2}=C$. Thus, the above equation simplifies to

$$-\lambda \frac{\partial f(x,y,z)}{\partial x}-\frac{x''(t)}{C}=0.$$

Again, the same equations for $y$ and $z$. Thus, we see that we just need to solve $-C\lambda \nabla f=\gamma''(t)$, where $\gamma(t)$ is the curve. This differential equation is still hard to solve, but at least we have a shot. Good luck!

You're looking for a minimizing geodesic. In particular (assuming $f$ is nice enough for $S$ to be a regular surface) the shortest curve must satisfy the geodesic equation, which in this case reduces to the condition that $\gamma''(t)$ is parallel to $\nabla f$.

This turns the problem into a second-order ODE with two boundary conditions, but you will only find closed form solutions in very special cases. Even once you have found all the solutions joining the two points, you still need to choose the shortest one - in general there are curves that satisfy the geodesic equation but are not minimizing, and there may not even be a unique shortest one.