# Computation of a few Lyapunov exponents for continuous and discrete dynamical systems

COMPUTATION OF A FEW LYAPUNOV EXPONENTS FOR CONTINUOUS AND DISCRETE DYNAMICAL SYSTEMS .

by Luca Dieci and Erik S. Van Vleck. Key Words: Lyapunov exponents, weakly skew-symmetric systems, orthogonalization techniques, error analysis.

ABSTRACT.

In this paper, an error analysis of QR based methods for computing the rst few Lyapunov exponents of continuous and discrete dynamical systems is given. Algorithmic developments are discussed. Implementation details, error estimators and testing are also given.

1. INTRODUCTION

In this paper, we consider the computation of a few Lyapunov exponents for continuous and discrete nite dimensional dynamical systems. In DRV2], we considered approximating all of the exponents of continuous dynamical systems, and gave error analysis and convergence results, as well as algorithmic details, for the class of QR based methods. Here we consider the case in which only the rst few exponents are desired. As it turns out, the extension to this case is not automatic, and new theoretical and algorithmic aspects need to be addressed. We still focus on continuous and discrete QR methods, give convergence results, algorithmic details and testing. For continuous dynamical systems, a proper implementation of the continuous QR method appears superior to the discrete QR algorithm, which has thus far been the usual choice. There are many nonlinear systems arising in applications where it is important to know only the rst few (often just the largest) Lyapunov exponents, for example to detect chaotic behavior. It is then natural to study algorithms which target only these few exponents, rather than the entire spectrum of the system. The standard implementations of the QR algorithms we consider will approximate the exponents in decreasing order of magnitude. Naturally, it would also be important to devise techniques to nd the exponents in a given range; for example, the negative and positive exponents closest to the origin. We are currently studying these techniques, and we plan to report on these methods in the future. Work supported in part under NSF Grant DMS-9306412 and NSERC Grant OGP0121873 School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332 U.S.A. Dept. of Mathematical & Computer Sciences, Colorado School of Mines, Golden, CO 80401 U.S.A. 1

For a continuous dynamical system, computation of some Lyapunov exponents by the continuous QR method results in a nonlinear matrix equation which has a weakly skew symmetric structure (see Section 2). Solutions of this type of equations are matrices with orthonormal columns, but the right hand side of the system cannot be written as a skew-symmetric matrix times an orthogonal matrix. As a consequence, we will see that the known automatic unitary integrators, i.e. Gauss Runge-Kutta (GRK) methods (see DRV1]), will not generally preserve orthogonality of the computed solution. Methods for preserving orthogonality in weakly skew systems are discussed in this paper. As it turns out, none of the standard schemes can automatically maintain orthogonality while integrating weakly skew systems (see Section 2). Besides QR based methods, which have been used extensively (see SN], BGGS], WSSV], ER], GSO], GPL], DRV2]), methods based on the Singular Value Decomposition (SV D) have been considered for computation of Lyapunov exponents (see ABK], GK], GPL]). However, these methods appear to be considerably more expensive, they present scaling di culties, and problems arise when singular values coalesce. For these reasons, they will not be considered here. The plan of this paper is as follows. In Section 2, we present the continuous and discrete QR methods for computing a few Lyapunov exponents. We introduce weakly skew-symmetric systems and their discretization. In section 3, we extend the error analysis of DRV2] to the case of a few Lyapunov exponents. It is emphasized that QR based methods are e ective for systems that are regular in the sense of Lyapunov. Section 4 is devoted to implementation considerations. This includes details of how the error in the computation of the exponents is estimated. Section 5 provides numerical results for several nonlinear examples and Sections 6 and 7 give conclusions and references.

2. COMPUTATION OF A FEW LYAPUNOV EXPONENTS Lyapunov Exponents of Discrete and Continuous Dynamical Systems. xi = fi (xi); i = 0; 1; : : : ; x given; xi 2 IRn ; _ x(t) = f (x; t); t 0 ; x(0) = x ; x(t) 2 IRn (t) ;

+1 0 0

We consider nonautonomous discrete and continuous dynamical systems of the following forms (2:1) (2:2)

where the ffi g and f are assumed to be continuously di erentiable. Small perturbations to the orbits fxig and x(t) of (2.1) and (2.2) evolve according to the dynamics of the respective linear variational equations

Yi+1 = Dfi(xi)Yi =: AiYi ; Yi 2 IRn n ; i = 0; 1; : : : ; Y0 = I;

f x x f x x

(2:3)

_ Y (t) = Df (x(t))Y =: A(t)Y (t); Y (t) 2 IRn n ; t 0; Y (0) = I; (2:4) @ where Ai = ( @ i ) i 2 IRn n is henceforth assumed to be full rank, and A(t) = ( @ ) (t) 2 IRn n (t). In @ this paper, we work with (2.3) and (2.4) directly, thus neglecting the errors which presumably went into the computation of fxig and x(t). There are a number of situations where this is legitimate (see DRV2, Remark 2.8]). 2

A0 ; i = 1; : : :, be the fundamental solution of (2.3) (and Y 0 = I ), and let Y (t) be the fundamental solution of (2.4). Then, the following symmetric positive de nite matrices exist

1

DEFINITION 2.1. O], JPS], ER]. Let Y i = Ai?

= ilim (Y iT Y i ) 21i ; !1

and

the logarithms of their eigenvalues are called Lyapunov exponents (LEs for short) of (2.3) and (2.4) respectively, and are denoted by 1 2 n. A Theorem of O] leads to an equivalent way to characterize LEs (see also BGGS], ER]). Let (1) > (2) > be the LEs of (2.3) or (2.4) not repeated by multiplicity. Let E (i) be the subspace of IRn corresponding to the eigenvalues of in (2.5) whose logarithm is less than (i) , so that IRn = E (1) E (2) . Let pk 2 E (k) n E (k+1) , then one has (k) (k) or = tlim 1 log kY (t)pk k ; (2:6) = ilim 1 log kY ipk k ; !1 t !1 i for (2.3) or (2.4) respectively. In (2.6), and throughout this work, k k is always the 2-norm. We are only interested in the rst p (1 p n) most dominant LEs, counted with their multiplicity: 1 2 p . A consequence of the above is that we can avoid nding the full matrices Y 's of (2.3)-(2.4), and it is enough to consider the following rectangular matrix systems

= tlim (Y (t)T Y (t)) 21t ; !1

(2:5)

Yi+1 = AiYi; Yi 2 IRn p ; i = 0; 1; : : : ; Y0T Y0 = I;

_ Y (t) = A(t)Y (t); Y (t) 2 IRn p ; t 0; Y (0) = Y0 ; Y0T Y0 = I :

(2:7) (2:8)

Since Yi and Y (t) are full rank matrices, then the k ; k = 1; ; p, are still de ned analogously to (2.6); thus, see also BGGS], one can in principle obtain the p most dominant LEs. The two standard di culties are: (i) computing Yi (or Y (t)) maintaining linearly independent columns, and (ii) choosing the pk in (2.6) to obtain all of the k ; k = 1; ; p, since any choice of a vector p will unavoidably lead to obtaining only 1 . To deal with these di culties, we need some further results.

