Abstract In this paper, an alternative method to compute the Lyapunov exponents of dynamical systems described by Ordinary Di¤erential Equations (ODEs) is introduced. The Lyapunov exponents are computed in terms of the solutions of two piecewise linear ODEs that approximate, respectively, the solutions of the original ODE and its associated variational equation. This approach is strongly connected with the Local Linearization method for ODEs and its major advantage is that these piecewise linear ODEs might be exactly integrated in a non simultaneous way. The performance of the method is illustrated with a numerical example.

Key words and phrases. Lyapunov Exponents, Local Linear Approximation, nonlinear dynamical systems.

1

Introduction

The increasing interest in the computation of Lyapunov exponents [1] is motivated by its fundamental role in the theory of dynamical systems. Lyapunov exponents allow for the generalization of the linear stability analysis from perturbations of steady state solutions to perturbations of time-dependent solutions, and also provide a meaningful way to characterize the asymptotic behavior of nonlinear dynamics [2]. In particular, in the last two decades the Lyapunov exponents have been widely used as a major tool to identify chaotic behaviors [3]. ¤

Departamento de Matematica Aplicada, Centro de Neurociencias de Cuba, Apartado 6880, La Habana, Cuba y Instituto de Cibernetica, Matematica y Fisica, Calle 15, e/ C y D, Vedado, La Habana 4, C.P. 10400, Cuba z Departamento de Matematica Aplicada, Facultad de Matematica y Computacion, Universidad de la Habana, San Lazaro y L, La Habana, C.P. 10400, Cuba

1

Whereas the exact computation of Lyapunov exponents is only possible for a few simple cases, the use of accurate and e¢cient numerical methods becomes very important. This paper concerns only with the computation of Lyapunov exponents for continuous time dynamical systems generated by ordinary di¤erential equations (ODEs). For this purpose, two types of numerical methods have been considered [4]: the QR decomposition based methods and the SVD decomposition based methods. They have a common main strategy, namely to perform some matrix decomposition of the solution of the associated variational equation. Wether the decomposition is performed after or before the numerical integration of the variational equation, the method is called discrete or continuous. In any case, the computation of Lyapunov exponents by these methods involves the simultaneous integration on a …nite interval of time of the given ODE and additional di¤erential equation-namely, the variational equation or the equation describing the trajectories of the matrix factors. In this paper, an alternative approach is introduced that can be described as follows. First, the given nonlinear ODE and its associated variational equation are approximated by two piecewise linear ODEs. Second, a QR or SVD decomposition based method (either discrete or continuous) is applied to these two new equations to obtain approximations to the Lyapunov exponents of the original system. A major advantage of this method is that it does not require the simultaneous integration of the above mentioned piecewise linear ODEs. In the present paper we focus on this approach in case of discrete decomposition based methods. Since the piecewise linear approximations involved are closely related with the Local Linearization (LL) method for the numerical integration of ODEs, section 2 gives a brief description of the LL method. Also, it contains some basic results concerning to Lyapunov exponents, which are used in the sequel. In section 3 the new method for the computation of Lyapunov exponents is introduced. Section 4 gives bounds for the approximation error of the Lyapunov exponents. Such bounds take into consideration both the error due to the truncation of the in…nite time interval in the de…nition of Lyapunov exponents and the error due to the piecewise linearization of the original ODE and its associated variational equation. Finally, section 4 illustrates the performance of an implementation of the method that includes a adaptive strategy to control the error due to truncation time.

2 2.1

Preliminaries Local Linearization method

Let :

x= f(t; x); t ¸ t0 2 R; x(t0) = x0 2 Rm 2

(1)

be an m-dimensional ODE that satis…es su¢cient conditions for the existence and uniqueness of solutions [5]. In addition, suppose that the function f : R £ Rm ! Rm is di¤erentiable. The Local Linearization (LL) method [6]-[9], also called exponentially …tted Euler method [11, 12] or Euler exponential method [13, 14], is a numerical method for the integration of ODEs. It is derived from the local linearization of the right member of (1) with respect to both variables x ; t; and the exact computation of the resulting linear ODEs. Speci…cally, at each point tn , the approximate solution of (1) given by the LL method is de…ned by the recursive relation [15]: ytn+1 = ytn + Á(tn; ytn ; hn );

(2)

