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.