LEMMA 2.2. (a) Let (2.7), (2.8), have LEs k ; k = 1; ; p. Let Qi 2 IRn n ; i = 0; 1; , be orthogonal e e ee e and let Ai := QT AiQi 2 IRn n . Then the system Yi = AiYi; Yi 2 IRn p; i = 0; 1; , has same LEs as i e _ (2.7). Likewise, let Q(t) 2 IRn n be di erentiable and orthogonal, and let A := QT AQ ? QT Q 2 IRn n . _ = AY ; Y (t) 2 IRn p has same LEs as (2.8). (b) Now, let Yi = QiZi , Qi 2 IRn p e ee e Then the system Y orthogonal: QT Qi = I , Zi 2 IRp p, i = 0; 1; , and let i

+1 +1

Zi+1 = Bi Zi ; Bi := QT+1 AiQi 2 IRp p : i

(2:9)

Then, (2.9) has same LEs as (2.7). Similarly, let Y = QZ , where Q 2 IRn p is any di erentiable and orthogonal matrix: QT Q = Ip , Z (t) 2 IRp p, and let _ _ Z = B (t)Z ; B := QT AQ ? QT Q 2 IRp p : Then (2.10) has same LEs as (2.8). 3 (2:10)

(or Y (t) = Q(t)Z (t)). Consider the systems (2.9) and (2.10). These systems are called regular ( Ly]) if one has, respectively:

p X k=1 p X k=1 k = ilim !1

e e Proof. Part (a) follows from letting Yi = QiYi (or Y (t) = Q(t)Y (t)), and part (b) from letting Yi = QiZi

1 log(jdet(B i?1 i

0

Zt = tlim 1 trace(B (s))ds = lim sup 1 log j det Z (t)j : k !1 t t

t!1 i?1 X = ilim 1 log j(Bj )kkj; k = 1; : : : ; p; k !1 i

B0 )j) = ilim 1 log(jdet(Zi)j) ; !1 i

If (2.9) and (2.10) are regular and upper triangular, one has ( Ly], SC])

j =0

(2:11) (2:12)

for the LEs of (2.9) and (2.10) respectively. Henceforth, we will assume that our systems are regular. Lemma 2.2, and (2.11) and (2.12), will then provide the basis of the two computational procedures we considered. Discrete QR Algorithm. This is the most widely used technique to approximate LEs (see BGGS], ER], GPL], JPS], SaSa], SN], WSSV]). The algorithm applies to problems of the type (2.7). It must be noticed that this includes any 1-step discretization of (2.8) as well. For example, with hi = ti+1 ? ti, these two 1st order discretizations t into (2.7): Yi+1 = I + hi A(ti)]Yi and Yi+1 = I ? hi A(ti+1)]?1Yi : The idea of the algorithm is the successive computation of QR factorizations. Given orthogonal Y0 : Y0T Y0 = I , let Q0 = Y0 : (2:13) Solve for Zi+1 and then decompose

Zt = tlim 1 Bkk (s)ds ; k = 1; : : : ; p; k !1 t

0

Zi+1 = AiQi; i = 0; 1; : : : Zi+1 = Qi+1Ri+1 ;

(2:14) (2:15)

where Ri+1 is upper triangular with positive diagonal entries. Since we obtain Qi+1Ri+1 = AiQi, then QT+1AiQi = Ri+1 is upper triangular, it is Bi in (2.9). Thus, we can obtain the LEs from (2.11), that is i 1 k = ilim i log((Ri)kk !1

i i?1 1 X log((R ) ) = lim 1 X log((B ) ) : (R1 )kk) = ilim i j kk i!1 i j kk !1 j=1 j =0

(2:16)

REMARKS 2.3.

4

(1) Without loss of generality we may set Y0 = I0p . (2) The QR factorizations can be performed in two di erent ways: (i) so to obtain also the orthogonal complement of Zi+1, that is Zi+1 = Qi+1Ri+1 ; Qi+1 2 IRn n ; Ri+1 2 IRn p, or (ii) simply representing the column space of Zi+1 , that is Zi+1 = Qi+1Ri+1; Qi+1 2 IRn p (QT+1Qi+1 = I ), Ri+1 2 IRp p. We have i adopted this second choice. In the paper of Goldhirsch, Sulem and Orszag GSO] a continuous method for computing a few Lyapunov exponents is presented. In exact arithmetic their method is a continuous QR method. A matrix formulation of the method is given by Geist, Parlitz and Lauterborn GPL]. The basic idea is to avoid directly solving (2.8), by requiring Y (t) = Q(t)R(t), where Q(t) is orthogonal and R(t) is upper triangular with positive diagonal entries, and then solving a di erential equation for Q and determining the upper triangular coe cient matrix for R. As before this can be done in two ways. (i) Here, Q(t) 2 IRn n ; QT Q = QQT = In , and R(t) 2 IRn p is upper triangular. Naturally, upon recovering Q, we also obtain a representation for the orthogonal complement of Y (t). This seems to be an expensive choice since no advantage is taken when nding Q of the fact that typically p << n. (ii) This is a more e cient choice when p << n. We still require Y (t) = Q(t)R(t), but now demand _ Q(t) 2 IRn p ; QT Q = Ip and R(t) 2 IRp p upper triangular. In this case, from the relation Y = _ _ QR + QR = AQR, we obtain _ R = (QT AQ ? S )R =: B (t)R ; and the matrix S is skew-symmetric with _ S := QT Q 2 IRp

p

Continuous QR Algorithm.

(2:17)

8 (QT AQ) ; k > j > kj > < Skj = > 0; k = j ; k; j = 1; : : : ; p: > T :

?(Q AQ)jk ; k < j

(2:18)

_ _ Since R is invertible (because Y is full rank) from the relations RR?1 = QT AQ ? S and Q = AQ ? _ QRR?1 , we now obtain _ Q = (A ? QQT AQQT + QSQT )Q ; or _ Q = H (Q; t)Q ; H (Q; t) = A ? QQT A + QSQT ; (2:19)

where we notice that the matrix H (Q; t) is not necessarily skew-symmetric. The equation (2.19) for the orthogonal factor looks like a quintic in Q, not a cubic as in DRV1-2]. But it can be easily rewritten as the cubic _ Q = AQ ? QQT AQ + QS : (2:20) _ Of course, in case p = n everything reduces to the case treated in DRV1-2], i.e. Q = QS . The writing in (2.19) emphasizes the role of the term QQT , which is also responsible for the lack of skew-symmetry 5

in H , since H T + H = AT (I ? QQT ) + (I ? QQT )A, which is of course 0 if and only if QQT = In . From now on, we only consider this second choice, that is Y (t) = Q(t)R(t); Q(t) 2 IRn p . _ Now, suppose we have solved Q = Q H (t; Q) for Q(t). Observe that B (t) in (2.17) is upper triangular and that Skk = 0; k = 1; ; p. Then, we immediately have _ Rkk = (QT AQ)kkRkk; k = 1; Thus, from (2.12) we can obtain the LEs as

k = tlim !1

; p:

(2:21)

1 log R (t) = lim 1 Z t B (s)ds ; kk t!1 t 0 kk t

k = 1;

;p:

(2:22)

