Lyle Noakes
email: lyle.noakes@uwa.edu.au

Tudor S Ratiu
email: tudor.ratiu@epfl.ch

BI-JACOBI FIELDS AND RIEMANNIAN CUBICS FOR LEFT-INVARIANT SO(3)

Keywords: Lie group, Riemannian manifold, Jacobi field, trajectory planning, mechanical system, rigid body, nonlinear optimal control, Riemannian cubic

Bi-Jacobi fields are generalized Jacobi fields, and are used to efficiently compute approximations to Riemannian cubic splines in a Riemannian manifold M. Calculating bi-Jacobi fields is straightforward when M is a symmetric space such as bi-invariant SO(3), but not for Lie groups whose Riemannian metric is only left-invariant. Because left-invariant Riemannian metrics occur naturally in applications, there is also a need to calculate bi-Jacobi fields in such cases. The present paper investigates bi-Jacobi fields for left-invariant Riemannian metrics on SO(3), reducing calculations to quadratures of Jacobi fields. Then left-Lie reductions are used to give an easily implemented numerical method for calculating bi-Jacobi fields along geodesics in SO(3), and an example is given of a nearly geodesic approximate Riemannian cubic.

Introduction

For T > 0 the velocity field x (1) of a C 1 curve x : [0,T ] → M in a C ∞ m-manifold M is defined by x (1) (t) := d dt x(t). Take M to be a Riemannian manifold, with Riemannian metric , and associated Levi-Civita covariant derivative ∇. The Jacobi operator on C ∞ vector fields X along a geodesic γ : [0,T ] → M is given by H(X) := ∇ 2 t X + R X,γ (1) γ (1) where R is Riemannian curvature. A vector field X along γ is a Jacobi field when HX = 0. As is well-known [START_REF] Postnikov | The Variational Theory of Geodesics[END_REF], Jacobi fields comprise a 2mdimensional real vector space of fields along γ, and are the variational fields of variations of γ through geodesics.

Whereas geodesics are critical points of the energy integral T 0

x (1) ,x (1) dt, Riemannian cubics are critical points of the higher order functional J(x) := T 0 ∇ t x (1) ,∇ t x (1) dt where x : [0,T ] → M and its velocity field x (1) are prescribed at t = 0 and t = T . Equivalently ( [START_REF] Gabriel | Spline interpolation in curved space[END_REF][START_REF] Noakes | Cubic splines on curved spaces[END_REF][START_REF] Noakes | Null cubics and Lie quadratics[END_REF]), a Riemannian cubic is a C ∞ curve satisfying ∇ 3 t x (1) + R ∇ t x (1) ,x (1) x (1) = 0.

(1.1)

This Euler-Lagrange equation for Riemannian cubics resembles the Jacobi equation HX = 0 for X = ∇ t x (1) , except that now X is defined along the Riemannian cubic x which is usually not a geodesic (geodesics comprise a codimension 2m class of Riemannian cubics). Most studies of Riemannian cubics focus on cases where M is highly symmetrical, especially S 3 and bi-invariant SO(3) ( [START_REF] Noakes | Cubic splines on curved spaces[END_REF][START_REF] Noakes | Null cubics and Lie quadratics[END_REF][START_REF] Noakes | Non-null Lie quadratics in E 3[END_REF][START_REF] Noakes | Duality and Riemannian Cubics[END_REF][START_REF] Noakes | Lax constraints in semisimple Lie groups[END_REF][START_REF] Noakes | Asymptotics of null Lie quadratics in E 3[END_REF][START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF][START_REF] Ramamoorthi | Fast construction of accurate quaternion splines[END_REF][START_REF] Camarinha | On the geometry of Riemannian cubic polynomials[END_REF][START_REF] Crouch | The dynamic interpolation problem: on Riemannian manifolds, Lie groups, and symmetric spaces[END_REF]). Geodesics are by far the simplest kind of Riemannian cubic, and the easiest to compute. On the other hand all curves, including Riemannian cubics, are nearly geodesic over sufficiently small time intervals. So in §3 of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] it turns out to be useful to study variations of geodesics through Riemannian cubics, and our starting point is the following consequence of Lemma 3 in [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF].

Proposition 1.1. Let Y be the variational field of a variation through Riemannian cubics of the geodesic γ. Then Y is a bi-Jacobi field, namely H 2 Y = 0. This explains our interest in bi-Jacobi fields, as infinitesimal variations of geodesics through Riemannian cubics. The bi-Jacobi fields comprise a 4m-dimensional vector space containing the Jacobi fields.