where tn = t0 + hn (0 · hn · h; n = 0; 1; 2; :::; ) are the discretization times, h < 1 is the maximum step size and yt0 = x0 is the given initial condition. Here, Á(s; u; hn ) = R0 (fx0 (s; u); hn )(f(s; u) + hn ft0 (s; u)) ¡ R1 (fx0 (s; u); hn )ft0 (s; u) and Rk (M; h) =

Z

h

exp(Mu)uk du 0

for all s 2 R, u 2 R and M 2 R : fx0 (s; u), ft0 (s; u) denote, respectively, the derivatives of the function f with respect to x and t evaluated at the point (s; u). See [11],[15]-[17] for di¤erent ways to compute the expression (2). At any point t ¸ t0 ; the Local Linear Approximation of the solution of (1) is de…ned as ([18]): m

m£m

y(t) = ytnt + Á(tnt ; ytnt ; t ¡ tnt );

(3)

sup kx(t) ¡ y(t)k · Mh2 :

(4)

where nt = max fn = 0; 1; 2; ::: : tn · tg. This continuous time approximation will be useful in obtaining bounds on the accuracy of the approximation of Lyapunov exponents. It can be shown [18] that LL approximation (3) converges uniformly to the solution of (1) with order 2, i.e., for any T < 1 there exists a constant M > 0 such that t2[t0 ;T ]

Furthermore, (3) is the solution of the (piecewise) linear ODE :

y = f(tn ; ytn ) + fx0 (tn ; ytn )(y(t) ¡ ytn ) +f 0t (tn ; ytn )(t ¡ tn ) y(tn ) = ytn ; where t ¸ t0 ; t 2 [tn ; tn+1 ]. 3

(5)

2.2

Lyapunov exponents for ODEs

Let :

X= fx0 (t; x(t))X;

(6)

X(t0 ) = I

be the variational equation for a given trajectory x of (1), whose solution X is called a fundamental solution matrix. De…nition 1 ([2]) The Lyapunov exponents with respect to the trajectory x of (1) are: ° ° 1 ¸ix = log(¹ix ) = lim log °X(t)pix ° ; t!1 t where ¹ix (i = 1; :::; r) are the distinct eigenvalues of the matrix ¤X = lim (X(t)t X(t))1=2t; t!1

and pix are the corresponding eigenvectors.

