Solving the geodesics on the ellipsoid as a boundary value problem

Abstract

The geodesic between two given points on an ellipsoid is determined as a numerical solution of a boundary value problem. The secondorder ordinary differential equation of the geodesic is formulated by means of the Euler-Lagrange equation of the calculus of variations. Using Taylor’s theorem, the boundary value problem with Dirichlet conditions at the end points is replaced by an initial value problem with Dirichlet and Neumann conditions. The Neumann condition is determined iteratively by solving a system of four first-order differential equations with numerical integration. Once the correct Neumann value has been computed, the solution of the boundary value problem is also obtained. Using a special case of the Euler-Lagrange equation, the Clairaut equation is verified and the Clairaut constant is precisely determined. The azimuth at any point along the geodesic is computed by a simple formula. The geodesic distance between two points, as a definite integral, is computed by numerical integration. The numerical tests are validated by comparison to Vincenty’s inverse formulas.