In §4 of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] it is shown how bi-Jacobi fields can be used to efficiently calculate piecewise nearly cubic interpolants in M , and when M is a symmetric space (as in §5 and §6 of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF]) the bi-Jacobi fields are easily found. However it is not so easy to find the bi-Jacobi fields in other cases of interest, such as SO(3) with a Riemannian metric that is only left-invariant. The present paper focuses on this case, because of its interest in applications [START_REF] Bloch | Nonholonomic Mechanics and Control[END_REF][START_REF] Holm | Geometric Mechanics and Symmetry: From Finite to Infinite Dimensions[END_REF][START_REF] Noakes | Geometry for robot path planning[END_REF]. The paper proceeds as follows.

In § 2 M is taken to be a general m-dimensional left-invariant Lie group G, and the Jacobi operator H is left-Lie reduced to a composite I -1 • K • L. Here I is the inertia transform corresponding to the left-invariant Riemannian metric, and K,L are first order linear operators on curves in the Lie algebra g of G. More precisely, the operator L is a Lax operator and K := L • I -ad(IV ) where V is the left-Lie reduction of the velocity field γ (1) . Then ω = Ad(γ)IV is constant, generalizing inertial angular momentum from the case of left-invariant SO(3) (the well-known theorem from classical mechanics stating that the spatial angular momentum is conserved).

There are m + 1 elementary linearly independent Jacobi fields along γ, leaving another 3m -1 linearly independent bi-Jacobi fields to be found. We find that bi-Jacobi fields X correspond to solutions X,W 1 ,W 2 ,W 3 : [0,T ] → g of the 4m-dimensional system of ODEs:

L X = W 1 , KW 1 = IW 2 , LW 2 = W 3 , KW 3 = 0. (1.2)
Nonhomogeneous Lax equations of the form LY = W are solvable by quadrature. So finding bi-Jacobi fields reduces via (1.2) to solving differential equations for Y : [0,T ] → g of the form KY = W , where W : [0,T ] → g is given. This reduces (Theorem 2.9) to a quadrature involving the fundamental solution Φ of a nonautonomous homogeneous linear ODE on the orthogonal complement J of the kernel K of ad(ω). Then (Corollary 2.10), Φ is computed in terms of Jacobi fields, and substitution for Φ in Theorem 2.9 gives the bi-Jacobi fields as quadratures in terms of Jacobi fields. So computation of bi-Jacobi fields reduces to quadratures in terms of m -1 nonelementary Jacobi fields. Now there are formulae in terms of elliptic functions for geodesics in left-invariant SO(3) [START_REF] Whittaker | A Treatise on the Analystical Dynamics of Particles and Rigid Bodies[END_REF], and so the Jacobi fields (variational fields of γ through geodesics) are computable by differentiation of families of geodesics with respect to parameters. It therefore follows from Corollary 2.11 that the bi-Jacobi fields along geodesics in left-invariant SO(3) are given by quadratures in terms of elliptic functions (we thank Michael Pauley for pointing this out).

Unfortunately the resulting formulae seem much too complicated for efficient implementation in a practical algorithm. So in § 4 we calculate bi-Jacobi fields along geodesics in left-invariant SO(3) using numerical solutions of ODEs instead of quadratures. The first step is to find a basis {X i : 1 ≤ i ≤ 6} of Jacobi fields, as follows. The Lie reduction V of the elementary Jacobi field γ (1) is determined as a numerical solution of (2.4), and then γ is a numerical solution of γ (1) (t) = γ(t)V (t). The other three elementary Jacobi fields are expressed simply in terms of γ, and Lie reductions of nonelementary Jacobi fields X 5 ,X 6 are particular numerical solutions of K • L X = 0, where Lemma 4.1 is used to compute K • L.

The basis of Jacobi fields then extends to a basis of bi-Jacobi fields by numerically solving K • L Xi+6 = I Xi for 1 ≤ i ≤ 6. In Example 4.1 the basis of bi-Jacobi fields is calculated in some 0.25 seconds on a 1.7GHz MacBook Air with 4GB of memory. In § 5 these computations are applied in the framework of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] to approximate Riemannian cubics near a geodesic γ : [0,T ] → SO [START_REF] Crouch | The dynamic interpolation problem: on Riemannian manifolds, Lie groups, and symmetric spaces[END_REF]. Given a bi-Jacobi field X we have a curve x X : [0,T ] → SO(3). At least when X is suitably small, x X approximately satisfies the ODE (1.1) for a Riemannian cubic. In Example 5.1, X is only moderately small, and x X is roughly similar but not extremely close to γ. The initial and terminal values of x X and its derivative are also interpolated by a chart-based interpolant x ch , with very different results, as illustrated in Figure 5.2. Although x ch agrees with x X to first order at t = 0 and t = 5, x X appears to take a more roundabout route. An approximation to a Riemannian cubic should normally produce a smaller value of J, and indeed J(x X ) = 1.45 < J(x ch ) = 5.25.

