Geodesics on a hyperbolic paraboloid

It is standard differential geometry to find the differential equation for the geodesics on this surface. (But I could easily have made a mistake in the calculation anyway.) Since it is a complete negatively curved surface, there is exactly one geodesic connecting any two points. You have a curve $\vec{p}(t) = (x(t),y(t),z(t))$ on the surface $z = xy$. The geodesic equation is $$\vec{p}''(t) \propto \vec{\nabla}z \oplus -1$$ (the acceleration is perpendicular to the surface), which expands to $$(x'',y'',x''y+2x'y' + xy'') \propto (y,x,-1).$$ Thus $$\frac{x''}{y} = \frac{y''}{x} = -(x''y+2x'y'+xy'').$$ That is an equation that Mathematica can solve. The only tricky part is to solve it with boundary conditions at both ends, where you may as well assume that the endpoints are at $t=0$ and $t=1$. For that purpose it could be better to minimize the curve's energy, by definition $$E[\vec{p}] = \int_0^1 |\vec{p}'(t)|^2 dt = \int_0^1 (x'^2+y'^2 + (x'y+xy')^2) dt.$$ I don't know the most convenient way to do this numerically, but somehow it should be possible. Again, since the surface is negatively curved, this energy functional is well-behaved.

I don't know whether there is a closed form solution in elementary functions. There is a closed form solution for a round paraboloid, but it's messy.

I was working on answering another question involving integrating the geodesic equations on a surface, and the links there lead me back to this question, which I hadn't noticed before. In case anyone is interested, there is a more direct way to get to the geodesics of the induced metric on the hyperbolic paraboloid, so I thought that would go ahead and input this approach.

