- Email: [email protected]
An asymptotically exact stopping rule for the numerical computation of the Lyapunov spectrum
Chaos, Solitons
& Fractals Vol. 7, No. 8, pp. 1213-1225, 1996 Coowieht CT, 19% Elsevier Science Ltd Printed I; &a;-Britain. All rights reserved o!+5om79/96 $15.00 + 0.00
09fio-0779(95)00107-7
An Asymptotically Exact Stopping Rule for the Numerical Computation of the Lyapunov Spectrum JELEL EZZINE Department of Systems Engineering, King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia (Accepted 9 November 1995)
is in general not possible to analytically compute the Lyapunov spectrum of a given dynamical system. This has been achieved for a few special cases only. Therefore, numerical algorithms have been devised for this task. However, one rnajor drawback of these numerical algorithms is their lack of stopping rules. In this paper, an asymptotically exact stopping rule is proposed to alleviate this shortcoming while computing the Lyapunov spectrum of linear discretetime random dynamical systems (i.e., linear systems with random parameters). The proposed stopping rule provides an estimate of the least number of iterations, for which the probability of incurring a prescribed error, in the numerical computation of the Lyapunov spectrum, is minimized. It exploits simple upper bounds on the Lyapunov exponents, along with some results from finite state Markov chains. The accuracy of the stopping rule, and the computational load, is proportional to the tightness of the bound. In fact, a series of increasingly tighter bounds are proposed, yielding an asymptotically exact stopping rule for the tightest one. It is demonstrated via an example, that the proposed stopping rule is applicable to nonlinear dynamics as well. Copyright @ 1996 Elsevier Science Ltd Abstract-It
1. INTRODUCTION
One century ago, Lyapunov [l] attempted to define an eigenvalue-like concept for linear time-varying dynamical systems, but his findings applied only to a very special class of systems he called regular. Regularity is a property hard to check. However, recently, Oseledest [l] showed that it is an almost sure (a.s.) property for random dynamical systems, and developed a useful eigenstructure-like theory for this class of systems. Though Lyapunov exponents are a very powerful tool for the analysis of dynamical systems, they suffer from an important drawback as regards their analytical as well as numerical computations. The Lyapunov exponents were analytically computed for a few simple dynamical systems only. The only available means to compute these numbers is via numerical algorithms. These numerical algorithms use the actual time evolution of the states of the system, thus requiring extended simulation time of the dynamics in question. Theoretically, an infinite time horizon is required to compute these numbers, thus, practically, there is an acute need for adequate stopping rules to alleviate this major limitation. Recently, this absence of a stopping rule was clearly acknowledged in Lichtenberg, and Lieberman’s book [2], and Talay’s recent work, such as [3], where it was stated in the former: “This points out a major difficulty in numerically evaluating the Lyapunov exponents: there is no a priori condition for determining the number of iterations it that must be used.” In addition, Eckmann and Ruelle [4], questioned the validity of a number of recently published estimates of dimensions of attractors, and Lyapunov exponents, based on rather short time series. They concluded that long time series, of high quality, are 1213
J. EZZINE
1214
necessary for such estimates, and that the number of points needed to estimate the Lyapunov exponents is about the square of that needed to estimate dimension. It is the objective of this paper to propose a stopping rule for the computation of the largest Lyapunov exponent for linear random dynamical systems, i.e., systems with random parameters. The computation of the remaining spectrum is similar to that of the largest Lyapunov exponent, and thus, will be omitted. As illustrated in the paper, the same stopping rule can be used to extract useful dynamical information from other classes of dynamical systems. The next section introduces the concept of discrete-time random dynamical systems, along with the definition of the largest Lyapunov exponent. In addition it provides the needed definitions, along with some preliminary results required to derive the main outcome of the paper. Section 3 states the main theorem, along with a short intuitive discussion of its key parameters. It also discusses the important scalar case of the proposed theorem with some numerical illustrations. Section 4 applies the result of the paper to a different problem, namely to iterated maps on the interval. Section 5 summarizes and concludes the paper. 2. DISCRETE-TIME RANDOM DYNAMICAL SYSTEMS
In this paper, attention is restricted to linear discrete-time systems with randomly varying parameters. This class of systems will be called linear discrete-time random (DTR) dynamical systems in contrast to linear discrete-time stochastic systems where randomness enters the system as an input and/or measurement noise. More precisely, the paper deals with the study of the following problem: Let {M,, n E N} be a sequence of random, IZ X n matrices. To each x0 E R” associate the process {X,3 IZ E N} with values in R”, which is the solution to X n+l
=
MJ,,
n E
N,
x0
=
x0.
(1)
We have Xn+l = M, * * . M,xo. The asymptotic behaviour of this process is addressed by the following theorem. Fzmtenberg-Kesten Theorem (FKT) [5]. Let {ML, i E N} be a stationary, metrically transitive stochastic process with values in the set of IZ x n matrices such that VI-a+ llMoll> < m. Then, with probability one
(2) := A17
where S, := M, *.. Ml and a1 E R U {-m}, is the largest Lyapunov exponent of the process. As is stated in the FKT, the largest Lyapunov exponent, A,, is an infimum of a particular set. Consequently, if any of the elements of this set is negative, one concludes, using the property of the i&mum, that the system is a.s. exponentially stable. However, by not knowing exactly A-,, it is not possible to tell how stable the system is! That is, by exploiting this property to alleviate the computational burden, some qualitative insight about the dynamical behaviour of the system is lost. This qualitative insight can be crucial for application purposes.
Stopping rule for the Lyapunov spectrum
1215
As an attempt to better exploit this theorem, we provide a stopping rule for the numerical computation of the Lyapunov spectrum. The stopping rule is based on simple upper bounds on hr. 2.1. More definitions and preliminary
results
In this subsection, definitions as well as basic results are introduced. In the sequel, it is assumed that the matrix of the homogeneous part of the random dynamical system is switching, or jumping, among a finite number k of matrices, following a finite state Markov chain (FSMC). The matrix Mj will be called the ith form matrix, and the finite set to which it belongs is called the matrix-form set, and is denoted by d [6]. This class of DTR dynamical systems will be called k-form DTR dynamical systems. The finite set, of all possible matrix products of length 1, whose factors are elements of J& is denoted by M 1’1. This set contains k’ elements. These elements will be referred to as Z-blocks, and denoted by M\‘], with i = 1, . . ., k’. At this point, simple, but key preliminary results are stated. A series of increasingly tighter bounds on the largest Lyapunov exponent is provided. These bounds provide the motivation for the main result of the paper which will be introduced in the next section. Proposition. The largest Lyapunov exponent, Al, of a homogeneous k-form DTR dynamical system, with a stationary irreducible FSMC satisfies the following inequality
(3) Moreover, $?‘I s #I 7
q, 1 E N\(O).
Proof. The first part is trivial. In addition, it follows from FKT that lim,,, last part, is a consequence of the fact that IA-[‘1is a subadditive sequence [l].
(4)
Al’1= hr. The n
The last part of the proposition is quite important, especially from a practical view point. Though the first part of the proposition asserts that the proposed sequences of upper bounds eventually converge to the largest Lyapunov exponent, they shed no light on the nature of this convergence! For instance, it would be of practical interest to identify, at least one monotonically convergent subsequence. This subsequence would allow one to systematically pick increasingly tighter bounds. For example, the following is such a subsequence: ~[@“I G i1[@‘1,i E N. The last result asserts that there is a countably infinite number of monotonically converging subsequences, which are easily identifiable, as shown in the last paragraph. Moreover, all the examples used in this work exhibited monotonically convergent A[‘] sequences, thus suggesting to conjecture that this sequence is indeed monotonically convergent, in the context of this paper.
3. THE KATZ-THOMASIAN
(KT) STOPPING RULE
Since the proposed stopping rule is motivated by the work of Katz and Thomasian [7], we shall call it the Katz-Thomasian (KT) stopping rule. It provides a probabilistic stopping criterion for the numerical computation of A[‘], using time-averaging. In other words, it supplies an estimate of the least number of matrix multiplications, needed in the
J. EZZINE
1216
time-average computation of il 1’1, for which the probability minimized. An earlier version of this result appeared in [8].
of having a large error is
Given a homogeneous k-form DTR dynamical system with a stationary and irreducible FSMC, a positive integer m and E > 0. Then Theorem.
(5) CC=
8k’ ~~‘(1 - e-L@)
(6)
(7) CT= fj~wdlll/
- mi;LogII#/(},
(8)
where p is the smallest entry in the FSMC probability transition matrix, and pi” := with pi, the entry of the FSMC steady state probability vector, corresponding PilIl~~:Pi.i.+,, to the first’matrix of the ith l-block. Proof. The proof readily follows from applying the main theorem in [7], to the k’ states Markov chain governing the random sequence of the Z-block matrices (see Doob [9, p. 1851. n
The KT stopping rule provides a probabilistic stopping criterion for the numerical computation of an upper bound on A-i while still obtaining a good approximation. This might seem useless since we have a simple expression for Al’]. However, the All1is an upper bound on A,, therefore, this stopping rule can be used as an approximate one in computing ill. It is important to note that Al’1converges to A1for sufficiently large 1, which makes the stopping rule asymptotically exact. To fully exploit the above result, one needs to understand the role of its decision parameters! The main decision parameter is the error tolerance denoted by E. The error tolerance, E, reflects the user’s willingness to allow for an acceptable mismatch between the exact, but unknown, result and the approximate one. The next decision parameter is the exit probability denoted by PL!Jt. After specifying the error tolerance, the user has to decide how often should this error be allowed, and this is accomplished by imposing a value on @‘I. em* It is obvious that the user would like to choose quite a small error tolerance, and exit probability. Unfortunately, this choice will always result in an unacceptably large number of iterations. To avoid such outcomes, the user has to strike a balance between E and P$t, and this can be accomplished by choosing these decision parameters in an inversely proportional manner. That is, allowing one of these parameters to be small, should automatically let the other take reasonably large magnitudes. To help remember this decision strategy, it is worth keeping in mind that, as a rule of thumb, the number of iterations is similar to l(l/E2) log { PLJii,}1. It is also possible to use the theorem by committing oneself to a chosen number of iterations, and then use the theorem to either estimate PL:i,, by fixing E, or conversely. Which route to follow depends on the problem at hand, as well as the user’s preference.
Stopping rule for the Lyapunov spectrum
3.1.
1217
The scalar case
As seen above, the proposed stopping rule is exact only asymptotically in 1. This limitation is eliminated when dealing with scalar dynamics, for in this case A1= A[‘]. Though we are not directly interested in this case, there is more than one application where scalar dynamics are of interest, e.g., one-dimensional maps [2], such as the logistic map. Evidently, the proposed result is applicable after minor manipulations, as will be seen in the sequel. Applying the stopping rule to scalar systems is a good opportunity to check its effectiveness as well as conservativeness. Moreover, the interrelations among the different parameters involved in the stopping rule can be investigated to better understand the strength and weaknesses of such a rule. Probably, the most important parameter in the stopping rule is 6. In fact, this parameter measures how close the process under investigation is to a one-valued, or constant process, and subsequently alters the estimate of the number of iterations needed to achieve a required accuracy accordingly. For instance, for 6 = 0, i.e., constant process PF.it = 0, implying instantaneous equality between the time and the ensemble averages, which is expected. For large 6 values, the opposite scenario is observed. A typical interrelationship between the number of iterations needed to achieve a given Posit and 6 is depicted in Fig. l.*
Fig. 1. Probability of exit vs number of iterations and delta +In Figs 1-3, E = 0.1, 6 = 0.5, and p = 0.5.
1218
J. EZZINE
In a like manner, the number of forms, i.e., k’, among which the process is switching, is also very important in the estimation of the number of iterations required for a target accuracy. Again, as expected, the more forms, the more iterations are required to achieve a given accuracy. That is, sufficient time is needed by the process to exhibit all these values according to their respective probabilities. Figure 2 shows such dependence between P$, k’, and the required number of iterations. The figure suggests that the required number of iterations grows quite fast as a function of the number of forms, thus substantiating the idea that a large number of iterations are always needed for such computations. As for the remaining parameter, p, used in the stopping rule, it complements the number of forms issue. Recall that this entry is the smallest entry in the probability transition matrix governing the switching between the forms, therefore, it reflects how infrequent the switching is between a given form or value, and the remaining ones. Obviously, the smaller this value is, the larger is the number of iterations required to achieve a desired accuracy, for sufficient time is to be allocated for all states to be sufficiently visited. This property is illustrated in Fig. 3. Finally, to take advantage of the applicability of the law of large numbers and illustrate the evolution of the time average, and compare it to its exact ensemble average, a scalar case where k’ = 2, 6 = 2, and p = 0.5 is used. To achieve a Z’exit= 0.1 and an E = 0.1, the KT stopping rule indicated that at least lo8 iterations are needed, which is in concordance with Figs l-3. However, as shown in Fig. 4, this estimate is quite conservative, since it took about 1500 iterations to achieve an E < O.l!
Fig. 2. Probability of exit vs number of iterations and number of forms.
Stopping rule for the Lyapunov spectrum
8.
lo-
1. log Fig. 3. Probability of exit vs number of iterations and p.
0.6
-
0.4
-
0.2
-
O-
Fig. 4. Lyapunov exponent for scalar dynamics.
1219
1220
J. EZZINE
Moreover, this run reveals the same typical qualitative dynamical behaviour of the successive iterates of the time average. The transient portion of the graph, i.e., the first few hundred or so, are very oscillatory. Nevertheless, the subsequent iterates approach the neighbourhood of the exact result quite fast, but keep fluctuating without ever seeming to converge. This last attribute attests to the complex dynamical nature of this numerical process, and could partly explain the conservativeness of the KT stopping rule. It is important to note that the above discussion is possible only in the scalar case, when the time and ensemble averages are identical as a result of the law of large numbers. Nevertheless, the argument remains valid for the upper bounds, where the tighter ones sketch an accurate picture of the Lyapunov exponent computational algorithm with respect to the mechanics indicated above. 3.2.
Numerical illustrations
After having developed a feel for the KT stopping rule, and its parameters, the criterion will be used in this section to accurately estimate the largest Lyapunov exponent of a two-form, two-dimensional DTR dynamical system, with 6 = 2, and p = 0.5. Then, the result will be compared with simulation findings. As depicted by the lower curve of Fig. 5, it took about 4 X ld iterations to converge to the largest Lyapunov exponent. Despite this acceptable convergence, as shown in the same figure, the iterates keep fluctuating very tightly about the correct value! These seemingly never fading fluctuations could be the reason behind the rather conservative number of iterations, i.e., lo’, estimated from the FK stopping rule. The above variations are typical of such computations, and more interestingly, evolve exactly the same way as the successive estimates of the l-block upper bound. Of course, both numerical iterations are for the same realization of the switching process. Figure 6
0.06
-0.06 Fig. 5. Largest Lyapunov exponent and the l-block upper bound.
Stopping rule for the Lyapunov spectrum
Fig. 6. Gap between the largest Lyapunov exponent and the l-block upper bound.
.shows clearly this synchronization, yielding an almost, and certainly asymptotically, constant difference between both numerical processes. The above observation is quite important from a practical viewpoint, for it suggests using I-block approximations to estimate the needed number of iterations to yield a desired accuracy. In turn, this will reduce from the conservativeness of the stopping rule.
4. APPLICATION TO ITERATED MAPS
The class of systems for which the stopping rule was derived is a special class of dynamical systems. Nevertheless, this stopping rule can be applied to a multitude of other systems, for which Lyapunov exponents and related dynamical attributes are of interest. To illustrate this possibility, scalar iterated maps on the interval are used. 4.1.
The skewed-tent map (STM)
To illustrate the usefulness of the proposed stopping rule beyond DTR dynamical systems, a variation on the tent map (TM) is proposed, and the Lyapunov exponent of the dynamics induced by iterating the STM is computed using time averaging. To estimate the least number of iterations needed to achieve a prescribed accuracy the KT stopping rule could be used. As will be shown next, few preparatory steps are needed for the proposed dynamical system is not described in a suitable format for the application of the stopping rule. First, the invariant probability distribution of the induced dynamics is computed, then a probabilistic description of these dynamics is introduced.
J. EZZINE
1222
1
The skewed-tent map is defined on [0, 11, for r E (-i, k), as follows:
2
f(x) :=
ifxe
TX
[
2
p-x+----1 + 22
ifxe 1 + 2t
I
1 - 22 O,---
2
lB2r -9
2
1
1. I
For z = i, the map is reduced to the usual tent map. So, to compute the Lyapunov exponent of the dynamical system induced by iterating the STM, one proceeds along the same lines as for the TM. First, by using a classical theorem by Lasota and York [lo], one asserts that there does exist an invariant measure for the dynamical system under consideration, then by using the well-known Frobenius-Perron operator [lo], one concludes that the invariant measure is nothing but the Lebesgue measure. As a matter of fact, this is the same invariant measure for all rs. Now, the unit interval is divided into two subintervals, I1 = [0, (1 - 2r/2)], and I, = [(l - 2x/2), 11, which will be referred to, sometimes, as ‘macrostates’. In addition, a stochastic matrix P is introduced, where pii is the fraction of interval Ii which is mapped into interval Ii by f(x), i.e., pij :=
m(zi
n
f-‘Czj))
m(zi)
’
where m(.) denotes the Lebesgue measure. The matrix P turned out to be
P= i
1 - 2t ___
1 + 2t ___
2
2
-1 - 2t
-1 + 22
2
2
The above manipulations have produced a ‘macroscopic’ probabilistic description of the ‘microscopic’ deterministic dynamics generated by iterating the STM. As a matter of fact, the resulting dynamical system is known as an iterated function system (IFS) [ll]. The exact Lyapunov exponent of the chaotic dynamics induced by the iterates of the STh4, readily follows from the probabilistic ‘macro’ model,
A(z)= i+bg
1 +2 22I [21 - 2t I +L++Jg [~)
where, as expected, n(O) = Log2, the Lyapunov exponent of the TM. As is shown in Fig. 7, the Lyapunov exponent assumes its largest value for t = 0, making the dynamics induced by TM most chaotic among the STMs. In contrast, they all induce the Lebesgue measure as their invariant measure. Notice that the invariance of the invariant measure as a function of the control parameter r, as well as the continuous dependence of the Lyapunov exponent on the same parameter, makes the STM an appropriate device to generate uniformly distributed numbers with different Lyapunov exponents! The dependence of the Lyapunov exponent on t, reflects the effect of the closeness of seeds, or initial conditions, on the asymptotic closeness of random number realizations. At this point one has all the necessary elements to investigate the dependence of the key
1223
Stopping rule for the Lyapunov spectrum
0.7 0.6 0.5 0.4 0.3 0.2
-~
I.-----
-----
0.1 0
I
-0.4
-0.2
0
0.2
0.4
Fig. 7. Lyapunov exponent for STh4 vs tau.
parameters of the KT stopping rule as a function of r. For instance, the separation measure, 6(r), between the minimum and maximum slopes is given by, a(z)=
1ogE. I I It is important to note that liml,l,cl/zj s(r) = ~0, and lim,O s(r) = 0, providing. the greatest variation range for such a parameter, thus yielding maximum information regarding its effect on the stopping rule; see Fig. 8. The next important parameter is the smallest entry p(r) in the probability transition matrix P(r), p(z) = (1 f r/2) E [0, $1. As a result, for 1rj + i, this probability takes very small values, concentrating the presence of the state of the system in one ‘macro’ state, making the visits to the other ‘macro’ state very scarce. Thus, to accurately estimate Lyapunov exponents via time averages in such situations, a substantial amount of time is indispensable! The approach sketched above for the STM is applicable to a more general class of systems, for instance to a continuous map since it can be approximated to any desired level of accuracy using piecewise linear maps. Moreover, using techniques as provided in [12], makes nonlinear multidimensional systems, yet another class of systems for which the same approach is applicable. 4.2.
Adaptive KT stopping rule
To stress the importance of stopping rules for the numerical or experimental, i.e., based on time series, computation of Lyapunov exponents for nonlinear dynamics, two numerical results, apropos the well-known Lorenz system with identical parameter values, are compared. In Wolf et al. [13], it was reported that (A,, 4, &) = (2.16, 0.00, -32.4), where
J. EZZINE
1224
-0.4
-0.2
0
0.2
0.4
Fig. 8. Separation measure delta (tau) vs tau.
in Aronson [14], it was found that (Ai, 4, &) = (2.164, 0.000, -32X1)! The discrepancy between these two numerical results is impo~~t, especially for 4, and thus is of theoretical as well as practical concern. Besides the different algorithms used to compute the Lyapunov spectrum, which is a secondary issue in this context, the number of points or iterations used in both references were substantially different. Wolf et al. [13] used 8192 (- 8.2 1s) points, where Aronson [14] used at least lo6 iterations.’ It is worth noting that in 1131, a rule of thumb was used to estimate the necessary number of points needed to achieve a ‘reasonable’ accuracy. This is the only instance, the author knows of where an attempt was made to devise an empirical stopping rule for this purpose. This rule says that about 10d - 30d points, where d is the Lyapunov dimension, are needed. For the Lorenz dynamics d = 2.07, and the required number of points estimated by this stopping rule is 102 - 932 points, which is about ten times less than the actual number of points used.” One way to alleviate this problem is to use an on-line adaptive version of the proposed KT stopping rule while estimating Lyapunov spectrums either from dynamical equations or time series. The following paragraph gives a rough sketch of this adaptive scheme. First, a warm-up stage is needed during which crude estimates of S and p, are generated, say either using a few hundred iterations, or Wolf et af. rule of thumb. These estimates with k' = 2, are used to compute the required initial number of iterations. A second stage of computation is initiated with the latter estimate as the needed number of iterations. +The number of iterations used for the Lorenz system was not given in [14], but for the H&on map it was reported to be 106. However, since the former dynamics have larger dimension, it is safe to assume that at least l@ iterations were needed for the reported numerical results. *As a side remark, and in the light of the results of this paper, the Wolf et aE. rule of thumb can only be used as a crude lower bound on the needed number of points or iterations.
Stopping rule for the Lyapunov spectrum
1225
During this second stage 6 and p are refined, and a new stopping time is estimated. This procedure is repeated until the successive stopping times become comparable. 5. SUMMARY AND CONCLUSIONS
In this paper an approximate stopping rule for the computation of the Lyapunov spectrum for linear random dynamical systems was proposed. Though the KT stopping rule is approximate it is asymptotically exact. The stopping rule involves all basic parameters that influence the dynamics under consideration. These parameters were explored for the special case of scalar dynamics, and were found to corraborate with theoretical as well as experimental understanding of complex dynamics. Even though the KT stopping rule is proposed for linear dynamics it was shown via an example that it could be used for nonlinear dynamics as well, after some basic manipulations. Nonetheless, when these manipulations are not easy or possible, an adaptive scheme of the KT stopping rule could be devised to estimate stopping times on-line. This adaptive scheme needs much work until it becomes practical. Some experimental results were presented in this work, but much more is still needed to fully understand the KT stopping rule and hopefully find clues to reduce its conservativeness. Acknowledgement--The
author acknowledges King Fahd University of Petroleum and Minerals for its support.
REFERENCES 1. L. Arnold and V. Wihstutz (eds), Lyupunov Exponents, Lecture Notes in Mathematics, No. 1186. Springer, NY (1986). 2. A. J. Lichtenberg and M. A. Liebermann, Regular and Chaotic Dynamics, Second Edition, Applied Mathematics Sciences, Vol. 38. Springer, NY (1992). 3. D. Talay, Approximation of upper Lyapunov exponents of bilinear stochastic differential systems, SIAM, J. on Num. Anul. 28(4), 1141 (1991). 4. J. P. Eckmann and D. Ruelle, Fundamental limitations for estimating dimensions and Lyapunov exponents in dynamical systems, IHES/P/90/11 (1990). 5. H. Furstenberg and R. Kesten, Product of random matrices. Ann. Math. Stats. 31,457 (1960). 6. J. Ezzine, On the stabilization and control of hybrid systems, Ph.D. thesis, Georgia Institute of Technology, (June 1989). 7. M. Katz and A. J. Thomasian, An exponential bound for functions of a Markov chain, Ann. Math. Stats. 31, 470 (1960). 8. J. Ezzine, On a practical stopping rule for the numerical computation of the Lyapunov spectrum, Proc. IEEE ZZth ACC, Chicago, IL (June 1992). 9. J. L. Doob, Stochastic Processes. John Wiley and Sons, NY (1953). 10. A. La&a and J. A. York, On the existence of invariant measures for piecewise monotonic transformations, Trans. Am. Math. Sot. 186, 481 (1973). 11. M. Bamsley, Fractals Everywhere. Academic Press, NY (1988). 12. G. Nicolis, Finite coarse-graining and Chapman-Kohnogorov equation in conservative dynamical systems, Chaos, Solitons & Fractak l(l), 25 (1991). 13. A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Determining Lyapunov exponents from a time series, Physica D 16, 285 (1985). .. 14. J. W. Aronson, CHAOS, A SUN-based program for analyzing chaotic systems, Computers in Physics 408 Wm.
& Fractals Vol. 7, No. 8, pp. 1213-1225, 1996 Coowieht CT, 19% Elsevier Science Ltd Printed I; &a;-Britain. All rights reserved o!+5om79/96 $15.00 + 0.00
09fio-0779(95)00107-7
An Asymptotically Exact Stopping Rule for the Numerical Computation of the Lyapunov Spectrum JELEL EZZINE Department of Systems Engineering, King Fahd University of Petroleum and Minerals, Dhahran 31261, Saudi Arabia (Accepted 9 November 1995)
is in general not possible to analytically compute the Lyapunov spectrum of a given dynamical system. This has been achieved for a few special cases only. Therefore, numerical algorithms have been devised for this task. However, one rnajor drawback of these numerical algorithms is their lack of stopping rules. In this paper, an asymptotically exact stopping rule is proposed to alleviate this shortcoming while computing the Lyapunov spectrum of linear discretetime random dynamical systems (i.e., linear systems with random parameters). The proposed stopping rule provides an estimate of the least number of iterations, for which the probability of incurring a prescribed error, in the numerical computation of the Lyapunov spectrum, is minimized. It exploits simple upper bounds on the Lyapunov exponents, along with some results from finite state Markov chains. The accuracy of the stopping rule, and the computational load, is proportional to the tightness of the bound. In fact, a series of increasingly tighter bounds are proposed, yielding an asymptotically exact stopping rule for the tightest one. It is demonstrated via an example, that the proposed stopping rule is applicable to nonlinear dynamics as well. Copyright @ 1996 Elsevier Science Ltd Abstract-It
1. INTRODUCTION
One century ago, Lyapunov [l] attempted to define an eigenvalue-like concept for linear time-varying dynamical systems, but his findings applied only to a very special class of systems he called regular. Regularity is a property hard to check. However, recently, Oseledest [l] showed that it is an almost sure (a.s.) property for random dynamical systems, and developed a useful eigenstructure-like theory for this class of systems. Though Lyapunov exponents are a very powerful tool for the analysis of dynamical systems, they suffer from an important drawback as regards their analytical as well as numerical computations. The Lyapunov exponents were analytically computed for a few simple dynamical systems only. The only available means to compute these numbers is via numerical algorithms. These numerical algorithms use the actual time evolution of the states of the system, thus requiring extended simulation time of the dynamics in question. Theoretically, an infinite time horizon is required to compute these numbers, thus, practically, there is an acute need for adequate stopping rules to alleviate this major limitation. Recently, this absence of a stopping rule was clearly acknowledged in Lichtenberg, and Lieberman’s book [2], and Talay’s recent work, such as [3], where it was stated in the former: “This points out a major difficulty in numerically evaluating the Lyapunov exponents: there is no a priori condition for determining the number of iterations it that must be used.” In addition, Eckmann and Ruelle [4], questioned the validity of a number of recently published estimates of dimensions of attractors, and Lyapunov exponents, based on rather short time series. They concluded that long time series, of high quality, are 1213
J. EZZINE
1214
necessary for such estimates, and that the number of points needed to estimate the Lyapunov exponents is about the square of that needed to estimate dimension. It is the objective of this paper to propose a stopping rule for the computation of the largest Lyapunov exponent for linear random dynamical systems, i.e., systems with random parameters. The computation of the remaining spectrum is similar to that of the largest Lyapunov exponent, and thus, will be omitted. As illustrated in the paper, the same stopping rule can be used to extract useful dynamical information from other classes of dynamical systems. The next section introduces the concept of discrete-time random dynamical systems, along with the definition of the largest Lyapunov exponent. In addition it provides the needed definitions, along with some preliminary results required to derive the main outcome of the paper. Section 3 states the main theorem, along with a short intuitive discussion of its key parameters. It also discusses the important scalar case of the proposed theorem with some numerical illustrations. Section 4 applies the result of the paper to a different problem, namely to iterated maps on the interval. Section 5 summarizes and concludes the paper. 2. DISCRETE-TIME RANDOM DYNAMICAL SYSTEMS
In this paper, attention is restricted to linear discrete-time systems with randomly varying parameters. This class of systems will be called linear discrete-time random (DTR) dynamical systems in contrast to linear discrete-time stochastic systems where randomness enters the system as an input and/or measurement noise. More precisely, the paper deals with the study of the following problem: Let {M,, n E N} be a sequence of random, IZ X n matrices. To each x0 E R” associate the process {X,3 IZ E N} with values in R”, which is the solution to X n+l
=
MJ,,
n E
N,
x0
=
x0.
(1)
We have Xn+l = M, * * . M,xo. The asymptotic behaviour of this process is addressed by the following theorem. Fzmtenberg-Kesten Theorem (FKT) [5]. Let {ML, i E N} be a stationary, metrically transitive stochastic process with values in the set of IZ x n matrices such that VI-a+ llMoll> < m. Then, with probability one
(2) := A17
where S, := M, *.. Ml and a1 E R U {-m}, is the largest Lyapunov exponent of the process. As is stated in the FKT, the largest Lyapunov exponent, A,, is an infimum of a particular set. Consequently, if any of the elements of this set is negative, one concludes, using the property of the i&mum, that the system is a.s. exponentially stable. However, by not knowing exactly A-,, it is not possible to tell how stable the system is! That is, by exploiting this property to alleviate the computational burden, some qualitative insight about the dynamical behaviour of the system is lost. This qualitative insight can be crucial for application purposes.
Stopping rule for the Lyapunov spectrum
1215
As an attempt to better exploit this theorem, we provide a stopping rule for the numerical computation of the Lyapunov spectrum. The stopping rule is based on simple upper bounds on hr. 2.1. More definitions and preliminary
results
In this subsection, definitions as well as basic results are introduced. In the sequel, it is assumed that the matrix of the homogeneous part of the random dynamical system is switching, or jumping, among a finite number k of matrices, following a finite state Markov chain (FSMC). The matrix Mj will be called the ith form matrix, and the finite set to which it belongs is called the matrix-form set, and is denoted by d [6]. This class of DTR dynamical systems will be called k-form DTR dynamical systems. The finite set, of all possible matrix products of length 1, whose factors are elements of J& is denoted by M 1’1. This set contains k’ elements. These elements will be referred to as Z-blocks, and denoted by M\‘], with i = 1, . . ., k’. At this point, simple, but key preliminary results are stated. A series of increasingly tighter bounds on the largest Lyapunov exponent is provided. These bounds provide the motivation for the main result of the paper which will be introduced in the next section. Proposition. The largest Lyapunov exponent, Al, of a homogeneous k-form DTR dynamical system, with a stationary irreducible FSMC satisfies the following inequality
(3) Moreover, $?‘I s #I 7
q, 1 E N\(O).
Proof. The first part is trivial. In addition, it follows from FKT that lim,,, last part, is a consequence of the fact that IA-[‘1is a subadditive sequence [l].
(4)
Al’1= hr. The n
The last part of the proposition is quite important, especially from a practical view point. Though the first part of the proposition asserts that the proposed sequences of upper bounds eventually converge to the largest Lyapunov exponent, they shed no light on the nature of this convergence! For instance, it would be of practical interest to identify, at least one monotonically convergent subsequence. This subsequence would allow one to systematically pick increasingly tighter bounds. For example, the following is such a subsequence: ~[@“I G i1[@‘1,i E N. The last result asserts that there is a countably infinite number of monotonically converging subsequences, which are easily identifiable, as shown in the last paragraph. Moreover, all the examples used in this work exhibited monotonically convergent A[‘] sequences, thus suggesting to conjecture that this sequence is indeed monotonically convergent, in the context of this paper.
3. THE KATZ-THOMASIAN
(KT) STOPPING RULE
Since the proposed stopping rule is motivated by the work of Katz and Thomasian [7], we shall call it the Katz-Thomasian (KT) stopping rule. It provides a probabilistic stopping criterion for the numerical computation of A[‘], using time-averaging. In other words, it supplies an estimate of the least number of matrix multiplications, needed in the
J. EZZINE
1216
time-average computation of il 1’1, for which the probability minimized. An earlier version of this result appeared in [8].
of having a large error is
Given a homogeneous k-form DTR dynamical system with a stationary and irreducible FSMC, a positive integer m and E > 0. Then Theorem.
(5) CC=
8k’ ~~‘(1 - e-L@)
(6)
(7) CT= fj~wdlll/
- mi;LogII#/(},
(8)
where p is the smallest entry in the FSMC probability transition matrix, and pi” := with pi, the entry of the FSMC steady state probability vector, corresponding PilIl~~:Pi.i.+,, to the first’matrix of the ith l-block. Proof. The proof readily follows from applying the main theorem in [7], to the k’ states Markov chain governing the random sequence of the Z-block matrices (see Doob [9, p. 1851. n
The KT stopping rule provides a probabilistic stopping criterion for the numerical computation of an upper bound on A-i while still obtaining a good approximation. This might seem useless since we have a simple expression for Al’]. However, the All1is an upper bound on A,, therefore, this stopping rule can be used as an approximate one in computing ill. It is important to note that Al’1converges to A1for sufficiently large 1, which makes the stopping rule asymptotically exact. To fully exploit the above result, one needs to understand the role of its decision parameters! The main decision parameter is the error tolerance denoted by E. The error tolerance, E, reflects the user’s willingness to allow for an acceptable mismatch between the exact, but unknown, result and the approximate one. The next decision parameter is the exit probability denoted by PL!Jt. After specifying the error tolerance, the user has to decide how often should this error be allowed, and this is accomplished by imposing a value on @‘I. em* It is obvious that the user would like to choose quite a small error tolerance, and exit probability. Unfortunately, this choice will always result in an unacceptably large number of iterations. To avoid such outcomes, the user has to strike a balance between E and P$t, and this can be accomplished by choosing these decision parameters in an inversely proportional manner. That is, allowing one of these parameters to be small, should automatically let the other take reasonably large magnitudes. To help remember this decision strategy, it is worth keeping in mind that, as a rule of thumb, the number of iterations is similar to l(l/E2) log { PLJii,}1. It is also possible to use the theorem by committing oneself to a chosen number of iterations, and then use the theorem to either estimate PL:i,, by fixing E, or conversely. Which route to follow depends on the problem at hand, as well as the user’s preference.
Stopping rule for the Lyapunov spectrum
3.1.
1217
The scalar case
As seen above, the proposed stopping rule is exact only asymptotically in 1. This limitation is eliminated when dealing with scalar dynamics, for in this case A1= A[‘]. Though we are not directly interested in this case, there is more than one application where scalar dynamics are of interest, e.g., one-dimensional maps [2], such as the logistic map. Evidently, the proposed result is applicable after minor manipulations, as will be seen in the sequel. Applying the stopping rule to scalar systems is a good opportunity to check its effectiveness as well as conservativeness. Moreover, the interrelations among the different parameters involved in the stopping rule can be investigated to better understand the strength and weaknesses of such a rule. Probably, the most important parameter in the stopping rule is 6. In fact, this parameter measures how close the process under investigation is to a one-valued, or constant process, and subsequently alters the estimate of the number of iterations needed to achieve a required accuracy accordingly. For instance, for 6 = 0, i.e., constant process PF.it = 0, implying instantaneous equality between the time and the ensemble averages, which is expected. For large 6 values, the opposite scenario is observed. A typical interrelationship between the number of iterations needed to achieve a given Posit and 6 is depicted in Fig. l.*
Fig. 1. Probability of exit vs number of iterations and delta +In Figs 1-3, E = 0.1, 6 = 0.5, and p = 0.5.
1218
J. EZZINE
In a like manner, the number of forms, i.e., k’, among which the process is switching, is also very important in the estimation of the number of iterations required for a target accuracy. Again, as expected, the more forms, the more iterations are required to achieve a given accuracy. That is, sufficient time is needed by the process to exhibit all these values according to their respective probabilities. Figure 2 shows such dependence between P$, k’, and the required number of iterations. The figure suggests that the required number of iterations grows quite fast as a function of the number of forms, thus substantiating the idea that a large number of iterations are always needed for such computations. As for the remaining parameter, p, used in the stopping rule, it complements the number of forms issue. Recall that this entry is the smallest entry in the probability transition matrix governing the switching between the forms, therefore, it reflects how infrequent the switching is between a given form or value, and the remaining ones. Obviously, the smaller this value is, the larger is the number of iterations required to achieve a desired accuracy, for sufficient time is to be allocated for all states to be sufficiently visited. This property is illustrated in Fig. 3. Finally, to take advantage of the applicability of the law of large numbers and illustrate the evolution of the time average, and compare it to its exact ensemble average, a scalar case where k’ = 2, 6 = 2, and p = 0.5 is used. To achieve a Z’exit= 0.1 and an E = 0.1, the KT stopping rule indicated that at least lo8 iterations are needed, which is in concordance with Figs l-3. However, as shown in Fig. 4, this estimate is quite conservative, since it took about 1500 iterations to achieve an E < O.l!
Fig. 2. Probability of exit vs number of iterations and number of forms.
Stopping rule for the Lyapunov spectrum
8.
lo-
1. log Fig. 3. Probability of exit vs number of iterations and p.
0.6
-
0.4
-
0.2
-
O-
Fig. 4. Lyapunov exponent for scalar dynamics.
1219
1220
J. EZZINE
Moreover, this run reveals the same typical qualitative dynamical behaviour of the successive iterates of the time average. The transient portion of the graph, i.e., the first few hundred or so, are very oscillatory. Nevertheless, the subsequent iterates approach the neighbourhood of the exact result quite fast, but keep fluctuating without ever seeming to converge. This last attribute attests to the complex dynamical nature of this numerical process, and could partly explain the conservativeness of the KT stopping rule. It is important to note that the above discussion is possible only in the scalar case, when the time and ensemble averages are identical as a result of the law of large numbers. Nevertheless, the argument remains valid for the upper bounds, where the tighter ones sketch an accurate picture of the Lyapunov exponent computational algorithm with respect to the mechanics indicated above. 3.2.
Numerical illustrations
After having developed a feel for the KT stopping rule, and its parameters, the criterion will be used in this section to accurately estimate the largest Lyapunov exponent of a two-form, two-dimensional DTR dynamical system, with 6 = 2, and p = 0.5. Then, the result will be compared with simulation findings. As depicted by the lower curve of Fig. 5, it took about 4 X ld iterations to converge to the largest Lyapunov exponent. Despite this acceptable convergence, as shown in the same figure, the iterates keep fluctuating very tightly about the correct value! These seemingly never fading fluctuations could be the reason behind the rather conservative number of iterations, i.e., lo’, estimated from the FK stopping rule. The above variations are typical of such computations, and more interestingly, evolve exactly the same way as the successive estimates of the l-block upper bound. Of course, both numerical iterations are for the same realization of the switching process. Figure 6
0.06
-0.06 Fig. 5. Largest Lyapunov exponent and the l-block upper bound.
Stopping rule for the Lyapunov spectrum
Fig. 6. Gap between the largest Lyapunov exponent and the l-block upper bound.
.shows clearly this synchronization, yielding an almost, and certainly asymptotically, constant difference between both numerical processes. The above observation is quite important from a practical viewpoint, for it suggests using I-block approximations to estimate the needed number of iterations to yield a desired accuracy. In turn, this will reduce from the conservativeness of the stopping rule.
4. APPLICATION TO ITERATED MAPS
The class of systems for which the stopping rule was derived is a special class of dynamical systems. Nevertheless, this stopping rule can be applied to a multitude of other systems, for which Lyapunov exponents and related dynamical attributes are of interest. To illustrate this possibility, scalar iterated maps on the interval are used. 4.1.
The skewed-tent map (STM)
To illustrate the usefulness of the proposed stopping rule beyond DTR dynamical systems, a variation on the tent map (TM) is proposed, and the Lyapunov exponent of the dynamics induced by iterating the STM is computed using time averaging. To estimate the least number of iterations needed to achieve a prescribed accuracy the KT stopping rule could be used. As will be shown next, few preparatory steps are needed for the proposed dynamical system is not described in a suitable format for the application of the stopping rule. First, the invariant probability distribution of the induced dynamics is computed, then a probabilistic description of these dynamics is introduced.
J. EZZINE
1222
1
The skewed-tent map is defined on [0, 11, for r E (-i, k), as follows:
2
f(x) :=
ifxe
TX
[
2
p-x+----1 + 22
ifxe 1 + 2t
I
1 - 22 O,---
2
lB2r -9
2
1
1. I
For z = i, the map is reduced to the usual tent map. So, to compute the Lyapunov exponent of the dynamical system induced by iterating the STM, one proceeds along the same lines as for the TM. First, by using a classical theorem by Lasota and York [lo], one asserts that there does exist an invariant measure for the dynamical system under consideration, then by using the well-known Frobenius-Perron operator [lo], one concludes that the invariant measure is nothing but the Lebesgue measure. As a matter of fact, this is the same invariant measure for all rs. Now, the unit interval is divided into two subintervals, I1 = [0, (1 - 2r/2)], and I, = [(l - 2x/2), 11, which will be referred to, sometimes, as ‘macrostates’. In addition, a stochastic matrix P is introduced, where pii is the fraction of interval Ii which is mapped into interval Ii by f(x), i.e., pij :=
m(zi
n
f-‘Czj))
m(zi)
’
where m(.) denotes the Lebesgue measure. The matrix P turned out to be
P= i
1 - 2t ___
1 + 2t ___
2
2
-1 - 2t
-1 + 22
2
2
The above manipulations have produced a ‘macroscopic’ probabilistic description of the ‘microscopic’ deterministic dynamics generated by iterating the STM. As a matter of fact, the resulting dynamical system is known as an iterated function system (IFS) [ll]. The exact Lyapunov exponent of the chaotic dynamics induced by the iterates of the STh4, readily follows from the probabilistic ‘macro’ model,
A(z)= i+bg
1 +2 22I [21 - 2t I +L++Jg [~)
where, as expected, n(O) = Log2, the Lyapunov exponent of the TM. As is shown in Fig. 7, the Lyapunov exponent assumes its largest value for t = 0, making the dynamics induced by TM most chaotic among the STMs. In contrast, they all induce the Lebesgue measure as their invariant measure. Notice that the invariance of the invariant measure as a function of the control parameter r, as well as the continuous dependence of the Lyapunov exponent on the same parameter, makes the STM an appropriate device to generate uniformly distributed numbers with different Lyapunov exponents! The dependence of the Lyapunov exponent on t, reflects the effect of the closeness of seeds, or initial conditions, on the asymptotic closeness of random number realizations. At this point one has all the necessary elements to investigate the dependence of the key
1223
Stopping rule for the Lyapunov spectrum
0.7 0.6 0.5 0.4 0.3 0.2
-~
I.-----
-----
0.1 0
I
-0.4
-0.2
0
0.2
0.4
Fig. 7. Lyapunov exponent for STh4 vs tau.
parameters of the KT stopping rule as a function of r. For instance, the separation measure, 6(r), between the minimum and maximum slopes is given by, a(z)=
1ogE. I I It is important to note that liml,l,cl/zj s(r) = ~0, and lim,O s(r) = 0, providing. the greatest variation range for such a parameter, thus yielding maximum information regarding its effect on the stopping rule; see Fig. 8. The next important parameter is the smallest entry p(r) in the probability transition matrix P(r), p(z) = (1 f r/2) E [0, $1. As a result, for 1rj + i, this probability takes very small values, concentrating the presence of the state of the system in one ‘macro’ state, making the visits to the other ‘macro’ state very scarce. Thus, to accurately estimate Lyapunov exponents via time averages in such situations, a substantial amount of time is indispensable! The approach sketched above for the STM is applicable to a more general class of systems, for instance to a continuous map since it can be approximated to any desired level of accuracy using piecewise linear maps. Moreover, using techniques as provided in [12], makes nonlinear multidimensional systems, yet another class of systems for which the same approach is applicable. 4.2.
Adaptive KT stopping rule
To stress the importance of stopping rules for the numerical or experimental, i.e., based on time series, computation of Lyapunov exponents for nonlinear dynamics, two numerical results, apropos the well-known Lorenz system with identical parameter values, are compared. In Wolf et al. [13], it was reported that (A,, 4, &) = (2.16, 0.00, -32.4), where
J. EZZINE
1224
-0.4
-0.2
0
0.2
0.4
Fig. 8. Separation measure delta (tau) vs tau.
in Aronson [14], it was found that (Ai, 4, &) = (2.164, 0.000, -32X1)! The discrepancy between these two numerical results is impo~~t, especially for 4, and thus is of theoretical as well as practical concern. Besides the different algorithms used to compute the Lyapunov spectrum, which is a secondary issue in this context, the number of points or iterations used in both references were substantially different. Wolf et al. [13] used 8192 (- 8.2 1s) points, where Aronson [14] used at least lo6 iterations.’ It is worth noting that in 1131, a rule of thumb was used to estimate the necessary number of points needed to achieve a ‘reasonable’ accuracy. This is the only instance, the author knows of where an attempt was made to devise an empirical stopping rule for this purpose. This rule says that about 10d - 30d points, where d is the Lyapunov dimension, are needed. For the Lorenz dynamics d = 2.07, and the required number of points estimated by this stopping rule is 102 - 932 points, which is about ten times less than the actual number of points used.” One way to alleviate this problem is to use an on-line adaptive version of the proposed KT stopping rule while estimating Lyapunov spectrums either from dynamical equations or time series. The following paragraph gives a rough sketch of this adaptive scheme. First, a warm-up stage is needed during which crude estimates of S and p, are generated, say either using a few hundred iterations, or Wolf et af. rule of thumb. These estimates with k' = 2, are used to compute the required initial number of iterations. A second stage of computation is initiated with the latter estimate as the needed number of iterations. +The number of iterations used for the Lorenz system was not given in [14], but for the H&on map it was reported to be 106. However, since the former dynamics have larger dimension, it is safe to assume that at least l@ iterations were needed for the reported numerical results. *As a side remark, and in the light of the results of this paper, the Wolf et aE. rule of thumb can only be used as a crude lower bound on the needed number of points or iterations.
Stopping rule for the Lyapunov spectrum
1225
During this second stage 6 and p are refined, and a new stopping time is estimated. This procedure is repeated until the successive stopping times become comparable. 5. SUMMARY AND CONCLUSIONS
In this paper an approximate stopping rule for the computation of the Lyapunov spectrum for linear random dynamical systems was proposed. Though the KT stopping rule is approximate it is asymptotically exact. The stopping rule involves all basic parameters that influence the dynamics under consideration. These parameters were explored for the special case of scalar dynamics, and were found to corraborate with theoretical as well as experimental understanding of complex dynamics. Even though the KT stopping rule is proposed for linear dynamics it was shown via an example that it could be used for nonlinear dynamics as well, after some basic manipulations. Nonetheless, when these manipulations are not easy or possible, an adaptive scheme of the KT stopping rule could be devised to estimate stopping times on-line. This adaptive scheme needs much work until it becomes practical. Some experimental results were presented in this work, but much more is still needed to fully understand the KT stopping rule and hopefully find clues to reduce its conservativeness. Acknowledgement--The
author acknowledges King Fahd University of Petroleum and Minerals for its support.
REFERENCES 1. L. Arnold and V. Wihstutz (eds), Lyupunov Exponents, Lecture Notes in Mathematics, No. 1186. Springer, NY (1986). 2. A. J. Lichtenberg and M. A. Liebermann, Regular and Chaotic Dynamics, Second Edition, Applied Mathematics Sciences, Vol. 38. Springer, NY (1992). 3. D. Talay, Approximation of upper Lyapunov exponents of bilinear stochastic differential systems, SIAM, J. on Num. Anul. 28(4), 1141 (1991). 4. J. P. Eckmann and D. Ruelle, Fundamental limitations for estimating dimensions and Lyapunov exponents in dynamical systems, IHES/P/90/11 (1990). 5. H. Furstenberg and R. Kesten, Product of random matrices. Ann. Math. Stats. 31,457 (1960). 6. J. Ezzine, On the stabilization and control of hybrid systems, Ph.D. thesis, Georgia Institute of Technology, (June 1989). 7. M. Katz and A. J. Thomasian, An exponential bound for functions of a Markov chain, Ann. Math. Stats. 31, 470 (1960). 8. J. Ezzine, On a practical stopping rule for the numerical computation of the Lyapunov spectrum, Proc. IEEE ZZth ACC, Chicago, IL (June 1992). 9. J. L. Doob, Stochastic Processes. John Wiley and Sons, NY (1953). 10. A. La&a and J. A. York, On the existence of invariant measures for piecewise monotonic transformations, Trans. Am. Math. Sot. 186, 481 (1973). 11. M. Bamsley, Fractals Everywhere. Academic Press, NY (1988). 12. G. Nicolis, Finite coarse-graining and Chapman-Kohnogorov equation in conservative dynamical systems, Chaos, Solitons & Fractak l(l), 25 (1991). 13. A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, Determining Lyapunov exponents from a time series, Physica D 16, 285 (1985). .. 14. J. W. Aronson, CHAOS, A SUN-based program for analyzing chaotic systems, Computers in Physics 408 Wm.