Geodesics on a Saddle Surface

Question

Does anyone know of any free and relatively easy-to-use math apps like Geogebra that can be programmed (or were explicitly designed) to draw the geodesic between any two selected points on a saddle surface like the one below?

Thanks in advance.

Answer 1

If $g(x) = (x,y(x),z(x))$ is such a geodesic over the surface $H(x,y(x),z(x))=0$ then considering the lagrangian

$$ L(x,y(x),z(x),y'(x),z'(x))= \sqrt{1+y'^2+z'^2}+\lambda H(x,y(x),z(x)) $$

we have

$$ \frac{d}{dx}\left(\frac{\partial L}{\partial \theta'}\right)-\frac{\partial L}{\partial \theta}=0 $$

where $\theta = (y,z)$, giving

$$ \left\{ \begin{array}{l} \frac{y''(x)}{\sqrt{y'(x)^2+z'(x)^2+1}}-\frac{y'(x) \left(2 y'(x) y''(x)+2 z'(x) z''(x)\right)}{2 \left(y'(x)^2+z'(x)^2+1\right)^{3/2}}-\lambda H_y = 0\\ \frac{z''(x)}{\sqrt{y'(x)^2+z'(x)^2+1}}-\frac{z'(x) \left(2 y'(x) y''(x)+2 z'(x) z''(x)\right)}{2 \left(y'(x)^2+z'(x)^2+1\right)^{3/2}} -\lambda H_z=0\\ H(x,y(x),z(x))=0 \end{array} \right. $$

now assuming regularity on $H$ we have equivalently

$$ \frac{d^2}{d x^2}H = y''(x) H_y+2 z'(x) \left(y'(x) H_{y,z}+H_{x,z}\right)+y'(x)^2 H_{y,y}+2 y'(x) H_{x,y}+z''(x) H_z+z'(x)^2 H_{z,z}+H_{x,x}=0 $$

so the ode system reads

$$ \left\{ \begin{array}{l} \frac{y''(x)}{\sqrt{y'(x)^2+z'(x)^2+1}}-\frac{y'(x) \left(2 y'(x) y''(x)+2 z'(x) z''(x)\right)}{2 \left(y'(x)^2+z'(x)^2+1\right)^{3/2}}-\lambda H_y = 0\\ \frac{z''(x)}{\sqrt{y'(x)^2+z'(x)^2+1}}-\frac{z'(x) \left(2 y'(x) y''(x)+2 z'(x) z''(x)\right)}{2 \left(y'(x)^2+z'(x)^2+1\right)^{3/2}} -\lambda H_z=0\\ y''(x) H_y+2 z'(x) \left(y'(x) H_{y,z}+H_{x,z}\right)+y'(x)^2 H_{y,y}+2 y'(x) H_{x,y}+z''(x) H_z+z'(x)^2 H_{z,z}+H_{x,x}=0 \end{array} \right. $$

Solving this system for $(y'', z'',\lambda)$ we obtain finally the geodesic equations

$$ \cases{ y''=g_1(x,y,y',z,z')\\ z''=g_2(x,y,y',z,z') } $$

which computed between $(x_1,y(x_1),z(x_1))$ and $(x_2,y(x_2),z(x_2))$ give the geodesic orbit. Note that here $H(x_i,y(x_i),z(x_i))=0$. As an application, for $H(x,y,z) = x y-z = 0$ we have

$$ \left\{ \begin{array}{l} y''(x)+\frac{2 y'(x) \left(-y'(x) z'(x)+x y'(x)^2+x\right)}{x^2 y'(x)^2+x^2-2 x y'(x) z'(x)+z'(x)^2+1} = 0\\ z''(x)+\frac{2 y'(x) \left(x y'(x) z'(x)-z'(x)^2-1\right)}{x^2 y'(x)^2+x^2-2 x y'(x) z'(x)+z'(x)^2+1}=0 \\ \end{array} \right. $$

Attached a MATHEMATICA script which solves the ode

xmin = -2;
xmax = 2;
ymin = -2;
ymax = 2;
zmin = -2;
zmax = 2;
frac = 2;
H[x_, y_, z_] := x y - z
gr0 = ContourPlot3D[H[x, y, z] == 0, {x, xmin, xmax}, {y, ymin, ymax}, {z, zmin, zmax},ContourStyle -> Opacity[0.3]];
L = Sqrt[y'[x]^2 + z'[x]^2 + 1] + lambda H[x, y[x], z[x]];
equs = D[Grad[L, {y'[x], z'[x]}], x] - Grad[L, {y[x], z[x]}];
equstot = Join[equs, {D[H[x, y[x], z[x]], x, x]}];
odes = Thread[{y''[x], z''[x]} == ({y''[x], z''[x]} /. Solve[equstot == 0, {y''[x], z''[x], lambda}][[1]])];

path[p1_, p2_] := Module[{x1 = p1[[1]], y1 = p1[[2]], x2 = p2[[1]], y2 = p2[[2]], sol},
 sol = NDSolve[{odes, y[x1] == y1, z[x1] == x1 y1, y[x2] == y2, z[x2] == x2 y2}, {y, z}, {x, x1, x2}];
  Return[Evaluate[{x, y[x], z[x]} /. sol[[1]]]]
  ]

y1 = RandomReal[{ymin, ymax}/frac];
y2 = RandomReal[{ymin, ymax}/frac];
p1 = {xmin/frac, y1, xmin y1/frac};
p2 = {xmax/frac, y2, xmax y2/frac};
path0 = path[{xmin/frac, y1}, {xmax/frac, y2}];
grpts = {Graphics3D[{Red, PointSize[0.02], Point[p1]}], 
         Graphics3D[{Red, PointSize[0.02], Point[p2]}]};
gr1 = ParametricPlot3D[path0, {x, xmin/frac, xmax/frac}, PlotStyle -> {Thick, Blue}];
Show[gr0, gr1, grpts, Axes -> False]

Answer 2

The hyperbolic paraboloid is the surface given by $$ z=x^2-y^2\,. $$ Then $$ dz=2x\,dx-2y\,dy\,. $$ This leads to the metric $$ ds^2=dx^2+dy^2+dz^2=(1+4x^2)\,dx^2+(1+4y^2)\,dy^2-8xy\,dx\,dy\,, $$ or $$ g=\begin{pmatrix}1+4x^2&-4xy\\-4xy&1+4y^2\end{pmatrix}\,. $$ The python package sympy.diffgeom can be used to calculate the Christoffel symbols: \begin{align*} {\Gamma^x}_{xx} &= \frac{4 x}{4 x^{2} + 4 y^{2} + 1}\,,& {\Gamma^x}_{yy} &= - \frac{4 x}{4 x^{2} + 4 y^{2} + 1}\,,\\ {\Gamma^y}_{xx} &= - \frac{4 y}{4 x^{2} + 4 y^{2} + 1}\,,& {\Gamma^y}_{yy} &= \frac{4 y}{4 x^{2} + 4 y^{2} + 1}\,. \end{align*} These lead to the geodesic equations

\begin{align*} 0&=\ddot x+\frac{4x}{4x^2+4y^2+1}(\dot x^2-\dot y^2)\,,\\[2mm] 0&=\ddot y-\frac{4y}{4x^2+4y^2+1}(\dot x^2-\dot y^2)\,. \end{align*}

• Note that Cesareo's equations in the other answer are for $(y(x),z(x))$ with parameter $x$ while the ones obtained here are for the curve $t\mapsto(x(t),y(t))$ with a "time parameter" $t$ that we don't see in the graphs. Once $(x(t),y(t))$ is known we simply set $z(t)=x^2(t)-y^2(t)$ to get the curve on the saddle.

• The simpler looking geodesic equations here can be solved numerically given an initial position $(x(0),y(0))$ and an initial velocity $(\dot x(0),\dot y(0))\,.$

• To find the solution for given initial and final condition seems like a whole subject in itself.

• The graph below shows a geodesic ($\color{green}{green}$) with initial position $(3,-5)$ and initial velocity. $(-0.9,1)\,.$ Its length is approximately $24.309\,.$ The two neighbouring curves ($\color{blue}{blue}$ and $\color{red}{red}$) have approximate lengths $24.316\,.$

To keep things simple I give you only the python code that produces the green geodesic.

import matplotlib.pyplot as plt
from numpy import array, linspace, meshgrid, sqrt

Limit = 5
n = 1000
dt = 1/n

def solve( x0, y0, xp0, yp0 ):
    x = array([x0,x0 + xp0*dt])
    y = array([y0,y0 + yp0*dt])
    i = 2
    T = 0
    while 1:
        xp = (x[i-1]-x[i-2])/dt
        yp = (y[i-1]-y[i-2])/dt
        xpp = -4*x[i-1]*(xp**2-yp**2)/(4*x[i-1]**2+4*y[i-1]**2+1)
        ypp = +4*y[i-1]*(xp**2-yp**2)/(4*x[i-1]**2+4*y[i-1]**2+1)
        if i >= len(x):
           x.resize(i+100)
           y.resize(i+100)
        x[i] = (xpp*dt**2 + 2*x[i-1] - x[i-2])
        y[i] = (ypp*dt**2 + 2*y[i-1] - y[i-2])
        if x[i] < Limit and x[i] > -Limit and y[i] < Limit and y[i] > -Limit:
           i += 1
           T += dt
        else:
           break
        x.resize(i)
        y.resize(i)
      return x, y, T

fig = plt.figure(figsize=[10,12])
ax = fig.gca(projection='3d')
u = linspace(-Limit,Limit,100)
v = linspace(-Limit,Limit,100)

X, Y = meshgrid(u, v)
Z = X**2-Y**2

ax.plot_wireframe(X, Y, Z, linewidth=0.5, color='y')
x0 = 3
y0 = -5
xp0 = -0.9
yp0 = 1
x, y, T = solve( x0, y0, xp0, yp0 )    
ax.plot(x,y,x**2-y**2,linewidth=1,color='g')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.set_zlabel(r'$z$')

plt.show()