In the error analysis of DRV2] for the continuous QR algorithm, the essential component was maintaining the orthogonality of the computed Q. Motivated by this, we next investigate how to maintain orthogonality of the solution of (2.20) under discretization. We show that this cannot be generally achieved automatically, with the standard formulas, unless p = n, i.e. the square case. Consider equations of this type 8_ < Q = H (Q; t)Q =: F (Q; t); t 0; (2:23) : Q(0) = Q0 ; QT Q0 = I : 0 d These are a special case of conservative systems. Being conservative means dt (QT Q) = 0, and therefore QT Q = C 2 IRp p where C is a constant matrix. For us, QT Q0 = I , and then C = I . If we rewrite (2.23) 0 as the vector system _ q = (I H (Q; t))q =: B (t; q)q F(t; q) ; where is the Kronecker product and

Weakly Skew Systems.

Note that the matrix Q(t) with Q(0) = Y (0)R?1(0) and the resulting upper triangular B (t) are unique. As before, we can assume that Y (0) = Q(0) = I0p .

q = (q ; : : : ; qn ; : : : ; q k ; : : : ; qnk )T ;

11 1 1

then conservative means qT F (t; q) = 0. If we partition Q by rows, i.e. let qT = (QT ; : : : ; QT ), then 1 p

qT F (t; q) = QT H (t; q)Q +

1 1

+ QT H (t; q)Qp = ; p

and we need = 0 for the system to be conservative. The two special cases of interest to us are (i) H is skew, for all t and Q. In this case, obviously = 0. (ii) H = A ? QQT A + QSQT , where S is skew (p p), A 2 IRn n , Q 2 IRn p , n > p, and QT Q = I . In this case, let j = QT H Qj . Then j

j

= 1 QT (A ? QQT A + QSQT + AT ? AT QQT ? QSQT )Qj 2 j 6

so that

j

and thus = 0. Case (i) has been treated previously (see DRV1], C]). Case (ii) is (2.19), or (2.20), and as we saw it arises in the computation of few LEs, which is treated in this work. We propose the following

= 1 ((QT ? QT )A + AT (Qj ? Qj )) = 0 j 2 j

DEFINITION 2.4. The system (2.23) is called weakly skew symmetric if

QT Q = I; QT (H + H T )Q = 0 :

Notice that De nition 2.4 allows for H T 6= ?H . If H T = ?H , the system is just skew-symmetric. The key di erence between cases (i) and (ii) above is that in case (i) = 0 for every vector q 2 IRkn , whereas in case (ii) = 0 only for vectors q exact solutions of the di erential system. This fact precludes application of the results in DRV1], C], SS]. In these works, it was shown that if = 0 for every vector q 2 IRkn, then under discretization with GRK (Gauss-Runge-Kutta) schemes the computed solution also stays orthogonal. However, for case (ii) this is not necessarily true. In fact, GRK schemes, or any symplectic integrator, cannot generally be expected to maintain orthogonality for case (ii). We reason as follows. Consider the general s-stage RK scheme for (2.23):

Qi+1 = Qi + h Qi;j = Qi + h

s X

s X j =1

bj H (Qi;j )Qi;j

ajl H (Qi;l )Qi;l; j = 1; : : : ; s l=1 and assume that Qi is column orthogonal: QT Qi = I . A bit of calculation shows that i s X QT+1Qi+1 = QT Qi + h bj QT (H (Qi;j )T + H (Qi;j ))Qi;j i;j i i j =1 s X ? h2 (bj ajl + bl alj ? bj bl )QT H (Qi;j )T H (Qi;l)Qi;l i;j j;l=1

s X j =1

(2:24)

If the RK scheme is symplectic, the h2 -term drops out (e.g., see SS]) and so QT+1Qi+1 = I if and only if i

bj QT (H (Qi;j )T + H (Qi;j ))Qi;j = 0; i;j

(2:25)

which just cannot be expected to hold, in general. The next example shows this fact clearly. ? ? ?p _ EXAMPLE 2.5. Consider y = Ay, y = y1 , A = 0 ?0 , > 0, y(0) = = p2=2 . Obviously, y2 2=2 ? e t , t 0. Now, orthogonalize y(t), to get Q(t) = (t) , Q(t) 2 IR2 1 . (2.20) in our case is y(t) = e? t kk

y y

_ Q = (I ? QQT )AQ; 7

Q(0) =

:

(2:26)

Consider the symplectic GRK scheme of order 2, the implicit midpoint rule, applied to (2.26). That is, QT + QT Qi+1 = Qi + h I ? Qi +2Qi+1 i 2 i+1 A Qi +2Qi+1 = Qi + hH Qi +2Qi+1 Qi +2Qi+1 :

T T ? ? Q We know from (2.25) that QT+1Qi+1 = 1 , Qi +2Qi+1 (H + H T ) Qi+2 i+1 = 0. Let Qi+1 = s and Qi = . i c To derive a contradiction suppose that QT+1Qi+1 = 1. A bit of algebra gives i

QT+1 + QT 1 i i (H + H T ) Qi + Qi+1 = (s + )2 1 ? 4 (s + )2 ? (c + )2 1 ? 1 (c + )2 ; 2 2 2 4 p and this must be 0. After simpli cation, and using = 2=2, for this right-hand side above we must have

(s ? c)(2 ?

3

? (1 + sc)) = 0:

Since s2 + c2 = 1, then we must have at once s = c and sc = 1=2 which forces s = c = . In other words, the only possible orthogonal solutions are Qi+1 = Qi or Qi+1 = ?Qi. But neither of these has any order, which means that we cannot have QT+1Qi+1 = 1, if Qi+1 is computed by the implicit midpoint rule. i Naturally, the next question is whether or not there are (nonsymplectic) schemes which automatically maintain orthogonality for weakly-skew problems, in particular of the type (ii). None of the standard choices will work. For argument sake, restrict attention to RK formulas with xed coe cients. From (2.24), to preserve orthogonality we must have

s X j =1 s X

bj QT (H (Qij )T + H (Qij ))Qij = 0 ij

(bj ajl + blal ? bj bl )QT H (QT )H (Qil)Qil = 0 ij ij

(2:27)

j;l=1

and i = 0; 1; : : :. But this cannot be achieved for xed values of the RK constants, for all i = 0; 1; : : :, unless H is skew, that is p = n. There exist, however, RK-like schemes with solution-dependent weights which maintain orthogonality. The construction can be based on an idea in Dekker-Verwer (see Example 10.3.8 in DeVe]). and orthogonality preserving. It is a modi the RK tableau 0 1=2 1=2 1 where bj = ej , = z= , z = b

4

EXAMPLE 2.6. The following tableau can be used to obtain a formula of order at least 2, fully explicit

j j j j ?? ?? j

1

cation of the classic 4th order RK tableau (see DeVe]). Consider 0 1=2 0 0 1=2 0 0 0 1

?? ?? ?? ??

b1 b2 b3 b4

4 4 4 3 2

0

X Xe T b bj Fi;j? Fi;j + e FT Fi; and = k ej Fi;j k , where F is the right-hand b i; j j side of the system rewritten in vector form, as after (2.23). From (2.27), to determine the ej we need to b

=2 =1

8

solve the underdetermined system

j T = qT (I (Hi;j + Hi;j ))qi;j ; i;j

0 @1

