Computing Lyapunov exponents on a Stiefel manifold

The problem of numerical computation of a few Lyapunov exponents (LEs) of finite-dimensional dynamical systems is considered from the viewpoint of the differential geometry of Stiefel manifolds. Whether one computes one, many or all LEs of a continuous dynamical system by time integration, discrete or continuous orthonormalization is essential for stable numerical integration. A differential-geometric view of continuous orthogonalization suggests that one restricts the linearized vectorfield to a Stiefel manifold. However, the Stiefel manifold is not in general an attracting submanifold of the ambient Euclidean space: it is a constraint manifold with a weak numerical invariant. New numerical algorithms for this problem are then designed which use the fiber-bundle characterization of these manifolds. This framework leads to a new class of systems for continuous orthogonalization which have strong numerical invariance properties and the strong skew-symmetry property. Numerical integration of these new systems with geometric integrators leads to a new class of numerical methods for computing a few LEs which preserve orthonormality to machine accuracy. This idea is also taken a step further by making the Stiefel manifold an attracting invariant manifold in which case standard explicit Runge-Kutta algorithms can be used. This leads to an algorithm which requires only marginally more computation than a standard explicit integration without orthogonalization. These class of methods should be particularly effective for computing a few LEs for large-dimension dynamical systems. The new schemes are straightforward to implement. A test case is presented for illustration, and an example from dynamical systems is presented where a few LEs are computed for an array of coupled oscillators. (C) 2001 Elsevier Science B.V. All rights reserved.