Some theoretical results

For any

C k -curve χ on an interval [a,b] ⊂ R, we use the notation χ (k) = d k dt k χ.
Let G be a connected Lie group of finite dimension m with Lie algebra g. Definition 2.1. The left-Lie reduction of a vector field X defined along a

C ∞ curve x : [a,b] → G is the curve X : [a,b] → g given by X(t) := dL(γ(t)) -1
e X(t) where (momentarily) L stands for left-multiplication, and e ∈ G is the identity.

The left-Lie reduction of the velocity vector field x (1)

along x is denoted by V : [a,b] → g. The Lax operator along x is the first order linear differential operator L on C ∞ curves Y : [a,b] → g given by (LY )(t) := Y (1) (t) + ad(V (t))(Y (t)) = Ad(x(t)) -1 d dt (Ad(x(t))Y (t)). ( 2
Y (t) = Ad(x(t)) -1 C + t 0 Ad(x(s))W (s) ds , where C ∈ g. (2.2)
Suppose 1 G has a bi-invariant Riemannian metric , BI , and also a left-invariant Riemannian metric , . The inertia transform is the positive-definite , BI -self-adjoint linear automorphism I of g given by v,I(w) BI := v,w for any v,w ∈ g. For the Levi-Civita covariant derivative ∇ of , , the left-Lie reduction of ∇ v w, with v,w left-invariant vector fields on G, is

1 2 ([v, w] -h(v, w)). (2.3)
where (for now) h : g × g → g is symmetric bilinear and given by h(v,w)

:= I -1 ([Iv,w] + [Iw,v])
for any v,w ∈ g. Formula 2.3 is found as in §3.2 of [START_REF] Michor | On the geometry of the Virasoro-Bott group[END_REF], with a sign change due to the fact that the calculations in [START_REF] Michor | On the geometry of the Virasoro-Bott group[END_REF] are for right-invariant Riemannian metrics. The extension of h in (2.3) to left invariant vector fields is denoted by the same letter and is defined by h(v, w)(x) := h dL(x) -1 e v(x),dL(x) -1 e w(x) . Consequently, denoting as usual,

∇ t := ∂ ∂t + ∇ x (1)
, the covariant derivative along the curve x, we obtain the following.

Lemma 2.2. The left-Lie reduction of ∇ t x (1) is V (1) -I -1 [IV,V ].
As in [START_REF] Michor | On the geometry of the Virasoro-Bott group[END_REF] §3.3 (but for left-invariant metrics), from (2.3) one also gets the curvature formula.

Lemma 2.3. The left-Lie reduction of the curvature is

R(u,v)w = - 1 4 [[u,v],w] + [u,h(v,w)] -[v,h(u,w)] + h(u,[v,w]) -h(v,[u,w] -h(u,h(v,w)) + h(v,h(u,w)) -2h([u,v],w)).
Take x to be a geodesic γ : [0,T ] → G, i.e., ∇ t γ (1) = 0. By Lemma 2.2, this left-reduces to the Euler equation

V (1) (t) = I -1 [IV (t),V (t)]. (2.4)
By Equation (2.1), d dt (Ad(γ(t))IV (t)) = Ad(γ(t))(IV (1) (t) + [V (t),IV (t)]) = 0 by (2.4).

So ω := Ad(γ(t))IV (t) is constant. If ω = 0 then γ is constant, which is of little interest. Suppose from now on that ω = 0. Consider the second order differential operator H := ∇ 2 t + S t defined on fields along γ, where S t is the self-adjoint linear endomorphism given by S t (X) := R(X,γ (1) )γ (1) . Lemma 2.4. For any vector field X defined along γ, the left-Lie reduction of HX is

I -1 • K • L X where K is the first order linear differential operator on C ∞ curves in g given by K := L • I -ad(IV ). Proof. By (2.3) the left-Lie reduction of ∇ t X is X(1) + 1 2 [V, X] - 1 2 h(V, X). So by (2.3) again, the left-Lie reduction of ∇ 2 t X is X(2) + 1 2 [V (1) , X] + 1 2 [V, X(1) ] - 1 2 h(V (1) , X) - 1 2 h(V, X(1) ) + 1 2 [V, X(1) ] + 1 4 [V,[V, X]] - 1 4 [V,h(V, X)] - 1 2 h(V, X (1)
)

- 1 4 h(V,[V, X]) + 1 4 h(V,h(V, X)) = X(2) + 1 2 [V (1) , X] + 1 2 [V, X(1) ] - 1 2 h(V (1) , X) - 1 2 h(V, X(1) ) + 1 2 [V, X(1) ] + 1 4 [V,[V, X]] - 1 4 [V,h(V, X)] - 1 2 h(V, X(1) ) - 1 4 h(V,[V, X]) + 1 4 h(V,h(V, X)) = X(2) + 1 4 [h(V,V ), X] + 1 2 [V, X(1) ] - 1 4 h(h(V,V ), X) - 1 2 h(V, X(1) ) + 1 2 [V, X(1) ] + 1 4 [V,[V, X]] - 1 4 [V,h(V, X)] - 1 2 h(V, X (1)
)

- 1 4 h(V,[V, X]) + 1 4 h(V,h(V, X)),
where (2.4) is used to substitute for V (1) in the first line. Combining this with Lemma 2.3, the left-Lie reduction of HX simplifies eventually to

X(2) + 1 2 [h(V,V ), X] + [V, X(1) ] -h(V, X(1) ) -h(V,[V, X]) =I -1 ( X(2) + I[I -1 [IV,V ], X] + I[V, X(1) ] -[IV, X(1) ] -[I X(1) ,V ] -[IV,[V, X]] -[I[V, X],V ]) =I -1 (L • I -ad(IV )) • L( X).
The operator H is also calculated in §3.4, §3.5 of [START_REF] Michor | On the geometry of the Virasoro-Bott group[END_REF] (for right-invariant Riemannian metrics and without factorisation). As we see next, the factorisation in Lemma 2.4 reduces computation of Jacobi fields to solving nonhomogeneous linear differential equations of the form

KY = W (2.5)
for W : [0,T ] → g.

A Jacobi field is a vector field X along the geodesic γ for which HX = 0. Equivalently, by the factorisation in Lemma 2.4, L X = W where KW = 0. So the Jacobi fields along γ comprise a 2m-dimensional real vector space and, by Equation (2.2), can be found by quadratures in terms of elements of ker(K). (Something similar is done for bi-Jacobi fields in Lemma 2.7.)

The following elementary examples span an (m + 1)-dimensional subspace of Jacobi fields.

Example 2.1. Corresponding to Y = 0 ∈ ker(K) we have X(t) = Ad(γ(t)) -1 C where C ∈ g, namely X(t) = dR(γ(t)) e C where R(x) denotes right-multiplication by x.

Example 2.2. Y = V (1) ∈ ker(K). Taking X to be the velocity field γ (1) defined along γ we have X = V and L X = Y . Proposition 2.5. For any Y ∈ ker(K), Y (t),V (t) and Y (t),IV (t) are independent of t.

Proof. Since I is self-adjoint with respect to the bi-invariant inner product, So, by Example 2.2, g e ,g m : g → R given by g e (v) := I -1 [Iv,v],v and g m (v) := [Iv,v],v are integrals of (2.4). Definition 2.6. A bi-Jacobi field is a vector field X along γ for which H 2 X = 0.

d dt Y (t),V (t) = d dt Y (t),IV (t) BI = IY (1) (t),V (t) BI + Y (t),IV (
Bi-Jacobi fields comprise a 4m-dimensional real vector space including the Jacobi fields.

Lemma 2.7. Computation of bi-Jacobi fields reduces by quadratures to solving equations of the form (2.5).

Proof.

By Lemma 2.4, X is a bi-Jacobi field when there are

C ∞ functions W 1 ,W 2 ,W 3 : ]0,T ] → g for which L X = W 1 , KW 1 = IW 2 , LW 2 = W 3 , KW 3 = 0.
The lemma follows, because the nonhomogeneous Lax equation is solved by (2.2). Definition 2.8. Denote by J the , BI -orthogonal complement in g of the kernel K of ad(ω) : g → g. Let r be the nullity of ad(ω). For any v ∈ g write