(Notice that the last 2 equations here would give us a 2nd-order RK method when bj = ej .) Solve this b system, to obtain, e.g., the min 2-norm solution ej , j = 1; 2; 3; 4. Since the j ! 0 with O(h2), then we get b h 2 ej which are only O(h ) apart from those of an order 2 method, and thus we get an order 2 scheme. (If we b let j = 0, j = 1; : : : ; 4, the min 2-norm solution is ej = 1=4, j = 1; : : : ; 4 which gives 2nd order.) b Despite the theoretical possibility of constructing this type of orthogonality preserving schemes, it is di cult for us to see how they could be competitive with projected schemes (see DRV1]). These result from orthogonalizing after each step (say, by a modi ed Gram-Schmidt QR-factorization) the numerical solution of (2.23) computed with any scheme, say a one-step scheme. This is a straightforward procedure, which we have adopted in this work. Clearly, the resulting hybrid scheme maintains orthogonality and inherits the order of the formulas used to integrate (2.23).

1 0e b Be b 1 1 1 A Be @b 0 1=2 1=2 1 e b

1 2 3 4

1 2 3 4

1 0 1 C = @ 0 A: C 1 A

1=2

3. ERROR ANALYSIS

In what follows, this elementary Lemma is useful.

LEMMA 3.1. Let Q ; Q 2 IRn p be di erentiable time dependent orthogonal matrices: QT Q = QT Q = Ip ; 8t. Then, there exists di erentiable orthogonal Q(t) 2 IRn n , such that

1 2 1 1 2 2

Q1 (t) = Q(t)Q2(t); 8t :

and let Qc (t) be the di

1

erentiable orthogonal complement of Q1 (t). Since Q1 and Q2 are di erentiable, their orthogonal complements can be taken so as well. Then, the matrix Q below does the trick

Proof. Based on the following construction. Let (QT )c be the di erentiable orthogonal complement of QT

2 2

Q Q(t) = (Q1; Qc ) (QT2)c : 1 2 REMARK 3.2. If p = n, then Q is simply Q = Q1QT . 2

In order to give the error analysis, we need to extend Theorem 1.6 of DRV2] to the present setting.

T

and such (Ui(t) performs the change of coordinates on the i-th subinterval b guaranteed by Lemma 3.1.) Consider Y (t) 2 IRn p , solution of the problem _ b b bb b Y = AY ; Y (0) = I0p ; Y (ti) = UiT (ti)Y (ti) ; 9

THEOREM 3.3. Let Y (t) = Q(t)R(t). Let Qi (t) be di erentiable and orthogonal on ti ; ti ); i = 0; 1; , and let hm = min(ti ? ti) h > 0. Let Ui (t) 2 IRn n (t) be orthogonal and di erentiable in ti; ti )

+1

that Qi(t) = UiT (t)Q(t) there.

+1

+1

with

b are b b _ Then we have that the LEs of Y Z t the same as those of Y . Moreover, if B (t) := (Qi)T AQi ? (Qi)T Qi is b upper triangular, then k = tlim Bkk (s)ds; k = 1; ; p. !1

0

b b _ A(t) 2 IRn n ; A(t) = UiT (t)A(t)Ui(t) ? UiT (t)Ui(t) ; t 2 ti ; ti+1 ) :

b let Y (t) = Q(t)R(t) be the unique QR decomposition of Y . Then, Y = UiT (t)Q(t)R(t) = Qi(t)R(t) is the b b QR decomposition of Y . Therefore, the R factor of Y and Y is the same and thus the di erential equation b for the R factor of Y is

_ _ b _ R = (QT AQ ? QT Q)R = ((Qi)T AQi ? (Qi)T Qi))R; and the result follows from Lemma 2.2 and (2.22). Finally, we need to realize that we must necessarily truncate the in nite procedures (2.16) and (2.22) to a nite index J and a nite time T , respectively. Thus, for both cases, the best we can hope is to accurately compute truncated LEs

J X 1 X log((R ) ) = 1 J ?1 log((B ) ) and j kk j kk k (J ) = J J j=0 j =1

b b b Proof. Since Y = UiY on each given subinterval, then Y T Y = Y T Y , so that the LEs are the same. Now,

1 Z T B (s)ds = 1 log(R (T )): (3:1) k (T ) = T kk kk T 0

apparent that for given J or T , the truncated LEs di er from the LEs and their relative ordering might be di erent as well. Moreover, if (2.7) and (2.8) are linearizations of (2.1) and (2.2) respectively, then the truncated LEs might also depend on the initial conditions x0. Because of these considerations, determining J and T can be rather tricky. To quantify the error which is made when one works with the truncated LEs of (3.1), rather than with the

j ?1

REMARK 3.4. Clearly, as J or T approach 1, the truncated LEs approach the LEs. However, it is

LEs, we adapt the arguments of DRV2] to the present case. Let

, for the discrete and continuous cases, respectively. Then (see and k (t) = e? k t k (t), k (t) = e 0 Section 2 and (4.1)-(4.3) of DRV2]), for any k > 0 and k > 0, there are constants Kk 1 and Lk 1 so that we have:

Rt Bkk s ds

( )

k(

j ) = e?j

k

k (j ), k (j ) = e l=0

P

log((

Bl )kk )

,

j k k (j ) ?k k (m)j Kk ; j m; and j k k (t) ?k k (s)j Kk ; t s; j k ? k (m) ?k? k (j )j Lk; j m; and j k ? k (s) ?k? k (t)j Lk ; t s :

+ 1 + 1 + 1 + 1

Setting m = 0; s = 0, we get

je? k k j k (j )j Kk ; j 0; and je? k k t k (t)j Kk ; t 0; je? k ? k j k (j )j L? ; j 0; and je? k ? k t k (t)j L? ; t 0; k k

( + ) ( + ) ( ) 1 ( ) 1

10

and so we obtain the following estimates for J 1; T > 0, respectively:

X 1 1 J? ? k ? J log(Lk) J log((Bj )kk) ? j

1 =0

k = k (J ) ? k

k+J k+T

1 log(K ); k

(3:2) (3:3)

1 ? k ? T log(Lk)

1 Z B (s)ds ? = (T ) ? kk k k k T

T

0

1 log(K ): k

In other words, we have that For any given b > 0, there exist nite values of J and T such that

j k (J ) ? k j < b ;

and

j k (T ) ? k j < b ; k = 1; : : : ; p;

(3:4)

for the discrete and continuous cases respectively. How to estimate Kk ; Lk ; k ; k , and in turn choose J and T so to achieve a certain error is discussed in Section 4.

Discrete QR: Few LEs case, discrete or continuous dynamical systems.

Consider the problem (2.7), and the discrete QR algorithm of (2.13)-(2.15) to compute its LEs.

THEOREM 3.5. Suppose that for all i = 1; 2; ; J , for the Zi's of (2.14) we have Zi = Xi + Wi , where Xi is a full rank approximation to Zi, and kWi k < . Then, for su ciently small, the computed LEs bk (J )

satisfy

jbk (J ) ? k (J )j < C ; k = 1; ; p;

where the form of the constant C is shown in the proof below.

Proof. Let Zi = PiUi and Xi = QiRi , where Pi; Qi 2 IRn p and Ui; Ri 2 IRp p denote the unique QR

factorizations of Zi and Xi. By St, p. 516] we have that

