# Practical drawing of geodesics

I want to use a computer to draw geodesics on a known parameterized surface of revolution, starting from a known point and at a known angle to the meridian.

What would be the easiest way of doing this?

Some methods I've considered:

• Using Clairaut's relation ($r \cos \theta = C$), and using small increments, but I'm afraid that this is not a general solution, as this will cause the geodesic to "get stuck" on parallels whenever the angle crosses zero.
• Use the relation that for a surface given by: $$\left(\varphi(v)\cos u, \varphi(v)\sin u, \psi(v)\right)$$ The geodesics are given by: $$\frac{du}{dv}=\frac{\sqrt{\varphi'^2+\psi'^2}}{\varphi\sqrt{\varphi^2-c^2}}\ \longrightarrow\ u(v_1) = u_0 + \int_{v_0}^{v_1} \frac{\sqrt{\varphi'^2+\psi'^2}}{\varphi\sqrt{\varphi^2-c^2}} dv$$ This presents the problem of finding the correct value of $c$ given the starting angle, and as user levap pointed out, has the same "getting stuck" problem, since $u(v)$ is no longer a function.
• Various numerical methods which work for any triangulated surface, but I feel it would be a shame to use these when I know the exact form of the surface.

Are there any easier methods or modifications to the above?

• On a surface of rotation, a geodesic that becomes parallel to a "latitude" at a point $p$ doesn't "cross" the latitude, but instead "turns back", tracing its own reflection in the meridian ("longitude") through $p$. Does that help? – Andrew D. Hwang Dec 12 '15 at 16:03
• @AndrewD.Hwang - That is very helpful - I didn't think of the symmetry involved. However, from a numerical point of view, it may be difficult to "hit" the inflection point exactly, so knowing it's position analytically would be helpful. – nbubis Dec 13 '15 at 8:14

Assume that we are given a curve $\gamma \colon J \rightarrow \mathbb{R}^3$ of the form $\gamma(v) = (\varphi(v), 0, \psi(v))$ such that $\frac{d \gamma}{dv}(v) \neq 0$ and $\varphi(v) > 0$ for all $v \in J$. The surface of revolution obtained by revolving $\gamma$ around the $z$-axis is

$$S = \{ (\varphi(v) \cos u, \varphi(v) \sin u, \psi(v)) \, | \, v \in J, u \in [0,2\pi] \}.$$

Since $S$ cannot be covered with a single coordinate chart, it is more useful to consider $S$ as the image of the map $X \colon \mathbb{R} \times J \rightarrow \mathbb{R}^3$ given by

$$X(u,v) = (\varphi(v) \cos(u), \varphi(v) \sin u, \psi (v)).$$

The map $X$ is a local diffeomorphism and the pullback of the first fundamental form from $S$ to $\mathbb{R} \times J$ is given by

$$g(u,v) = \begin{pmatrix} \varphi(v)^2 & 0 \\ 0 & \varphi'(v)^2 + \psi'(v)^2 \end{pmatrix} = \begin{pmatrix} E(v) & 0 \\ 0 & F(v) \end{pmatrix}.$$

Geodesics on $\mathbb{R} \times J$ with the pullback metric will map under $X$ to geodesics of $S$. Using $g$, we can calculate the Christoffel symbols and write explicitly the geodesic equation $\nabla_{\dot{\alpha}(t)}{\dot{\alpha}} = 0$ for a curve $\alpha(t) = (u(t), v(t))$ in $\mathbb{R} \times J$. Written shortly, the equations are

$$\ddot{u} + \frac{E'}{E} \dot{u} \dot{v} = 0, \\ \ddot{v} - \frac{E'}{2G} \dot{u}^2 + \frac{G'}{2G} \dot{v}^2 = 0.$$

More explicitly, we have

$$\ddot{u}(t) + \frac{2 \varphi(v(t)) \varphi'(v(t))}{\varphi(v(t))^2} \dot{u}(t) \dot{v}(t) = 0, \\ \ddot{v}(t) - \frac{\varphi(v(t)) \varphi'(v(t))}{\varphi'(v(t))^2 + \psi'(v(t))^2} \dot{u}(t)^2 + \frac{\varphi'(v(t)) \varphi''(v(t)) + \psi'(v(t)) \psi''(v(t))}{\varphi'(v(t))^2 + \psi'(v(t))^2} \dot{v}(t)^2 = 0. \label{equation1}\tag{1}$$

This is a system of (coupled) non-linear second order differential equations for $(u(t), v(t))$ which can be solved numerically given initial conditions at some (arbitrary) time $t = t_0$ of the form

$$\alpha(t_0) = (u_0, v_0), \,\,\, \dot{\alpha}(t_0) = (u_0', v_0')$$

and you can use them to find and draw the geodesic $X(\alpha(t))$. The geodesic will be of constant speed (with respect to the metric $g$) and won't be generally a graph (so $\alpha(t)$ won't be of the form $\alpha(t) = (t, v(t))$ or $\alpha(t) = (u(t), t)$).