v K + v J where v K ∈ K and v J ∈ J . For t ∈ [0,T ] set • f (t) := π • Ad(γ(t)) • I -1 • Ad(γ(t)) -1 : g → J ,
where π : g → J is , BIorthogonal projection; • g(t) := ad(ω) • f (t) • ι : J → J where ι is the inclusion of J in g. Denote by Φ : [0,T ] → End R (J ) the solution of the nonautonomous homogeneous linear ODE Φ (1) (t) = -Φ(t) • g(t) with Φ(0) = 1 J .

(2.6)

Theorem 2.9. The general solution of (2.5) is

Y (t) = I -1 Ad(γ(t)) -1 κ + t 0 (Ad(γ(s))W (s)) K ds + Φ(t) -1 ρ + t 0 Φ(s)h(s) ds ,
where κ ∈ K, ρ ∈ J , h : [0,T ] → J and

h(t) = ad(ω) • f (t) κ + t 0 (Ad(γ(s))W (s)) K ds + (Ad(γ(t))W (t)) J .
Proof.

Decompose Z(t) := Ad(γ(t))IY (t) into its K-and J -orthogonal components, namely

Z(t) = Z K (t) + Z J (t). If Y satisfies (2.5) then Z (1) (t) = ad(ω) • f (t)Z(t) + Ad(γ(t))W (t). So, Z K (t) = Z K (0) + t 0 (Ad(γ(s))W (s)) K ds, and
Z (1) J (t) = g(t)Z J (t) + h(t) (2.7)
where