kUi ? Rik kWik + b(kRik + kWi k)

with and therefore

b

3kR?1k kWik i 1 ? 2kR?1 k kWik i

1 1 1

( R? kUi ? Rik < (1 + 3 1Ri)2+Rk? ik k ] ) =: C ; ? k i where (Ri) := kRik kR? k is the condition number of Ri . If fZigJ is the exact solution sequence in i i

1

(2.14), then

=1

J J 1X 1X (J ) = J log((Uj?1)kk) = J log((Rj?1)kk + C (i) ) ; k j =1 j =1

11

where jC (i) j C1 ; i = 1;

J . Instead, we are computing

bk (J ) = 1 X log((Rj? )kk) :

J

J j=1

1

Thus, to rst order, we have

j k (J ) ? bk (J )j C

1

J 1X 1 J j=1 (Rj?1)kk k = 1;

;p:

(3:5)

J REMARK 3.6. For J large, we expect that P j =1 (Rj ?1 )kk

<< 1 for itive exponents.

1

J P

j =1 (Rj ?1 )kk

1

>> 1 for

k

< 0, while we expect that

k

> 0. Then, with this method, we should obtain more accurate results for pos-

In the above Theorem 3.5, we can think of as being due to the roundo errors which occur during the matrix multiply of (2.14). Then, would be easily bounded by nkAik EPS (EPS is the machine unit roundo ), which is just as good as one might hope. In essence, then, Theorem 3.5 says that for a given discrete dynamical system (2.7), the discrete QR algorithm delivers accurate approximations of the truncated LEs. But, when the discrete system (2.7) is the result of some discretization of the continuous problem (2.8), naturally we wanted good approximations to the LEs of (2.8). In other words, we would have liked Ai in (2.7) to be the transition matrix of the problem from ti to ti+1, whereas in practice we only have an approximation to it. The question then is whether accurate LEs of a sensible discretization of (2.8) are also accurate LEs of the continuous problem (2.8). We claim that the answer is yes, and in fact that Theorem 3.5 still holds in this case, and that is a bound on the LTE (local truncation error) committed at each step. To see it, consider the problems _ Y i = A(t)Y i ; t 2 ti; ti+1 ] ; Y i (ti) = Qi ; i = 0; 1; ; and let Xi be the computed approximation to Y i?1(ti), that is Xi = AiQi ; i = 0; 1; , where Ai is the discretized transition matrix. Then, Xi = Y i?1(ti)+ Wi , kWi k . If we had Xi = Y i?1(ti), then we would be getting the same LEs as those of the original problem, as a consequence of Theorem 3.3. But then the claim is veri ed, as long as Xi is full rank. This latter fact is obviously achieved for small enough step-size ti ? ti?1.

Continuous QR: Few LEs case, continuous dynamical systems.

Consider now (2.8) and the continuous QR algorithm (2.20)-(2.22) to compute its LEs. We give the error analysis for this algorithm by allowing for the numerical integrations to be done with a variable stepsize. 12

tions to Q(t) at the points 0 = t0 < t1 < : : :, computed with a unitary integrator of order l. Let Qi (t) be the unique orthogonal factor of Z i(t) 2 IRn p : _ Z i(t) = (UiT A(t)Ui) Z i ; Z i(ti) = I0p ; t 2 ti ; ti+1 ] ; i = 0; 1; : : : ; (3:6)

LEMMA 3.7. Assume that A(t) in (2.8) is a C l matrix function. Let Qi; = 0; 1; : : :, be the approxima+1

where Ui 2 IRn n are orthogonal such that UiT Qi = I0p . Suppose that the LTE (local truncation error) in computing Qi is always less then . Then we have

kQi ? Qi(ti )k

+1 +1 +1

:

(3:7)

Proof. Let ti t ti . Let Qi (t) be the solution to _ Qi = (UiT AUi)Qi ? Qi(Qi)T (UiT AUi)Qi + QiS (Qi) ;

_ Qi(ti) = I0p ; S (Qi) = (Qi)T Qi : Let Y i(t) solve

_ Y i = A(t)Y i ; Y i(ti) = Qi :

What we compute over the step from ti to ti+1 is the approximation Qi+1 to this problem _ Qc = AQc ? QcQT AQc + QcSc ; c Qc(ti) = Qi Sc = QT Qc : c _ Now, we have that Z i (t) = Qi (t)Ri(t) and so we have that

Qi (t) = UiT Qc(t) ; t 2 ti ; ti+1 ) ; UiT Qi = I0p

from which it follows that

kQi ? Qi(ti )k

+1 +1

: t 2 ti ; ti+1 ):

With same notation as in Lemma 3.7, consider the ODE

b _ b z = A(t)z ; A(t) = UiT A(t)Ui ;

Let Z i(t) 2 IRn p be the rst p columns of the transition matrix associated to it. Then Qi(t) brings Z i(t) into upper triangular form; in other words:

b b _ B := (Qi)T AQi ? (Qi)T Qi

is upper triangular. By Theorem 3.3,

0

(3:8)

Zt b = tlim 1 Bkk(s)ds ; k !1 t

1 Z T B (s)ds ; k = 1; bkk (T ) = T k

0

;p:

(3:9)

13

b Now, let again be the LTE in integrating for the orthogonal factors Qi, and let B c be the computed b b approximation to B in (3.8). Let the composite trapezoidal rule be used for the integral in (3.9) with B c b replacing B , and let J : tJ = T . Then, what we compute is

c (T ) = k J 1 X hj B c (t ) + B c (t )] : b bkk j T j=1 2 kk j?1

Rewrite and notice that

bc b j bc b j Bkk(tj ) = Bkk (t?) + (Bkk(tj ) ? Bkk (t?)) ;

j

J J 1 X hj (B (t? ) + B (t?)) + X(? h3 )B 00 ( )] ; t j b bkk j?1 bkk j (T ) = T k j ?1 2 12 kk j j =1 j =1

< tj ; j = 1; : : : ; J :

Also, notice that if we have two orthogonal matrices U1 and U2 such that kU1 ? U2 k , then we also have T T T T T (for any bounded matrix M ) kU1 MU1 ? U2 MU2k = k(U1 ? U2 )MU1 + U2 M (U1 ? U2 )k 2 kM k. Then, by summing up the integration errors, and using Lemma 3.7, we get 1 b 00 j k (T ) ? c (T )j T j (? hj )Bkk( j )j + T k 12 j

3 =1

J X

J X j =1

b b hj (kB(tj?1 )k + kB(tj )k) :

(3:10)

Thus, the error between exact and computed truncated LEs is proportional to (i) quadrature errors, and (ii) LTE in integration of Q. We summarize this result in a theorem Let the assumptions of Lemma 3.7 hold, and let the integral in (3.9) be approximated by the composite trapezoidal rule. Then, for the truncated LEs, the error estimate (3.10) holds.

THEOREM 3.8. Consider the continuous QR algorithm (2.20)-(2.22) to approximate the LEs of (2.8).

4. IMPLEMENTATION ISSUES

In this section we give implementation details: on the approximation of LEs, on the estimation of the error due to truncating to a nite time the in nite limits processes, and on the estimation of the error due to discretization. The algorithms will produce both estimates of the Lyapunov exponents and estimates of the error in their computation.