In terms of the QR factorization of X; the Lyapunov exponents can be expressed as ([19]): Z t 1 1 i e x (s)]ii ds; ¸x = lim log[Rx (t)]ii = lim [A (7) t!1 t t!1 t ¡ t0 t 0

where Qx (t) and Rx (t) are, respectively, orthogonal and upper triangular matrie x (t) = Qx (t)| fx0 (t; x(t))Qx (t): ces such that X(t) = Qx (t)Rx (t) and A In practice, the expression (7) can be computed only on a …nite interval of time. This leads to the following de…nition.

De…nition 2 ([20]) On the basis of the QR factorization, the truncate Lyapunov exponents with respect to the trajectory x(t) of (1) are:

¸ix (T )

1 1 = log[Rx (T )]ii = T ¡ t0 T ¡ t0

Z

T

t0

f (s)]ii ds; [A x

(8)

where T < 1. Through this paper we will assume that the system (1) has exponential dichotomy [20] so the following lemma holds. Lemma 3 (see lemma 4.2 in [19]) Given positive numbers »i ; ½i ; there exist constants Ki ; Li ¸ 1 (i = 1; 2; :::; r); not depending on T , such that ¡½i ¡

1 1 log(Li ) · ¸ix (T ) ¡ ¸ix · »i + log(Ki ): T ¡ t0 T ¡ t0 4

3

Approximate computation of the Lyapunov Exponents via Local Linearization

Let x be a trajectory of the ODE (1) and ¸ix the i-th Lyapunov Exponents of this trajectory. We de…ne the approximation to the i-th Lyapunov exponent ¸ix as ¸iy (T ) =

° ° 1 log °Y(T )piy ° ; T

(9)

where piy denotes the i-th eigenvector of the matrix ¤Y = (Y(T )t Y(T ))1=2T and Y(T ) is the solution of the piecewise linear equation :

Y = fx0 (tn ; y(tn ))Y; for t 2 [tn ; tn+1 ] Y(tn ) = Yn

(10)

at the instant of time T < 1. Here, t0 · tn < tn+1 · T with n = 1; ::; N , Y(t0 ) = I and y is the solution of the piecewise linear equation (5), which is the LL approximation to the solution of (1). The expression (9) can be computed in a similar way as ¸ix via the expression (7), i.e., by rewriting (9) in terms of the QR or SVD decomposition of the solution of the equation (10). The values tn and T can be …xed in advance or adaptively selected by strategies that control the error in the numerical integration of the equation (1) by the LL scheme (2) and the truncation time error due to the integration on a …nite interval of time. It is worth to note that the piecewise linear approximating equations (10) and (5) may be used to set a di¤erential equation for the factors of the QR or SVD decomposition of Y(t), so leading to a continuous decomposition method (following the terminology in [19]), though this will not be pursued in the present paper.

4

Accuracy of the approximate Lyapunov Exponents

In this section, the relation between ¸iy (T ) and ¸ix will be studied. With this aim the relation between X(t) and Y(t); as well as that of ¸iy (T ) and ¸ix (T ), are also analyzed. Theorem 4 Suppose that the following conditions are satis…ed: for all t; u; v there are positive constants k0 ; k1 ; l1 such that kf(t; u)k · k0 ;

(11)

kfx0 (t; u)k · k1 ;

(12)

5

kfx0 (t; u) ¡ fx0 (t; v)k · l1 ku ¡ vk :

(13)

Let X(t) and Y(t) be, respectively, the solutions of the equations (6) and (10). Then sup kX(t) ¡ Y(t)k · Ch; t2[t0 ;T ]

where the constant C does not depend on h. The following three lemmas will be used in the proof of this theorem. Lemma 5 (Gronwall’s Lemma, [5]) Let u; v; c ¸ R t0 be positive functions de…ned on [t0; t] : If c is di¤erentiable and v(t) · c(t) + t0 u(s)v(s)ds then

.

Z

v(t) · c(t0 ) exp(

t

t0

u(s)ds) +

Z

t

t0

0

·

c (s) exp

µZ

t

u(¿ )d¿ s

¶¸

ds:

(14)

Lemma 6 Let A1 ; B1; C1; A2; B2 ; C2 be matrices such that the products involved in the following inequalities are de…ned. Then: kA1 B1 ¡ A2 B2 k · kA1 k kB1 ¡ B2 k + kB2 k kA1 ¡ A2 k ;

(15)

kA1B1 C1 ¡ A2 B2 C2 k · kA1 ¡ A2 k kB2 k kC2 k + kB1 ¡ B2 k kA1 k kC2k + kC1 ¡ C2 k kA1 k kB1k : (16)

Proof. These inequalities follow, respectively, from the equalities : A1B1 ¡A2 B2 = A1 (B1 ¡B2 ) + (A1 ¡A2 )B2 ; A1 B1 C1¡A2 B2 C2 = A1 B1(C1¡C2) + A1(B1 ¡B2)C2+(A1¡A2 )B2 C2 :

Lemma 7 It holds that

6

sup t2[tn ;tn+1 ]

kfx0 (t; x(t)) ¡ fx0 (tn ; y (tn ))k · l1(k0 + M )h;

(17)

where M is the constant de…ned in (4), and l1; k0 are the constants in the assumptions of Theorem 8. Proof. From the condition (13) it follows that sup t2[tn ;tn+1 ]

kfx0 (t; x(t)) ¡ fx0 (tn ; y (tn ))k · l1

t2[tn ;tn+1 ]

· l1

t2[tn ;tn+1 ]

kx(t) ¡ y(tn )k

sup

kx(t) ¡ x(tn )k

sup

+l1 kx(tn ) ¡ y(tn )k : Considering the equation (1) in its integral form and using the condition (11) we obtain that sup t2[tn ;tn+1 ]

kx(t) ¡ x(tn )k · k0 hn · k0 h:

From this, the inequality (4) and considering that h < 1; the assertion of the lemma is obtained. Proof. of Theorem 4. From equations (6) and (10) it is obtained that Z t X(t) = X(tnt ) + F(s; X(s))ds tnt

and Y(t) = Y(tnt ) +

Z

t

G(s; Y(s))ds tnt

for all t 2 [t0 ; T ], where F(t; X) = fx0 (t; x(t))X(t) and G(t; Y) = fx0 (tnt ; ytnt )Y(t). Using these equalities in a recursive way we obtain X(t) = X(t0 ) +

n t ¡1 Z tn+1 X n=0

F(s; X)ds +

tn

Z

t

F(s; X)ds tnt

and Y(t) = Y(t0 ) +

n t ¡1 Z tn+1 X n=0

G(s; Y)ds +

tn

Therefore, 7

Z

t

G(s; Y)ds: tnt

kX(t) ¡ Y(t)k ·

n t ¡1 Z tn+1 X tn

n=0

+

Z

kF(s; X) ¡ G(s; Y)k ds

t

tnt

kF(s; X) ¡ G(s; Y)k ds:

(18)

Using (15) we obtain ° ° kF(t; X) ¡ G(t; Y)k · °fx0 (t; x(t)) ¡ fx0 (tnt ; ytnt )° kXk ° ° + °f 0 (tnt ; ytn )° kX ¡ Yk ; x

t

From (17) and (12) it follows that

kF(t; X) ¡ G(t; Y)k · ®0¯0 h + k1 kX(s) ¡ Y(s)k ; where ¯0 = sup kX(t)k and ®0 = l1 (k0 + M): t2[t0 ;T ]

Substituting the last inequality in (18) we have

kX(t) ¡ Y(t)k ·

n t ¡1 Z tn+1 X tn

n=0

+

Z

(®0 ¯0 h + k1 kX(s) ¡ Y(s)k)ds

t

tnt

(®0 ¯0 h + k1 kX(s) ¡ Y(s)k)ds

· ®0 ¯0 h(t ¡ t0 ) + k1

Z

t t0

kX(s) ¡ Y(s)k ds:

By the Gronwall’s Inequality (14) we have

kX(t) ¡ Y(t)k · ®0¯0 h =

Z

t t0

exp (k1 (t ¡ s)) ds

®0¯0 h (exp(k1 (t ¡ t0 )) ¡ 1): k1

This implies that sup kX(t) ¡ Y(t)k · Ch; with C =

t2[t0 ;T ]

l1(k0 + M )¯0 (exp(k1 (T ¡ t0 )) ¡ 1); k1

which is the assertion of the Theorem. The next theorem states a relation between the truncated Lyapunov exponents and ¸iy (T ).

¸ix (T )

8

Theorem 8 Under conditions of the Theorem 4 and Lemma 3 there exist positive constants h0 ; b and ai (not depending on h and ai = ai (T ) converging to zero as T tends to in…nite) such that ¯ i ¯ ¯¸x (T ) ¡ ¸iy (T )¯ · bh

and

¯ i ¯ ¯¸x ¡ ¸iy (T )¯ · ai + bh

hold for all h < h0 and i = 1; ::; r: The following lemma will be used to prove this theorem. Lemma 9 ( Theorem 3.1 in [21]) Let A = QR and A + E = QE RE be, respectively, the QR decompositions of the square matrices A and A + E. If kA¡1k kEk < 12 then QE = Q + W and RE = R + F, where kWk ·

3 kA¡1k kEk : 1 ¡ 2 kA¡1 k kEk

Proof. of Theorem 8 In terms of the QR factorization of Y; the expression (9) can be rewritten as Z T 1 1 i e y (s)]ii ds; ¸y (T ) = log[Ry (t)]ii = [A (19) T ¡ t0 T ¡ t0 t0 e y (t) = Qy (t)| fx0 (tnt ; ytn )Qy (t) and Y(t) = Qy (t)Ry (t): where A t Using this expression and (8) we have ¯ i ¯ ¯¸x (T ) ¡ ¸iy (T )¯ ·

N ¡1 ¯ ¯ 1 X ¯ f ¯ f h sup ¯[Ax (s)]ii ¡ [Ay (s)]ii ¯ : T ¡ t0 n=0 s2[tn ;tn+1 ]

From (16) it follows that

¯ ¯ ° ° ¯ f fy (s)]ii ¯¯ · ° fx (s) ¡ A fy (s)° ¯[Ax (s)]ii ¡ [A °A °

° ° · kQx (s) ¡ Qy (s)k °fx0 (tnt ; ytnt )° kQy (s)k + kQx (s) ¡ Qy (s)k kQx (s)k kfx0 (s; x(s))k ° ° + °fx0 (s; x(s)) ¡ fx0 (tnt ; ytnt )° kQx (s)k kQy (s)k ;

where, without lost of generality, k:k denotes the Frobenius norm. 9

Since (17) holds and Qx (s); Qy (s) are orthogonal matrices for all s 2 [t0 ; T ]; it is obtained that ¯ ¯ p ¯ f ¯ f ¯[Ax (s)]ii ¡ [Ay (s)]ii ¯ · 2k1 m kQx (s) ¡ Qy (s)k + l1 (k0 + M)mh:

(20)

Solving the equation (6) it can be shown that there is a positive constant P such that kX(t)¡1k · P for all t 2 [t0 ; T ] : Now, from Lemma 9 and Theorem 4 it follows that, for all h < h0 and h0 < 2P1C , 3 kX(s)¡1 k kX(s) ¡ Y(s)k 1 ¡ 2 kX(s)¡1 k kX(s) ¡ Y(s)k 3P Ch · 1 ¡ 2P Ch · C1 h;

kQx (s) ¡ Qy (s)k ·

(21)

3P C where C1 = 1¡2P . Ch From (21) and (20) we obtain that

¯ i ¯ p ¯¸x (T ) ¡ ¸iy (T )¯ · (2 mk1C1 + l1(k0 + M )m)h;

(22)

which is just the …rst assertion of the theorem. Given ½i ; »i > 0; from Lemma 3 there exist constants Ki ; Li > 1 such that ¯ i ¯ ¯¸x ¡ ¸ix (T )¯ · ai ; ¯ ¯ ¯o n¯ ¯ ¯ ¯ ¯ 1 1 » + where ai = max ¯½i + T ¡t ; log(L ) log(K ) i ¯ ¯ i i ¯ . Then, the second asT ¡t0 0 sertion of the theorem follows directly from ¯ ¯ ¯ ¯ ¯ ¯ i ¯ ¯ ¯¸x ¡ ¸iy (T )¯ · ¯¸ix ¡ ¸ix (T )¯ + ¯¸ix (T ) ¡ ¸iy (T )¯ · ai + ¯¸ix (T ) ¡ ¸iy (T )¯

and (22). ¯ ¯ Note that in the last inequality the errors ¯¸ix ¡ ¸iy (T )¯ are split into two components: ai and bh. The former represents the error introduced by the truncation of the in…nite time interval in the de…nition of Lyapunov exponents, while the later represents the error due to the approximation of X by Y.

5

Numerical experiments

In this section, the performance of the numerical method introduced above is illustrated by means of an example. The numerical implementation of the discrete QR decomposition based method follows that described in [19]. 10

5.1

Algorithm description

Let E0 be the error tolerance in the computation of the approximate Lyapunov exponents. Let T0 ; Tmax and ¢T · Tmax ¡ T0 be, respectively, the minimum, maximum and increment truncation times for which the truncated Lyapunov exponents will be computed. Steps of the algorithm: I. Set initial values for T0 , Tmax ; ¢T and E0 . II. Set n = j = 0, S = Li1 = 0, Q0 = I and Ea = Er = 1 + E0: III. Repeat the computation of: 1) The solution y of (5) at the point tj+1 2) The solution Z of :