As Greg pointed out, the geodesic equations in the graphing coordinates are $$ x'' = -\frac{2yx'y'}{(1+x^2+y^2)}\qquad\text{and}\qquad y'' = -\frac{2xx'y'}{(1+x^2+y^2)}. $$ These equations have two first integrals that are quadratic in velocity: $$ Q_1 = (1+y^2)\,x'^2 + 2xy\,x'\,y' + (1+x^2)\,y'^2, $$ i.e., the induced metric itself, which is, of course, the statement that geodesics have constant speed, and $$ Q_2 = 2(1+x^2+y^2)\,x'\,y', $$ which is easily verified. (It's a classical fact that the geodesic equations on any (non-flat) surface of degree $2$ in Euclidean $3$-space have a second first integral that is quadratic in velocities, namely $Q_2 = |K|^{-3/4}I\!I$ in addition to the obvious first integral $Q_1 = I$.)

Since there are two independent first integrals quadratic in velocity, this is a Liouville metric and hence can be put in Liouville form. This is an algorithmic procedure; following it leads to the result that, if one establishes new coordinates $z$ and $w$ on the surface by $$ x = \sinh \left(\frac{z+w}2\right)\quad\text{and}\quad y = \sinh\left(\frac{z-w}2\right), $$ then we have the first integrals $$ \begin{aligned} Q_1 &= \tfrac14\bigl(\cosh(z)+\cosh(w)\bigr)\bigl(\,\cosh(z)\,z'^2 +\cosh(w)\,w'^2\,\bigr)\\ Q_2 &= \tfrac14\bigl(\cosh(z)+\cosh(w)\bigr)\bigl(\,\cosh(w)\cosh(z)\,z'^2 -\cosh(z)\cosh(w)\,w'^2\,\bigr) \end{aligned} $$ (The actual Liouville coordinates would be $(u,v)$ where $\mathrm{d}u = \sqrt{\cosh z}\,\mathrm{d}z$ and $\mathrm{d}v = \sqrt{\cosh w}\,\mathrm{d}w$, but these are elliptic integrals, and it seems pointless to change to these coordinates.) In particular, for a unit speed geodesic, i.e., one for which $Q_1 = 1$ and $Q_2 = c$ for some constant $c$, we have $Q_2-c\,Q_1 =0$, so $$ (\cosh(w)-c)\cosh(z)\,z'^2 -(\cosh(z)+c)\cosh(w)\,w'^2, $$ and one can thus separate variables, yielding $$ \frac{\cosh(z)\,\mathrm{d}z^2}{(\cosh(z)+c)} = \frac{\cosh(w)\,\mathrm{d}w^2}{(\cosh(w)-c)}. $$ Now, $c=0$ corresponds to $z\pm w$ being constant, i.e., $x$ or $y$ is constant, which are straight lines in the surface. When $c\not=0$, the geodesic has to stay in the region where $\cosh z + c$ and $\cosh w - c$ are non-negative, and we have an equation that can be integrated 'by quadratures' to yield two foliations by geodesics of the region where $\cosh z + c$ and $\cosh w - c$ are positive.
$$ \sqrt{\frac{\cosh z}{\cosh z + c}}\,\mathrm{d}z = \pm\, \sqrt{\frac{\cosh w}{\cosh w - c}}\,\mathrm{d}w $$ This will include a family that envelopes either $\cosh z + c = 0$ (if $c<-1$) or $\cosh w - c = 0$ (if $c>1$). (The curves $z = 0$ and $w=0$ are geodesics.)

In order explicitly to compute the distance between two points using these formulae, one would need to compute the corredponding $z$ and $w$ coordinates of the two points (easy), find the $c$ belonging to the (unique) geodesic joining those two points (nontrivial), and then compute the elliptic integrals as above.

I expect that the likelihood that one could actually carry this out and find a 'closed form' expression for the distance function on the surface is rather low.

I haven't worked out the details, but this should be doable in paraboloidal coordinates, in the same way that the geodesics on an ellipsoid were computed by Jacobi using his famous ellipsoidal coordinates (invented for that very purpose).

In the notation of the Wikipedia article, take the parameters to be $A=1$ and $B=-1$; then the coordinate surface $\mu=0$ is your paraboloid $2z=y^2-x^2$, and the geodesics on any coordinate surface should be possible to integrate (more or less) explicitly, although the details might be a bit messy. The idea is to rewrite the Hamiltonian for geodesic motion, $H=(p_x^2+p_y^2+p_z^2)/2$, in the new coordinates and apply the Hamiltion–Jacobi method (the Hamilton–Jacobi equation is separable in these coordinates).

I haven't got time to look for a good reference explaining Jacobi's work right now. I'll update the answer later if I find something.

Edit: Here's what I think is the most convenient, and least error-prone, way of setting up the equations for the geodesics. Change to paraboloidal coordinates in $R^3$; I'll call them $(u_1, u_2, u_3)$ instead of $(\lambda, \mu, \nu)$. Since this is an orthogonal coordinate system, the Euclidean metric tensor is diagonal, $ds^2=\sum_{k=1}^3 h_k^2 du_k^2$, where $h_1$, $h_2$, $h_3$ are the scale factors given in the Wikipedia article. The hyperbolic paraboloid that you are interested in is the coordinate surface $u_2=0$, which is a Riemannian manifold in itself, with coordinates $(u_1,u_3)$ and metric tensor given by $ds^2 = h_1^2 du_1^2 + h_3^2 du_3^2$ (where of course $u_2=0$ should be substituted into the expressions for the scale factors). On any Riemannian manifold, the geodesic equations are the canonical Hamiltonian equations given by the Hamiltonian function $H=\frac{1}{2} g^{ij} p_i p_j$, where $g^{ij}$ is the inverse metric tensor. In this case, we get $H(u_1,u_3,p_1,p_3) = \frac{1}{2} \left( \frac{p_1^2}{h_1(u_1,u_3)^2} + \frac{p_3^2}{h_3(u_1,u_3)^2} \right)$. So just feed this function to Mathematica, and numerically integrate the equations $$\dot{u}_1 = \partial H/\partial p_1,$$ $$\dot{u}_3 = \partial H/\partial p_3,$$ $$\dot{p}_1 = -\partial H/\partial u_1,$$ $$\dot{p}_3 = -\partial H/\partial u_3,$$ with suitable initial conditions. This will give you a geodesic emanating from a given point in a given direction. The result is in terms of paraboloidal coordinates, of course, but it is trivial to express it in terms of Cartesian coordinates (for plotting) using the formulas defining the change of variables. Finding a geodesic between two given points seems more complicated; perhaps use some "shooting" algorithm?

As I wrote above, it should be possible to integrate the equations by hand, but numerical integration seems to suffice for your purposes.

For the explicit integration of the geodesics on the ellipsoid, Jacobi's own lectures are probably as good a source as anything else (if you can read German). They are available at the Internet Archive. Elliptic coordinates are described in Lecture 26, the geodesics on the ellipsoid in Lecture 28 (p. 212).