Computing LEs.

The basic algorithms have been discussed in Section 2. For the nonlinear problems (2.1)-(2.2), the solution of the nonlinear equations and the linearization step are carried out simultaneously. For example, for continuous dynamical systems and discrete QR algorithm, we integrate numerically the equations ( y = f (t; y); y(ti) = yi; _ (4:1) _ Z = f (t; y)Z := A(t)Z; Z (ti) = Qi:

y

14

from ti to ti+1 and then decompose (the computed) Z (ti+1) as Qi+1Ri+1 . For continuous dynamical systems and continuous QR algorithm, we end up integrating ( y = f (t; y); y(0) = y0 ; _ (4:2) _ Q = H (Q; t); Q(0) = Q0: The one-step routine RKF45 is used to perform the integrations, and after each integration step the computed value of Q in (4.2) is reorthogonalized by the modi ed Gram-Schmidt procedure (see GvL]). This same procedure is used to perform the decomposition of Zi+1 in (4.1) (see (2.15)). For the continuous algorithm, the integral in (3.1) must be approximated. We use the composite trapezoidal rule. We have also used Richardson's extrapolation to increase the accuracy, but this did not give sizable improvements. In detail, let a(t) := Bkk (t), and let T = tJ . Let hj := tj+1 ? tj , and Ij := h2j a(tj ) + Z tJ J ?1 X X J ?1 hj Ij = a(tj ) + a(tj+1 )], a(tj+1 )], j = 0; ; J . The composite trapezoidal rule is just a(s)ds 0 j =0 2 j =0 J ?1 h3 X j (2) which has error term a ( j ) where j 2 tj ; tj+1 ]. This is a 2nd order rule, globally. To extrapolate, j =0 12 we let b I2j := h2j +2h2j+1 a(t2j ) + a(t2j+2 )]; j = 0; ; J=2; and then obtain a globally fourth order rule by computing 3 b j +1 I I2j := C (h2j ; h2(h )(h2j + )I2j+1 ) ? I2j ; C (h2j ; h2j+1 ) := ((h23j + h23j+1 )) ; j = 0; ; J=2: C 2j ; 2j+1 ? 1 h2j + h2j+1 For both the continuous and discrete QR algorithms we have also tried to accelerate the convergence of the sequence of approximate LEs k;j ; k = 1; ; p; j = 0; ; J by the Aitken 2 process. This did not produce appreciable gains.

Finite Time Error Estimates.

The error analysis of Section 3 allows us to infer that the Lyapunov exponents of the original system will be close to the Lyapunov exponents obtained by our method, in in nite time. What we need to assess is the error due to truncating the in nite time process. In order to do this, we use (3.2) and (3.3), and thus we need to estimate the dichotomy constants. The construction below is used for this purpose, and it is more robust than a similar one in DRV2]. It is enough to discuss the case of a given k, k = 1; ; p.

proximation of the original dynamical system we obtain numerically. Given global dichotomy constants (dichotomy constants that hold on the in nite time interval) we could determine the nite time error for the Lyapunov exponents of the discrete approximation. By combining the nite time error of the Lyapunov exponents for the discrete approximation with the error analysis in section 3 we obtain the error in the nite time approximation of the Lyapunov exponents we obtain numerically as compared with the exact Lyapunov exponents of the original dynamical system. However, we do not have the dichotomy constants for the innite time interval, and we use the information from the nite time interval to approximate these global 15

REMARK 4.1. We use the procedure below to determine the dichotomy constants for the discrete ap-

dichotomy constants. For this reason, we can only say that we are obtaining estimates of the nite time error and not error bounds. If the dichotomy constants obtained from the nite time interval are indicative of those on on the in nite time interval, then our estimates will be indicative of bounds for the nite time error. Let b be the approximate truncated LE k (T ), computed on the (possibly uneven) mesh 0 = t0 < t1 < < tJ , and T = tJ . Let N be a positive integer and let h > 0 be a reference step size (we choose h = 1 for discrete dynamical systems and h approximately equal the average step size for continuous dynamical systems) with Nh T . For i = 1; 2, we set

i) M (j;j+1) := hj ?t<h(j +1)

sup

(?1)iI ( ; t) ; 0 j N ? 1 ; < t T ;

(4:3)

where I ( ; t) is the computed approximation to Bkk (s)ds for continuous dynamical systems, and to

j =t

P log((B ) ) for discrete dynamical systems. Set t j kk

R

i L(i) := M (N)?1;N ) ? (?1)iNhb ;

and

( )

i

:= L : Nh

( ) ( )

i

(4:4)

i

Now, let

( )

i

:= (?1)ib + (i) ; we have

and

i) K (i) := sup M (j;j+1) ?

0

so that for arbitrary t

j N ?1

jh ;

(4:5)

?

(1)

? K? t

(1)

1 I ( ; t) ? b ?t

(2)

+ K? t ;

(2)

(4:6)

which is then used in (3.2) or (3.3). (i) The (K (i) ; (i) ) pairs, i = 1; 2 are optimal with respect to the given value of T . To see this, consider the perturbation (i) := (i) + (i) with corresponding constants K (i) , then for (i) < 0, K (i) ! 1 as t ? ! 1 and for (i) > 0, K (i) K (i) . (ii) This method of computing dichotomy constants also provides us with alternative expressions for the LEs, (i) , which depend on the value of N . Thus, ideally, we could choose N to minimize the di erence between the computed LEs, (i) and b. This reduces to minimizing the absolute value of (i) for i = 1; 2. (iii) In our implementation of the numerical computation of the dichotomy constants we only sample the i) growth/decay factors M (j;j+1) in (4.3) for j = 0; 1; 3; 7; 15; : : : ; N ? 1, where N = 2k for some nonnegative integer k, while modifying the computation of K (i) in (4.5) appropriately. This is much more e cient in terms of both storage and time than sampling for j = 0; 1; 2; 3; 4; : : : ; N ? 1, and from our experience appears to be nearly as accurate. 16

REMARKS 4.2.

To compute the estimate (3.5) is easy, because the norm and condition number of R?1 are easily i obtainable, and so is , either from LTE estimates if we are discretizing a continuous problem, or as roundo error. See the discussion below Remark 3.6. In the estimate of the right-hand side of (3.10), the LTE component is estimated by the local error estimator provided by the numerical integration scheme, kB (tj )k is computed in the 1-norm, and the quadrature error bit is estimated by taking the di erence between the second order approximation of the integral and the fourth order approximation obtained with Richardson's extrapolation. To estimate the total error for the discrete QR method we combine (3.5) and (4.6) to obtain

Computing the Error Estimates (3.5) and (3.10).

j k ? bk (J )j imaxf ;

=1 2

( )

i

J 1X + K (i) =J g + C1 J (R 1 ) j =1 j ?1 kk

(4:7)

and for the continuous QR method we combine (3.10) and (4.6) to obtain

j k ? bk (T )j imaxf ;

=1 2

( )

i

J J X b 1 X h3 b 00 j b + K (i) =T g + T j (? 12 )Bkk( j )j + T hj (kB(tj?1 )k + kB(tj )k) : j =1 j =1

(4:8)

5. NUMERICAL EXPERIMENTS