Z = fx0 (tj ; y(tj ))Z; Z(tj ) = Qj ; t 2 [tj ; tj+1 ] at the point tj+1 3) The QR factorization Z(tj+1 ) = Qj+1 Rj+1 4) The expression for all i = 1; :::; m

S = S + log j[Rj+1 ]ii j

5) The expressions S t ¡t ¯j i 0 i¯ i Ea = ¯L1 ¡ L2 ¯ Ea Eri = i (for Li1 6= 0) jL1 j i L1 = Li2 n=n + 1 Li2 =

for all i, if tj ¸ T0 + n¢T

6) Update j = j + 1

until Eri < E0 or Eai < E0 for all i or tj > Tmax IV. Set T = tj and ¸iy (T ) = Li2 for all i: 11

(23)

In this algorithm, each value tj could be …xed in advance or adaptively chosen in the step III.1 by using some criterium to control the error in the numerical computation of y. In the case of autonomous ODEs, the solutions y and Z in the steps III.1 and III.2 of the algorithm are computed by using the results obtained in [22]. Namely, y(tj+1) = y(tj ) + g (y(tj ));

(24)

Z(tj+1) = E(y(tj )) Qj ;

(25)

and where the vector g(y(tn )) and the matrix E(y(tj )) are de…ned by the expression · ¸ E(y(tj )) g(y(tj )) = exp((tj+1 ¡ tj )C) (26) 0 1 with