h(t) : = ad(ω) • f (t)Z K (t) + (Ad(γ(t))W (t)) J = ad(ω) • f (t) κ + t 0 (Ad(γ(s))W (s)) K ds + (Ad(γ(t))W (t)) J
on substitution for Z K with κ := Z K (0). The general solution of (2.7) is hence

Z J (t) = Φ(t) -1 Z J (0) + t 0 Φ(s)h(s) ds . (2.8)
Therefore, if ρ := Z J (0), we get

Y (t) = I -1 Ad(γ(t)) -1 κ + t 0 (Ad(γ(s))W (s)) K ds + Φ(t) -1 ρ + t 0 Φ(s)h(s) ds .
Corollary 2.10. For ρ ∈ J , we have Φ(t) -1 ρ = Ad(γ(t))I(L Xρ )(t) where X ρ is the Jacobi field along γ satisfying X ρ (0) = 0 and (L Xρ )(0) = I -1 Ad(γ(0)) -1 ρ.

Proof. By Theorem 2.9,

Y ∈ ker(K) ⇐⇒ Ad(γ(t))IY (t) = κ + Φ(t) -1 ρ + t 0 Φ(s)ad(ω) • f (s)κ ds ,
where κ ∈ K and ρ ∈ J . So, for any ρ ∈ J , define Y ρ ∈ ker(K) by Ad(γ(t))IY ρ (t) := Φ(t) -1 ρ. Define a field X ρ along γ by Xρ (t) := Ad(γ(t)) -1 t 0 Ad(γ(s))Y ρ (s) ds. By (2.2) and Lemma 2.4, X ρ is a Jacobi field, indeed the unique Jacobi field satisfying the conditions in the corollary.

Corollary 2.11. The bi-Jacobi fields along the geodesic γ : [0,T ] → G are given by quadratures in terms of γ and its Jacobi fields.

Proof. By Corollary 2.10, Φ is given by quadratures in γ and the Jacobi fields. Therefore so are the bi-Jacobi fields, by Theorem 2.9.

The Euler-Poincaré formalism

The Euler-Poincaré formalism for G-invariant variational problems consists in replacing a G-invariant Lagrangian L on T G by a reduced Lagrangian defined on g, then studying a reduced variational problem for curves in g. This formalism extends to discretisations of the variational problem, as follows.

The Lagrangian, L, is approximated by a real-valued function on G × G, integration is replaced by summation, and there is a discrete Euler-Lagrange algorithm based on discrete Euler-Poincaré equations for sequences in G [START_REF] Marsden | Discrete Euler-Poincaré and Lie-Poisson equations[END_REF]. The discrete Euler-Lagrange equations generate a flow that is symplectic with respect to a discrete symplectic form on G × G constructed from the approximation to L. The discrete reduced Lagrangian is a real valued function defined on G, and the associated discrete Euler-Poincaré equations preserve a Poisson structure on G induced from the discrete symplectic form, namely solutions to the discretised problem are structure-preserving in the same way as solutions to the original continuous variational problem. This has advantages for numerical approximations of solutions to the original problem, especially for long-term behaviour where conservation laws are more or less enforced.

