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?

  • $\begingroup$ 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? $\endgroup$ – Andrew D. Hwang Dec 12 '15 at 16:03
  • $\begingroup$ @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. $\endgroup$ – 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.

  • $\begingroup$ 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? $\endgroup$ – nbubis Dec 9 '15 at 15:59
  • $\begingroup$ Any additional thoughts? $\endgroup$ – nbubis Dec 13 '15 at 12:29
  • $\begingroup$ I've added some details about possible methods I thought about. $\endgroup$ – levap Dec 14 '15 at 8:36
  • $\begingroup$ Well done! I'll try and link the code I'm writing here when I'm done. Thanks again! $\endgroup$ – 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.


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 $$


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:


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.