·

¸ fx0 (y(tj )) f(y(tj )) C= 2 R(m+1)£(m+1): 0 0

In the general case of non-autonomous ODEs, the solutions y and Z in the steps III.1 and III.2 are also computed by using the results obtained in [22]. That is, y(tj+1 ) = y(tj ) + q(y(tj )) and Z(tj+1) = E(y(tj )) Qj ; where the vectors g(y(tn )); q(y(tn )) and the matrix E(y(tj )) are de…ned by the expression 2 3 E(y(tj )) g(y(tj )) q(y(tn )) 4 0 1 p(y(tn )) 5 = exp((tj+1 ¡ tj )C) 0 0 1

with

5.2

2

3 fx0 (y(tj )) ft0 (y(tj )) f(y(tj )) 0 0 1 5 2 R(m+2)£(m+2): C=4 0 0 0

Numerical Example

Consider the non-autonomous ODE de…ned by the equation [19]: :

x= A(t)x;

(27)

where t ¸ t0 = 0; A(t) = Q(t)B(t)Q(t)> + Q0t (t)Q(t)> with B(t) = diag(1; cos(t); ¡ 2p1t+1 ; ¡10), and Q(t) = diag(1; Q¯ ; 1)¢diag(Q® (t); Q® (t)): Here ® = 1; µ ¶ p cos (°t) sin (°t) ¯ = 2 and Q° (t) = : ¡ sin (°t) cos (°t) 12