The Euler-Poincaré formalism extends to higher order Lagrangians [START_REF] Gay-Balmaz | Invariant higher-order variational problems[END_REF][START_REF] Gay-Balmaz | Invariant higher-order variational problems II[END_REF] where Riemannian cubics are of central importance as examples. These papers focus on a variational problem for a reduced Lagrangian, rather than reductions of solutions of the Euler-Lagrange equation in [START_REF] Noakes | Cubic splines on curved spaces[END_REF][START_REF] Noakes | Null cubics and Lie quadratics[END_REF][START_REF] Noakes | Non-null Lie quadratics in E 3[END_REF][START_REF] Noakes | Duality and Riemannian Cubics[END_REF][START_REF] Noakes | Lax constraints in semisimple Lie groups[END_REF][START_REF] Noakes | Asymptotics of null Lie quadratics in E 3[END_REF][START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] and in the present paper. Although at some level the two points of view are equivalent, each has advantages depending on the context.

In the context of [START_REF] Noakes | Non-null Lie quadratics in E 3[END_REF][START_REF] Noakes | Asymptotics of null Lie quadratics in E 3[END_REF], it seems likely that structure-preserving numerical methods would be useful in studying asymptotic properties of Riemannian cubics. Considering the successes of the Euler-Poincaré formalism in constructing such methods for first-order variational problems, it seems likely that the Euler-Poincaré formalism might also be useful for long-term approximations of Riemannian cubics. This topic has not been explored in the literature so far as we are aware, but [START_REF] Gay-Balmaz | Invariant higher-order variational problems[END_REF][START_REF] Gay-Balmaz | Invariant higher-order variational problems II[END_REF] provide a starting point.

On the other hand the aim of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF], and of the present paper, is to provide short-term approximations of Riemannian cubics that are sufficiently explicit to satisfy boundary conditions and interior point conditions for interpolation problems arising in engineering applications. The simplest example is the boundary value problem described in § 1 with stepsize T sufficiently small for the Riemannian cubic x to be nearly geodesic. More generally, we may want to approximate near-geodesic natural cubic splines as in [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] for arbitrary left-invariant Riemannian metrics.

It seems unlikely that such explicit approximations would be useful for studying asymptotics of Riemannian cubics, or even for refined local estimates. Indeed, the difficulty of satisfying boundary and interior point conditions is such that a trade-off has to be made between accuracy of approximation and efficiency of computation, as seen already in [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF] for the case of bi-invariant Riemannian metrics. Since explicit approximations are only needed and valid over intervals of moderate stepsize, their long-term dynamics are of little importance. Considering the trade-off, standard numerical techniques for initial-value problems are perfectly fit for purpose in calculating approximations of arcs of Riemannian cubics, in the same way as in [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF].

The present situation is somewhat2 analogous to approximating Euclidean elastic splines by polynomial splines, where the calculation of polynomial segments is not a particular issue once the polynomials are found, and the advantage of polynomials over genuine elastica is due to their explicit nature which leads to manageable computations.

So, in the present context, the need is for sufficiently explicit approximations of arcs of Riemannian cubics with sufficient accuracy to construct approximations to Riemannian cubic splines. The present paper achieves this, enabling computations of splines in G along the lines of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF]. We shall see, in § 4, a practical method to compute the bi-Jacobi fields X needed to construct our approximations x X . Then examples are given in § 5 for G = SO(3), illustrating the quality of approximation by the x X to Riemannian cubics. K • LW =IW (2) + I[V,W (1) ] + [V,IW (1) ] -[IV,W (1) ]

Computations of Bi-Jacobi fields in left-invariant Lie groups

+ I[I -1 [IV,V ],W ] + [V,I[V,W ]] -[IV,[V,W ]].

Proof.

Substituting Y = LW = W (1)

+ [V,W ] in KY = L • IY -[IV,Y ] = IY (1) + [V,IY ] -[IV,Y ], we get K • LW = I(W (1) + [V,W ]) (1) + [V,I(W (1) + [V,W ])] -[IV,W (1) + [V,W ]]
= IW (2) + I[V (1) ,W ] + I[V,W (1) ] + [V,IW (1) ] + [V,I[V,W ]] -[IV,W (1) ] -[IV,[V,W ]] and the lemma follows by substituting for V (1) using (2.4).