The description above doesn't really use the fact that we are looking for geodesics on a surface of revolution that has a special symmetry. The special symmetry gives a conserved quality that must be constant along geodesics. This is called Clairaut's relation:

$$\varphi(v(t)) \cos \theta(t) = C$$

where $\theta(t)$ is the angle the geodesic $\alpha(t)$ makes with the parallel $v = v(t)$ and $C$ is some constant. In fact, the converse almost holds in the sense that if $\alpha$ is a regular curve that doesn't coincide with a parallel on any interval and satisfies Clairaut's relation then $\alpha$ is a reparametrization of a geodesic.

If we assume that $\alpha$ has the form of a graph $\alpha(u) = (u, v(u))$ (note that this is different than the equation you present in which $\alpha(v) = (v(u), u)$) and plug $\alpha$ into Clairaut's relation, we obtain the following equation for $v$:

$$\frac{ dv}{du} = \pm \frac{1}{C} \frac{\varphi(v(u)) \sqrt{\varphi^2(v(u)) - C^2}}{\sqrt{ \varphi'(v(u))^2 + \psi'(v(u))^2}} = F(v). \label{eq2}\tag{2}$$

Thus, instead of having two second order equations, we have reduced the problem to a single first order equation for $v$. Again, given initial conditions $v(u_0) = v_0$ with $\varphi(v_0) \neq \pm C$, the equation has a unique solution which can found numerically. However, this equation has the drawback that it cannot really used to find geodesics which, at $(u_0, v_0)$, are tangent to a parallel. The reason is that the problem is not well-posted - if $\alpha$ is tangent to a parallel at $(u_0, v_0)$ then $\cos (\theta) = 1$ and hence, $\varphi(v_0) = C$ and $F(v_0) = 0$. The right hand side is not a Lipschitz function and you don't have uniqueness of solutions. Indeed, the solution $v(u) \equiv v_0$ is a solution of ($\ref{eq2}$) but it is not necessarily a geodesic - it will be a geodesic if and only if $\varphi'(v_0) = 0$. If it is not a geodesic, then equation ($\ref{eq2}$) will have two solutions.

The advantage of seeking a parametrization of the form $(u, v(u))$ lies in the fact that any geodesic which is not a meridian can be parametrized in this way. A parametrization of the form $(u(v), v)$ is possible only for geodesics that doesn't become tangent to a parallel at some point. If we write the equation for $u(v)$ (which is the inverse function of $v(u)$), we get (more or less, up to a constant that seems to have been lost) the equation you bring in your second point. In some sense, this improves the situation because the equation you get for $u(v)$ is just a direct integral which can be evaluated numerically. However, it suffers from a similar problem as equation ($\ref{eq2}$) because if the geodesic becomes parallel to a tangent at some point, the integrand blows up and the derivative $\frac{du}{dv}$ approaches $\pm \infty$. This is not surprising, because one can't expect such a parametrization to hold near a point of tangency to a parallel.

Taking all this into account, I can offer a few different methods for drawing geodesics:

1. Just use equations ($\ref{equation1}$). They look more complicated but unless you will have performance issues, solving them numerically will get you want without worrying about whether the geodesic becomes tangent to a parallel or not.
2. Use equation ($\ref{eq2}$) and the symmetry involved. More explicitly, let us assume you want to find a geodesic that is not a meridian, emanating forward from $(u_0, v_0)$, whose angle with the parallel $v = v_0$ is $\cos \theta_0 \neq 1$. Set $C = \varphi(v_0) \cos \theta_0$ and solve equation ($\ref{eq2}$) numerically to get (some approximation of) $v(u)$. There are three possible cases:
• The geodesic won't become tangent to a parallel at any point.
• The geodesic won't become tangent to a parallel but will become asymptotic to a parallel. This means that the solution will satisfy $\lim_{u \to \infty} v(u) = v_1$ and $F(v(u)) \neq 0$ for all $u > u_0$. This is possible if and only if the parallel itself is a geodesic, which can be detected by checking whether $\varphi'(v_1) = 0$.
• The geodesic will become tangent to a parallel after some finite number of windings around the surface. This means that the solution will satisfy $\lim_{u \to u_1} v(u) = v_1$ and $F(v_1) = 0$ (so $v'(u_1) = 0$). To continue the solution beyond $u_1$, use the fact that the geodesic bounces back from the parallel tracing its reflection as Andrew D. Hwang notes. More formally, the solution will satisfy $v(u) = v(2u_1 - u)$ for $u_1 < u < u_1 + (u_1 - u_0)$. Then repeat the process.
3. The method described in the previous item does not allow you to find a geodesic that starts as tangent to a parallel. In order to do that, you can combine equations ($\ref{equation1}$) and ($\ref{eq2}$). Namely, assume that you want to find a geodesic $\alpha(t) = (u(t), v(t))$ satisfying $\alpha(t_0) = (u_0, v_0)$ and $\dot{\alpha}(t_0) = (1,0)$. Linearizing equations ($\ref{equation1}$), we have $$u(t) \approx u_0 + (t - t_0), \\ v(t) \approx v_0 + \frac{1}{2}\frac{\varphi(v_0) \varphi'(v_0)}{\varphi'(v_0)^2 + \psi'(v_0)^2} (t - t_0)^2$$

Use this to find $\alpha(t_1)$ for some small $t_1 > t_0$ and then use the method of the previous item to continue solving for the geodesic.

• Thanks! This kind of comes back to the same problem as using Clairaut's method - we can't cross parallels. What methods are there that solve this problem? – nbubis Dec 9 '15 at 15:59
• Any additional thoughts? – nbubis Dec 13 '15 at 12:29
• I've added some details about possible methods I thought about. – levap Dec 14 '15 at 8:36
• Well done! I'll try and link the code I'm writing here when I'm done. Thanks again! – nbubis Dec 15 '15 at 9:00

Triangulation need not any more be adopted once you know Clairaut's Law as operative and geodesic trajectory is found by quadrature for the surfaces of revolution.

EDIT2:

The parametrization you gave is not full, valid for all meridians.

Wlog you could consider

$$X(u,v) = (v \cos u, v \sin u, \psi (v) )$$

$C$ is determined by initial condition: $C_{initial}= v_0 \sin \theta_{0}.$

It is beneficial with respect to easily crossing the parallels in numerical computations to work on arc length basis (primes) using Liouville's theorem for zero geodesic curvature lines aka geodesics.

$$\theta^{'}(s) = -\sin \theta \sin \phi / v$$

where the meridian has slope:

$$\dfrac{dv}{dz} = \tan \phi$$

and polar coordinate rotation:

$$v\,u^{'} = \sin \theta$$

EDIT1:

Say you want to practically compute and construct a great circle on a sphere. Using Runge-Kutta_4 numerically integrate the set of equations to define the circle in space and supply boundary conditions.

$$\theta ^{'} = -\sin \theta \sin \phi / v ;\; \phi^{'} = 1/ (a \cos \phi ) ;\; v^{'} = \sin \phi \cos \theta;\; z^{'} = \cos \phi \cos \theta\; ; u{'}= \sin \theta / v ;$$

(The following is more advanced, after some experience gained in geodesic tracing:)

Geodesics appear to cross all parallels before the one corresponding to $$r_{min} = C, \theta = \pi/2$$

There are however, three ways geodesics move during turn-around at above extreme parallel:

The returning behavior types are commonly traced for:

1. Positive/Negative Gauss curvatures: Geodesics are returning after meeting parallel circle $r_{min}< C_{initial}$ at turn-around.

Non-returning types are for Negative Gauss curvature:

1. Geodesics shoot through after meeting $r_{min}= C$ if $r_{min} > C_{initial}$

2. Geodesics are asymptotic (never reach) $r_{min}= C$ if $r_{min} = C_{initial}$

The latter are demonstrated in WolframAlpha Demo:

ONE_sheet_Hyperboloid_geodesics

• Thanks,I'm not sure what your suggestion is on how to actually draw them. – nbubis Dec 10 '15 at 8:18
• Is it clear now with given example, what CAS are you using? – Narasimham Dec 12 '15 at 6:59
• The effort is much appreciated, but I still don't understand what you are trying to say. – nbubis Dec 13 '15 at 12:29