¸4x

In this system, the Lyapunov exponents are [19]: ¸1x = 1; ¸2x = 0; ¸3x = 0; = ¡10 for all initial points x:

5.3

Results

Table I shows the values of ¸iy (T ) computed by the algorithm for various values of T and h taking discrete times tj = t0 +jh: In this case, it was set E0 = ¢T = 0 and Tmax = T0 = T . Table II shows the errors E between the approximate and exact Lyapunov exponents. Note that the error decreases as T increases and h decreases. Fig 1 and 2 plot these errors vs h. For each ¸x and two values of T , Table III shows the estimates b a, bb and b c of the parameters a; b and c resulting from …tting the power function curve E¸x (h) = a + bhc

(28)

to the data presented in Table II. Figures 1 and 2 display the …tted curves obtained for each Lyapunov Exponent. Note that in both cases the estimates b c are in agreement with the theoretical result c = 1 stated in theorem 8. Besides, the estimates b a decreases as T increases, as is expected according to this theorem. For various values of E0 and the step size h, Table IV shows the values of i ¸y (T ) computed by the algorithm taking uniformly spaced times tj = t0 + jh: Here, we set T0 = 100, ¢T = 100 and Tmax = 30 000: In addition, this table shows the values of T reached by the algorithm in the computation of the Lyapunov exponents with the precision E0:

13