Example 4.1. Consider Euclidean 3-space E 3 as a Lie algebra with respect to the crossproduct. The adjoint action on E 3 defines an isomorphism of Lie algebras ad : E 3 → so(3) onto the Lie algebra so(3) of the rotation group G = SO(3), relative to the usual commutator bracket of matrices. The Euclidean inner product , E corresponds under ad to a bi-invariant inner product on so(3), which extends, by left-translation, to a bi-invariant Riemannian metric , BI on SO(3). Let Ĩ be the positive-definite self-adjoint endomorphism of E 3 whose matrix with respect to the standard basis is diagonal with entries 1,1/2,1/3. Let I be the endomorphism of so(3) corresponding to Ĩ under ad. As in § 2, I defines a left-invariant Riemannian metric , on SO(3), given by v,w = v,I(w) BI for v,w ∈ so [START_REF] Crouch | The dynamic interpolation problem: on Riemannian manifolds, Lie groups, and symmetric spaces[END_REF].

Although quadrature formulae are available, it was more convenient to use numerical methods to approximate the , -geodesic γ : [0,5] → SO(3) with γ(0) the identity and γ (1) (0) = ad(1/3,1,5/3). Mathematica's NDSolve was used to solve the nonlinear ODE (2.4) for V , and again to solve the first order linear system γ (1) (t) = γ(t)V (t) with variable coefficients. We have ω = ad(1/3,1/2,5/9).

As in Example 2.1 there are three Jacobi fields with left-Lie reductions Xi = Ad(γ(t)) -1 ad(e i ) ∈ so(3), with {e 1 ,e 2 ,e 3 } the standard basis of E 3 . Figure 4.1 shows the images of the ad -1 • Xi for i = 1,2,3, with labels near t = 0. Figure 4.2 shows the image of ad -1 X4 , where X4 = V , as in Example 2.2.

Besides these elementary examples we see, first in Figure 4.2, images of ad -1 X5 ,ad -1 X6 : [0,5] where X 5 ,X 6 are the Jacobi fields found from numerical solu-

of K • L X = 0, with initial conditions X5 (0) = [ω,V (0)], X (1) 5
(0) = [ω,V (1) (0)], X6 (0) = [ω, X5 (0)], X (1) 6
(0) = [ω, X (1) 5 (0)].
Then {X 1 ,X 2 ,X 3 ,X 4 ,X 5 ,X 6 } is a basis of Jacobi fields along γ.

This extends to a basis {X i : 1 ≤ i ≤ 12} of bi-Jacobi fields by solving K • L Xi+6 = I Xi with initial conditions Xi (0) = X(1) i (0) = 0 for 1 ≤ i ≤ 6. Figures 4.3 and 4.4 show the images of ad -1 Xi+6 for 1 ≤ i ≤ 3 and 4 ≤ i ≤ 6 respectively (with labels at t = 5 since these curves all begin at 0).

Applications to Riemannian cubics in left-invariant SO(3)

From now on take G = SO(3) with a left-invariant Riemannian metric , corresponding to a positive-definite diagonal 3 × 3 real matrix Ĩ, as in Example 4.1. Following Example 8 of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF], define E : T SO(3) → SO(3) by E(z,w) := z(1 + v/2)(1 -v/2) -1 where z ∈ SO(3), w ∈ T SO(3) and v := z -1 w. Let γ : [0,T ] → SO(3) be a geodesic with respect to the left-invariant Riemannian metric. Definition 5.1.

For any field X along a geodesic γ define x X : [0,T ] → SO(3) by x X (t) := E(γ(t),X(t)).

By Lemma 4 of [START_REF] Noakes | Approximating near-geodesic natural cubic splines[END_REF], for some constant L depending only on γ, we have

∇ 3 t x (1) X + R(∇ t x (1)
X ,x

X )x (1) X x X (t) ≤ L X 2 2,∞ , (1)
namely the ODE (1.1) for a Riemannian cubic is satisfied by x X with O( X 2 2,∞ ) errors. For any C ∞ curve x : [0,T ] → SO(3), Lemma 2.2 gives

J(x) = T 0 Ṽ (1) (t) -Ĩ-1 (( Ĩ Ṽ (t)) × Ĩ Ṽ (1) (t)), Ĩ Ṽ (1) (t) -( ĨV (t)) × Ṽ (t) E dt (5.1)
which should be small (or at least critical) when x is a Riemannian cubic.