Here, the algorithms are applied to the computation of LEs for several problems. The numerical results highlight the need of maintaining orthogonality, as well as the overall better performance of the continuous QR method for continuous problems. All computations have been performed in double precision on a Sun SPARCstation with unit roundo EPS 2:2E ? 16. (We use scienti c notation here below, e.g. 2:2E ? 16 = 2:2 10?16.) In the tables, \CPU" refers to CPU times normalized with respect to the least expensive computation.

EXAMPLE 5.1. The rst example we consider is the double rotor map (see KY]). This map models the

position and velocity of two attached rods and is given by

n+1) n+1) _ n+1) _ n+1)

( 1 ( 2 ( 1 ( 2

=M = M_

(n) ( _1n) + 1n) ( ( _2n) 2 ( (n+1) _1n) + c1 sin( 1n+1) ) (n) ( _2 c2 sin( 2 )

where M and M _ are 2 2 matrices and c1 ; c2 are scalars. We use the same values of KY]:

?5 800 ?6: 0 7496 0 1203 M = ?6::602 ?12602 ; M _ = 0::1203 0::8699 ; c1 = 0:3536 ; c2 = 0:5 : :40

(0) (0) (0) (0) For all of our experiments the initial condition is 1 = 2:2, 2 = 3:5, _1 = _2 = 0. In Figure 1 we plot the value of k where N in (4.3)-(4.6) satis es N = 2k versus the computed for value of 2 which shows that 2 decreases as k increases.

17

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0 0 2 4 6 8 10

Figure 1. (k vs.

2

for the rst LE for k = 0; :::; 10)

The rst three rows of Table 1 refer to what we obtained when computing only the rst exponent, while the last four rows refer to the results obtained when computing all four exponents. To compute the error due to discretization for this example we have assumed that the error due to roundo is of size 10?15. The cost of computing all four exponents is approximately four times the cost of computing just the rst exponent. Our results are in good agreement with the values found in KY]. Note in the column \LTE" in Table 1 that the error due to discretization increases as we increase the number of exponents that are computed. The estimate of the nal error, in the column \Error", is predominantly due to the nite time truncation, and suggests that we have obtained at least one signi cant digit for each exponent.

i 1 1 1 1 2 3 4

0.117D+1 0.117D+1 0.117D+1 0.117D+1 0.119D+0 -.488D+0 -.125D+1

i

Table 1. Example 1: K (1) 0.174D+0 0.145D+2 0.109D+0 0.150D+2 0.574D-1 0.153D+2 0.109D+0 0.150D+2 0.712D-1 0.103D+2 0.559D-1 0.105D+2 0.110D+0 0.616D+1

(1)

Discrete QR Method (T = 10; 000) (2) K (2) k 0.153D+0 0.625D+1 8 0.107D+0 0.682D+1 9 0.784D-1 0.767D+1 10 0.107D+0 0.682D+1 9 0.652D-1 0.124D+2 9 0.595D-1 0.107D+2 9 0.100D+0 0.139D+2 9 18

LTE 0.191D-14 0.191D-14 0.191D-14 0.127D-11 0.293D-11 0.522D-11 0.116D-10

Error 0.175D+0 0.111D+0 0.791D-1 0.111D+0 0.723D-1 0.606D-1 0.110D+0

EXAMPLE 5.2. Our next example is the Lorenz equation (see L]) 0 x 1 0 (y ? x) 1 _ @ y_ A = @ x ? xz ? y A : z _ xy ? z

We consider the parameter values = 16, = 4:0 and = 45:92 and the initial condition (x(0); y(0); z(0)) = (0; 1; 0). For this example we compare the discrete and continuous QR methods in terms of accuracy and e ciency. For both methods we integrate using RKF45 with tolerances set to 10?6 and one-step integration mode. The column \Iterates" records the number of time steps taken by RKF45 and \ i(Extrap)" denotes the value of the corresponding exponent we obtain using one step of Richardson extrapolation. From Table 2, we see that the continuous QR methods requires fewer time steps to satisfy the same local error tolerances. From Table 3, we see that not using an orthogonal integrator for the Q equation (in Table 3, MGS stands for Modi ed-Gram-Schmidt) can severely impact performance and accuracy when computing more than one Lyapunov exponent. It must be noticed that the choice of k, such that N = 2k is used in (4.3), has an e ect on the error estimate we get, as it does whether one is approximating all of the exponents or only one of them. There is a trade o to consider, but the basic fact is that the procedure based on the computation of the dichotomy constants can deliver insightful error estimates if one is willing to pay the price. For example, from Table 2, we see that for k = 14 the error estimate suggests that the computed exponent is within ten percent of the exact value. This error is due almost entirely to the error due to nite time truncation as we found the error due to discretization in Table 2 to be approximately 0:5333E ? 4 for the continuous QR method and 0:1655E ? 2 for the discrete QR method. This suggests that a strategy that balances the nite time error more closely with local error tolerances may be more cost e ective. Method Cont QR Cont QR Cont QR Cont QR Cont QR Disc QR Disc QR Table 2. Example 2: Discrete vs. Cont. QR (T = 1000) i Iterates CPU k i 1 0.1492D+01 86425 1 0 1 0.1492D+01 86425 6.30 8 1 0.1492D+01 86425 29.64 10 1 0.1492D+01 86425 145.32 12 1 0.1492D+01 86425 489.21 14 1 0.1501D+01 415273 22.01 8 1 0.1501D+01 415273 115.47 10 Table 3. Example 2: Continuous QR (T = 1000) i Iterates i i (Extrap) 1 0.1492D+01 0.1490D+01 86425 1 0.1489D+01 0.1485D+01 86357 1 0.1490D+01 0.1489D+01 110503 2 0.4767D-02 0.3462D-02 110503 3 -.2249D+02 -.2244D+02 110503 1 0.1546D+01 0.1545D+01 110416 2 -.5388D-01 -.5530D-01 110416 3 -.2250D+02 -.2245D+02 110416 19 Error 0.3308D+01 0.1649D+01 0.3885D+00 0.9641D-01 0.4449D+01 0.2607D+01

Method Cont QR Cont QR (no MGS) Cont QR Cont QR Cont QR Cont QR (no MGS) Cont QR (no MGS) Cont QR (no MGS)

k 10 10 8 8 8 8 8 8

Error 0.1649D+01 0.1644D+01 0.6041D+01 0.3469D+01 0.4370D+01 0.7010D+01 0.3572D+01 0.4874D+01

to the position component of the limit cycle of a van der Pol oscillator. In particular, we consider the following system, similar to that considered in D],

EXAMPLE 5.3. The nal example we consider is a ring of oscillators with an external force proportional

y + (y2 ? 1)y + !2 y = 0 _ 0(xi ? xi?1) ? 0 (xi+1 ? xi)] = y i1 ; i = 1; : : : ; n :

xi + dixi + _

