Finding geodesics joining given points

Finding a geodesic joining two given points in a complete path-connected Riemannian manifold requires much more effort than determining a geodesic from initial data. This is because it is much harder to solve boundary value problems than initial value problems. Shooting methods attempt to solve boundary value problems by solving a sequence of initial value problems, and usually need a good initial guess to succeed. The present paper finds a geodesic gamma : [0, 1] -> M on the Riemannian manifold M with gamma(0) = x(0) and gamma(1) = x(1) by dividing the interval [0,1] into several sub-intervals, preferably just enough to enable a good initial guess for the boundary value problem on each subinterval. Then a geodesic joining consecutive endpoints (local junctions) is found by single shooting. Our algorithm then adjusts the junctions, either (1) by minimizing the total squared norm of the differences between associated geodesic velocities using Riemannian gradient descent, or (2) by solving a nonlinear system of equations using Newton's method. Our algorithm is compared with the known leapfrog algorithm by numerical experiments on a 2-dimensional ellipsoid Ell(2) and on a left-invariant 3-dimensional special orthogonal group SO(3). We find Newton's method (2) converges much faster than leapfrog when more junctions are needed, and that a good initial guess can be found for (2) by starting with Riemannian gradient descent method (1).