Example 5.1.

In Example 4.1, let X be the Jacobi field along γ : [0,5] → SO(3), given by X = 12 i=1 a i Xi where a i = (i -6) × 10 -2 . Figure 4.4 shows t → x X (t)p r ∈ E 3 (blue) and t → γ(t)p r ∈ E 3 (red), where p r = (1,2,3)/ √ 14 ∈ S 2 is a reference point. The endpoints are labelled as 0 and 5. So x X appears roughly similar but not extremely close to γ. We find X (5) are ad -1 (0.312966,1.01025,1.66543) and ad -1 (0.204442,1.17616,2.24975), respectively. Chart based methods can also be used to interpolate this data as follows. Let φ : U → φ(U ) be a diffeomorphism, where U is open in SO(3) and φ(U ) is open in E 3 . Suppose that x X (0),x X (5) ∈ U , and let x E : [0,5] → E 3 be the cubic polynomial satisfying

x E (0) = φ • x X (0),x (1)
E (0) = dφ x X (0) x (1)
X (0),x E (5) = φ • x X (5),x Figure 5.2 shows the images of t → x X (t)p r (blue) and t → x ch (t)p r (red), labelled 0 and 5 at endpoints. The curves are strikingly different, and the image of x ch is apparently much shorter with a sharp turn near t = 5. Closer examination shows that the velocities at t = 0 do indeed agree. We find J(x X ) = 1.45407 and J(x ch ) = 5.24754. So, as should happen, by this measure the nearly cubic x X is much better than the chart-based interpolant.

Conclusion

Bi-Jacobi fields along geodesics in left-invariant Lie groups are shown to be computable by quadratures from Jacobi fields. A practical method of computing bi-Jacobi fields is given in terms of numerical solutions of ODEs for left-Lie reductions, and special attention is given to the case of left-invariant SO(3). Building on this and results of a previous paper, bi-Jacobi fields are used to construct approximations to nearly geodesic Riemannian cubics.

Figures 4. 1

 1

  , 4.2 show reductions of the Jacobi fields as curves in so(3) ∼ = E 3 and Figures 4.3, 4.4 show reductions of the 6 remaining bi-Jacobi fields.

. 1 )

 1 Given any W : [a,b] → g, the solutions Y of the nonhomogeneous Lax equation LY = W are found from

1 2 3 Fig. 4 . 1 .

 1341 Fig. 4.1. Transformed Reductions of Jacobi Fields X 1 ,X 2 ,X 3 in Example 4.1

5 6 Fig. 4 . 2 . 8 9Fig. 4 . 3 .

 5642843 Fig. 4.2. Transformed Reductions of Jacobi Fields X 4 ,X 5 ,X 6 in Example 4.1

Fig. 4 . 4 .

 44 Fig. 4.4. Transformed Reductions of Bi-Jacobi Fields X 10 ,X 11 ,X 12 in Example 4.1

Fig. 5 . 1 . 1 x

 511 Fig. 5.1. Images of t → x X (t)pr and t → γ(t)pr in Example 5.1

5 Fig. 5 . 2 .

 552 Fig. 5.2. Images of t → x X (t)pr and t → x ch (t)pr in Example 5.1

  BI = IY(1) (t),IV (t) BI + IY (t),IV(1) (t) BI = [IY (t),V (t)],IV (t) BI + IY (t),[IV (t),V (t)] BI = 0.

	d dt	Y (t),IV (t) =	d dt	IY (t),IV (t)

1) (t) BI = [IV (t),Y (t)],V (t) BI + Y (t),[IV (t),V (t)] BI = 0

by (2.4), bi-invariance, and because KY = 0. Similarly,

More generally G might have only a bi-invariant pseudo-Riemannian metric. For instance G might be semisimple, with bi-invariant pseudometric generated by the negative of the Killing form. The analysis proceeds similarly, but I is not necessarily positive-definite.

Although Euclidean natural cubic polynomial splines already satisfy a variational condition, the elastic variational condition is different, more difficult to satisfy, and preferable for some applications.

Acknowledgement. We gratefully acknowledge a careful and critical reading by the referee, including instructive and interesting remarks on Euler-Poincaré reduction and numerical integrators.

Partially supported by NCCR SwissMAP and grant 200021-140238, both of the Swiss National Science Foundation, and by the government grant of the Russian Federation for support of research projects implemented by leading scientists, Lomonosov Moscow State University under the agreement No. 11.G34.31.0054.