Computing Lyapunov exponents based on the solution expression of the variational system

Abstract

A simple discrete QR algorithm based on a solution expression of the variational equation of a dynamical system is presented for computing the Lyapunov exponents of n-dimensional continuous dynamical systems. The developed numerical scheme of study is based on a time integration using a constant time-step fourth-order Adams–Bashforth method. Numerical results are presented for a Lorenz system with known Lyapunov exponents, and higher dimensional dynamical systems. The algorithm proposed to compute the Lyapunov exponents is found to be robust, computationally efficient, and stable for a sufficiently small step-size h.

Introduction

Lyapunov exponents, the qualitative measures characterizing stability and instability dynamical behaviours of a dynamical system, allow assessment of the sensitivity of solutions of a dynamical system with respect to initial data [17], [18] and are used for measuring the fractal dimension of strange attractors based on the Kaplan–Yorke formula [10]. They provide fundamental information describing behaviour characteristics of the system and understanding of the underlying properties associated with strange attractors [7], [19].

For means of illustration, we assume that the dynamical system under examination is described by the ordinary differential equation system

$$\dot{x}=f(x)\text{,}t>0\text{,}x(0)=x_0\in R^n$$

for a smooth n-vector valued function f(x). Differentiating this system with respect to the initial condition x₀, we derive the n²-dimensional ordinary differential equation system or the variational system

$$\dot{Y}=J(x(t))Y\text{,}Y(0)=I_n\text{,}$$

where I_n denotes the n × n identity matrix and J(x) = ∂f(x)/∂x, the Jacobian matrix of f.

Consider the evolution of an infinitesimal n-parallelepiped [p₁(t), …, p_n(t)] [8], [20] with the axis p_i(t) = Y(t)p_i(0) for i = 1, …, n, where {p₁(0), …, p_n(0)} denotes an orthogonal basis of Rⁿ. The ith Lyapunov exponent, which measures the long-time sensitivity of the flow x(t) with respect to the initial data x₀ at the direction p_i(t), is defined by the expansion rate of the length of the ith axis p_i(t) and is given by

$${\lambda }_i=\underset{t\to \infty }{\lim }\log \frac{\Vert p_i(t){\Vert }_2}{\Vert p_i(0){\Vert }_2}(i=1\text{,}\dots \text{,}n)\text{,}$$

where log denotes the natural logarithm function $\log {\mathrm{e}}^t=t$ and the norm $\parallel x{\parallel }_2=(\sum x_i^2{)}^{1/2}$. Without loss of generality [2], the direction p_i(t) is chosen as the ith column of the identity matrix I_n for i = 1, …, n. Thus if we use the Gram–Schimdt reorthonormalization or QR reorthonormalization (see, for example, [8], [20]) for Y(t) in the form

$$Y(t)=Q(t)R(t)$$

(with uniquely defined orthonormal matrix Q(t) and an upper triangular matrix R(t) with positive diagonal entries R_i,i(t) for i = 1, …, n), we introduce the alternative definition [8]

$${\lambda }_i=\underset{t\to \infty }{\lim }\frac{1}{t}\log R_{i\text{,}i}(t)\text{,}i=1\text{,}\dots \text{,}n\text{.}$$

The existence of this limit or the limit defined in (3) (which is usually assumed to exist in computations) essentially depends on the fundamental theorem of ergodic dynamical systems [16]. Note that, the value ∣λ_i∣ can be very large, and its evaluation fails if we calculate directly the Lyapunov exponents from (5), due to the limitation of computing memory. The numerical computation of Lyapunov exponents are well developed as described by Benettin et al. [3], Wolf et al. [20] and Geist et al. [8], using a discrete QR algorithm. One may also refer to Christiansen and Rugh [6] and Lu et al. [14] for successful developments of continuous algorithms in computing Lyapunov exponents. In this paper, we discuss a new and simple discrete QR method based on the solution expression of (2).