References [1] A. Lyapunov, Problém Géneral de la Stabilité du Mouvement, Ann. Math. Stud, 17 Princeton University Press, Princeton, NJ, 1949. [2] V.I.Osedelec, A multiplicative ergodic theorem. Lyapunov characteristic exponents for dynamical systems,Trudy Moskov. Mat. Obshch,19 (1968), 197231. [3] D. Ruelle, Chaotic evolution and strange attractors. Cambridge University Press, 1995. [4] K. Geist, U. Parlitz and W. Lauterborn. Comparison of di¤erent methods for computing Lyapunov exponents, Progr. Theoret. Phys., 83 (1990), 875-893. [5] E.A.Coddington, N.Levinson, Theory of Ordinary Di¤erential Equations. McGraw-Hill, New York,1955. [6] T. Barker, R. Bowles and W. Williams, Development and application of a local linearization algorithm for integration of quaternion rate equations in real-time ‡ight simulation problems, NASA TN-D-7347, 1973. [7] M. Fjeld and S. Aam, An implementation of estimation techniques to a hydrological model for prediction of runo¤ to a hydroelectric power station, IEEE Trans. Autom. Control, AC-25 (1980), 151-163. [8] T. Ozaki, A bridge between nonlinear time series models and nonlinear stochastic dynamical systems: a local linearization approach, Statistica Sinica, 2 (1992), 113-135. [9] J. Smith, Mathematical modeling and digital simulations for engineers and scientists, John Wiley, 1977. [10] J.I. Ramos and C.M. Garcia-Lopez, Piecewise-Linearized methods for initialvalue problems, Appl. Math. Comput., 82 (1997) 273-302. [11] M. Hochbruck and C. Lubich, On Krylov subspace approximations to the matrix exponential operator, SIAM Numer. Anal., 34 (1997), 1911-1925. [12] M. Hochbruck, C. Lubich and H. Selhofer, Exponential integrators for large systems of di¤erential equations, SIAM J. Sci. Comp., 19 (1998), 1552-1574. [13] J.M. Bower and D. Beeman, The book of GENESIS: exploring realistic neural model with the general neural simulation system. Springer-Verlarg, 1995. [14] R.J. MacGregor, Neural and brain modeling, Academic Press, 1987. [15] J.C.Jimenez, C.Mora and R.Biscay, Properties of a Local Linearization scheme for the numerical integration of ODEs, Technical Report No. I-98-06 (MB/CIMAT). Submitted. [16] R. Biscay, J.C. Jimenez, J. Riera and P. Valdes, Local linearization method for the numerical solution of stochastic di¤erential equations, Annals Inst. Statist.. Math. 48 (1996), 631-644. [17] R.B. Sidje, EXPOKIT: software package for computing matrix exponentials, AMC Trans. Math. Software, 24 (1998), 130:156. 14

[18] J.C. Jimenez and R. Biscay, Converge of the Local Linearization method for stochatic di¤erential equations revisited. Technical Report ICIMAF 2000-97, Cuba, 2000. Submitted. [19] L. Dieci, R. D. Russell and E. S. Van Vleck, On the computation of Lyapunov exponents for continuous dynamical systems, SIAM J. Numer. Anal., 34 (1997) 402-423. [20] L. Dieci and E. S. Van Vleck, Computing of a few Lyapunov exponents for continuous and discrete dynamical systems, Appl. Num. Math., 17 (1995), 275-291. [21] G. W. Stewart, Perturbation bounds for the QR factorization of a matrix, SIAM J. Numerical Analysis,14 (1977), 509-518. [22] C.F. Van Loan, Computing integrals involving the matrix exponential, IEEE Trans. Autom. Control, AC-23 (1978), 395-404.

15

FIGURE CAPTION Figure 1. Plots of the errors E¸ix between Lyapunov exponents ¸ix associated to the ODE (27) and the truncate Lyapunov Exponents ¸iy (T = 500) of the LL approximation, for di¤erent values of h. Figure 2. Plots of the errors E¸ix between Lyapunov exponents ¸ix associated to the ODE (27) and the truncate Lyapunov Exponents ¸iy (T = 1000) of the LL approximation, for di¤erent values of h.

16

TABLE CAPTION Table I: Values of the approximate Lyapunov Exponents ¸iy (T ) computed for di¤erent values of h and T . Table II: Errors E¸ix between Lyapunov exponents ¸ix associated to the ODE (27) and the approximate Lyapunov Exponents ¸iy (T ) for di¤erent values of h and T . Table III: Estimated values b a, bb and b c resulting from …ts the model E¸ix (h) = c a + bh to the data of Table II. Table IV: Values of the approximate Lyapunov exponents ¸iy (T ) computed for various values of E0 and h and the values of T reached by the algorithm.

17

T = 500 h 2¡4 2¡5 2¡6 2¡7 2¡8

¸1y (T ) 1.00043466 1.00041209 1.00040070 1.00039499 1.00039212

¸2y (T ) -0.00283596 -0.00282044 -0.00281191 -0.00280745 -0.00280518

T = 1000 h 2¡4 2¡5 2¡6 2¡7 2¡8

¸1y (T ) 0.99977229 0.99977671 0.99977898 0.99978013 0.99978071

¸2y (T ) 0.00068418 0.00066417 0.00065455 0.00064985 0.00064752 Table I

18

¸3y (T ) -0.02221151 -0.02214130 -0.02210584 -0.02208802 -0.02207909 ¸3y (T ) -0.01621231 -0.01618184 -0.01616615 -0.01615818 -0.01615417