Above, (x) = (x2=2) + (x4=4) is the single well Du ng potential, ; !; ; are scalar parameters, xi is the displacement of the i-th particle, di is the damping coe cient, and we have periodic boundary conditions (x0 = xn and xn+1 = x1). For our experiments we set n = 5 and set = 1, ! = 1:82, = 1 and = 4. We set di = 0:25 for i odd and di = 0:15 for i even. Based on the simulation in D], we expect the largest LE to be positive, and the results of Table 4 con rm this. (The rst LE is in good agreement with the results in D].) The value of K (1) to be used in (4.6) we obtained was nearly identical regardless of whether we used all the information up to T = 2000 or just the information up to T = 200 (for the same k). Instead, we consistently had (1) (T = 2000) < (1) (T = 200). One of the outstanding problems is to achieve good estimates for (i) and K (i) using shorter time intervals. At present, our approach seems to be costly. However, if we are willing to pay the price, we can obtain good error estimates for this problem as well; see the error estimates for k = 10; 12; 14 in Table 4. Table 4. Example 3: Continuous QR (1) i T K (1) CPU k Error i 1 0.1007D-01 200 1 0 2 0.4223D-02 200 1 0 1 0.2386D-02 2000 9.66 0 2 0.1028D-02 2000 9.66 0 1 0.1007D-01 200 0.5374D+00 0.1500D+01 1.18 6 0.6666D+00 2 0.4223D-02 200 0.3773D+00 0.4797D+00 1.18 6 0.6099D+00 1 0.2386D-02 2000 0.5269D+00 0.1500D+01 11.21 6 0.6696D+00 2 0.1028D-01 2000 0.4250D+00 0.6382D+00 11.21 6 0.6098D+00 1 0.1007D-01 200 0.5646D-01 0.2385D+01 3.97 10 0.8190D-01 2 0.4223D-02 200 0.3723D-01 0.1098D+01 3.97 10 0.4380D-01 1 0.2386D-02 2000 0.4867D-01 0.2385D+01 38.78 10 0.8160D-01 2 0.1028D-02 2000 0.3404D-01 0.1098D+01 38.78 10 0.4056D-01 1 0.1007D-01 200 0.1858D-01 0.2436D+01 16.89 12 0.3084D-01 2 0.4223D-02 200 0.1016D-01 0.1107D+01 16.89 12 0.1574D-01 1 0.2386D-02 2000 0.1080D-01 0.2436D+01 172.97 12 0.1924D-01 2 0.1028D-02 2000 0.6967D-02 0.1107D+01 172.97 12 0.1096D-01 1 0.2386D-02 2000 0.2739D-02 0.2769D+01 736.49 14 0.4166D-02 2 0.1028D-02 2000 0.2000D-02 0.1208D+01 736.49 14 0.2938D-02

6. CONCLUSIONS

We have provided a uni ed error analysis for discrete and continuous QR based methods for approximating the rst few LEs of discrete and continuous dynamical systems. We have also discussed and implemented 20

ways to estimate the errors in the computed approximations. Based on our results, it appears that QR based methods are a numerically sound and e ective way to approximate LEs for systems of the type (2.7),(2.8). For continuous dynamical systems, our preference goes to the continuous QR method, implemented as described in this work. There are a number of important issues which remain to be addressed. It would be desirable (a) to approximate the LEs in a more selective way; (b) to obtain sharp estimates of the error due to the nite time truncation of the in nite process; (c) to adaptively decide where to truncate the in nite time process based upon the error estimates; (d) to extend the error analysis to include the full nonlinear problem, as our analysis really deals only with the linear case. We are working on these issues, but at present have no complete answers for them.

7. REFERENCES

ABK] Abarbanel, H. D. I., Brown R. and Kennel M. B.: \Local Lyapunov Exponents Computed from Observed Data," Journal of Nonlinear Science 2 (1992) pp. 343-366. BGGS] Benettin, G., Galgani, L., Giorgilli, A. and Strelcyn, J.-M.: \Lyapunov Exponents for Smooth Dynamical Systems and for Hamiltonian Systems; A Method for Computing All of Them. Part 1: Theory ", and \ : : : Part 2: Numerical Applications ", Meccanica 15 (1980), pp. 9-20, 21-30. C] Cooper, G. J.: \Stability of Runge Kutta Methods for Trajectory Problems," IMA J. Numer. Anal. 7 (1987), pp. 1-13. DeVe] K. Dekker, J.G. Verwer, Stability of Runge-Kutta methods for sti nonlinear di erential equations, CWI Monograph Vol. 2, North-Holland (1984). DRV1] L. Dieci, R. D. Russell, E. S. Van Vleck, \Unitary Integrators and Applications to Continuous Orthonormalization Techniques," to appear in SIAM J. Numer. Anal. 30, December 1993. DRV2] L. Dieci, R. D. Russell, E. S. Van Vleck, \On the Computation of Lyapunov Exponents for Continuous Dynamical Systems," submitted to SIAM J. Numer. Anal. (1993). D] U. Dressler, \Symmetry property of the Lyapunov spectra of a class of dissipative dynamical systems with viscous damping," Physical Review A, 38-4 (1988), pp. 2103-2109. ER] Eckmann, J.-P. and Ruelle D.: \Ergodic Theory of Chaos and Strange Attractors," Reviews of Modern Physics 57, (1985) pp. 617-656. GPL] Geist, K., Parlitz, U. and Lauterborn, W.: \Comparison of Di erent Methods for Computing Lyapunov Exponents," Prog. Theor. Phys. 83 (1990), pp. 875-893. GSO] Goldhirsch, I., Sulem, P.-L. and Orszag, S. A.: \Stability and Lyapunov Stability of Dynamical Systems: A Di erential Approach and a Numerical Method," Physica 27D, (1987) pp. 311-337. GK] Greene, J. M. and Kim, J-S.: \The Calculation of Lyapunov Spectra," Physica 24D (1987), pp. 213-225. 21

JPS] Johnson, R.A., Palmer, K.J. and Sell, G.: \Ergodic Properties of Linear Dynamical Systems," SIAM J. Math. Anal. 18 (1987), pp. 1-33. KY] Kostelich, E. J. and Yorke, J. A.: \Lorenz Cross Sections of the Chaotic Attractor of the Double Rotor," Physica 24D (1987) pp. 263-278. Ly] Lyapunov A.: Problem General de la Stabilite du Mouvement, Annals of Mathematics Studies 17, Princeton University Press (1949). O] Oseledec, V. I.: \A Multiplicative Ergodic Theorem. Lyapunov Characteristic Exponents for Dynamical Systems," Trudy Moskov. Mat. Obsc. 19, (1968) pp. 197-231. SaSa] Sano, M. and Sawada, Y.: \Measurement of the Lyapunov Spectrum from a Chaotic Time Series," Phys. Rev. Let. 55 (1985) pp. 1082-1085. SC] Sansone, G. and Conti, R.: Nonlinear Di erential Equations, Pergamon Press, New York (1964). SS] Sanz-Serna, J. M.: \Runge-Kutta Schemes for Hamiltonian Systems," BIT, 28 (1988), pp. 877-883. St] Stewart, G. W.: \Perturbation Bounds for the QR Factorization of a Matrix," SIAM J. Numer. Anal. 14 (1977) pp. 509-518. SN] Shimada, I. and Nagashima, T.: \A Numerical Approach to Ergodic Problem of Dissipative Dynamical Systems," Prog. Theor. Phys. 61 (1979) pp. 1605-1616. WSSV] Wolf, A. , Swift, J. B., Swinney H. L. and Vastano, J. A.: \Determining Lyapunov Exponents from a Time Series," Physica 16D (1985) pp. 285-317.

22