¸4y (T ) -10.00972202 -10.009217089 -10.008962469 -10.008834622 -10.008770563 ¸4y (T ) -10.00555336 -10.00544538 -10.00539009 -10.00536211 -10.00534805

T = 500 h 2¡4 2¡5 2¡6 2¡7 2¡8

E¸1x (h) 0.00043466 0.00041209 0.00040070 0.00039499 0.00039212

E¸2x (h) 0.00283596 0.00282044 0.00281191 0.00280745 0.00280518

E¸3x (h) 0.02221151 0.02214130 0.02210584 0.02208802 0.02207909

E¸4x (h) 0.00972202 0.009217089 0.008962469 0.008834622 0.008770563

T = 1000 h 2¡4 2¡5 2¡6 2¡7 2¡8

E¸1x (h) 0.00022771 0.00022329 0.00022102 0.00021987 0.00021929

E¸2x (h) 0.00068418 0.00066417 0.00065455 0.00064985 0.00064752

E¸3x (h) 0.01621231 0.01618184 0.01616615 0.01615818 0.01615417

E¸4x (h) 0.00555336 0.00544538 0.00539009 0.00536211 0.00534805

Table II

19

T = 500 b a 0.000389 ¸1 0.002803 ¸2 0.022073 ¸3 0.008716 ¸4 T = 1000

bb 0.000771 0.000543 0.002648 0.017368

b c 1.016272 1.019264 1.065600 1.027164

b a 0.000218 0.000645 0.016151 0.005335

bb 0.000153 0.000655 0.001178 0.003777

b c 1.012373 1.018690 1.062024 1.022229

¸1 ¸2 ¸3 ¸4

Table III

20

h = 0:1 E0 10¡1 10¡2 10¡3 10¡4 h = 0:01 E0 10¡1 10¡2 10¡3 10¡4 h = 0:001 E0 10¡1 10¡2 10¡3 10¡4

T 200 500 1300 8600

¸1y (T ) 0.99471224 0.99976105 1.00013928 1.00002411

¸2y (T ) -0.00432135 -0.00437102 -0.00030678 0.00006761

T 200 500 1400 16200

¸1y (T ) 0.99411104 0.99920122 0.99916796 1.00001302

¸2y (T ) -0.00549059 -0.00431259 -0.00048698 -0.00005923

T 200 400 1400 14200

¸1y (T ) ¸2y (T ) -0.00555405 0.99408486 -0.00427181 0.99925455 -0.00049768 0.99914875 -0.00008251 1.00000514 Table IV

21

¸3y (T ) -0.03681002 -0.02334343 -0.01405960 -0.00541085 ¸3y (T ) -0.03709679 -0.02590003 -0.01369815 -0.00394468 ¸3y (T ) -0.037114603 -0.025851580 -0.013698698 -0.004234950

¸4y (T ) -9.98139310 -10.0118446 -10.0043157 -10.0006295 ¸4y (T ) -9.98234285 -10.0127662 -9.99735114 -10.0002439 ¸4y (T ) -9.98243938 -10.0128825 -9.99761727 -10.0003795

-4

4.6

-3

x 10

2.86

λ 1x

x 10

λ2

2.85

x

4.4 Eλ (h)

Eλ (h)

2.84 4.2

2.83 2.82

4 2.81 3.8

0

0.02

0.04 h

0.06

2.8

0.08

0

0.02

0.04 h

0.06

0.08

0.02

0.04 h

0.06

0.08

-3

0.0223

10.5

x 10

λ4

λ3

0.0223

x

x

10 Eλ (h)

Eλ (h)

0.0222 0.0222

9.5

0.0221 9 0.0221 0.022

0

0.02

0.04 h

0.06

0.08

8.5

0

-4

2.32

-4

x 10

7 λ1

2.3

x 10

2

λx

6.9

x

2.28

6.8

Eλ (h)

Eλ (h)

2.26

6.7

2.24

6.6

2.22

6.5

2.2 2.18

0

0.02

0.04

0.06

6.4

0.08

0

0.02

h

0.04

0.06

0.08

0.06

0.08

h -3

0.0162

5.65 3

λx

0.0162

x 10

λ4

5.6

x

5.55

Eλ (h)

Eλ (h)

0.0162

5.5

5.45

0.0162

5.4 0.0162 0.0161

5.35 0

0.02

0.04

h

0.06

0.08

5.3

0

0.02

0.04

h