3-lyapunov指数计算方法_wolf算法_

Physica 16D (1985)285-317 North-Holland, Amsterdam DETERMINING LYAPUNOV EXPONENTS FROM A TIME SERIES Alan WOLF~-, Jack B. SWIFT, Harry L. SWINNEY and John A. VASTANO Department of Physics, University of Texas, Austin, Texas 78712, USA Received 18 October 1984 We present the first algorithms that allow the estimation of non-negative Lyapunov exponents from an experimental time series. Lyapunov exponents, which provide a qualitative and quantitative characterization of dynamical behavior, are related to the exponentially fast divergence or convergence of nearby orbits in phase space. A system with one or more positive Lyapunov exponents is defined to be chaotic. Our method is rooted conceptually in a previously developed technique that could only be applied to analytically defined model systems: we monitor the long-term growth rate of small volume elements in an attractor. The method is tested on model systems with known Lyapunov spectra, and applied to data for the Belousov-Zhabotinskii reaction and Couette-Taylor flow. Contents 1. Introduction 2. The Lyapunov spectrum defined 3. Calculation of Lyapunov spectra from differential equations 4. An approach to spectral estimation for experimental data 5. Spectral algorithm implementation* 6. Implementation details* 7. Data requirements and noise* 8. Results 9. Conclusions Appendices* A. Lyapunov spectrum program for systems of differential equations B. Fixed evolution time program for ~'1 1. Introduction Convincing evidence for deterministic chaos has come from a variety of recent experiments [1-6] on dissipative nonlinear systems; therefore, the question of detecting and quantifying chaos has become an important one. Here we consider the spectrum of Lyapunov exponents [7-10], which has proven to be the most useful dynamical di- agnostic for chaotic systems. Lyapunov exponents are the average exponential rates of divergence or tPresent address: The Cooper Union, School of Engineering, N.Y., NY 10003, USA. *The reader may wish to skip the starred sections at a first reading. convergence of nearby orbits in phase space. Since nearby orbits correspond to nearly identical states, exponential orbital divergence means that systems whose initial differences we may not be able to resolve will soon behave quite differently-predic- tive ability is rapidly lost. Any system containing at least one positive Lyapunov exponent is defined to be chaotic, with the magnitude of the exponent reflecting the time scale on which system dynamics become unpredictable [10]. For systems whose equations of motion are ex- plicitly known there is a straightforward technique [8, 9] for computing a complete Lyapunov spec- trum. This method cannot be applied directly to experimental data for reasons that will be dis- cussed later. We will describe a technique which for the first time yields estimates of the non-nega- tive Lyapunov exponents from finite amounts of experimental data. A less general procedure [6, 11-14] for estimat- ing only the dominant Lyapunov exponent in ex- perimental systems has been used for some time. This technique is limited to systems where a well- defined one-dimensional (l-D) map can be re- covered. The technique is numerically unstable and the literature contains several examples of its improper application to experimental data. A dis- cussion of the 1-D map calculation may be found 0167-2789/85/$03.30 © Elsevier Science Publishers (North-H01!and Physics Publishing Division) 286 A. Wolf et al. / Determining Lyapunov exponents from a time series in ref. 13. In ref. 2 we presented an unusually robust 1-D map exponent calculation for experi- mental data obtained from a chemical reaction. Experimental data inevitably contain external noise due to environmental fluctuations and limited experimental resolution. In the limit of an infinite amount of noise-free data our approach would yield Lyapunov exponents by definition. Our abil- ity to obtain good spectral estimates from experi- mental data depends on the quantity and quality of the data as well as on the complexity of the dynamical system. We have tested our method on model dynamical systems with known spectra and applied it to experimental data for chemical [2, 13] and hydrodynamic [3] strange attractors. Although the work of characterizing chaotic data is still in its infancy, there have been many ap- proaches to quantifying chaos, e.g., fractal power spectra [15], entropy [16-18, 3], and fractal dimen- sion [proposed in ref. 19, used in ref. 3-5, 20, 21]. We have tested many of these algorithms on both model and experimental data, and despite the claims of their proponents we have found that these approaches often fail to characterize chaotic data. In particular, parameter independence, the amount of data required, and the stability of re- suits with respect to external noise have rarely been examined thoroughly. The spectrum of Lyapunov exponents will be defined and discussed in section 2. This section includes table I which summarizes the model sys- tems that are used in this paper. Section 3 is a review of the calculation of the complete spectrum of exponents for systems in which the defining differential equations are known. Appendix A con- tains Fortran code for this calculation, which to our knowledge has not been published elsewhere. In section 4, an outline of our approach to estimat- ing the non-negative portion of the Lyapunov exponent spectrum is presented. In section 5 we describe the algorithms for estimating the two largest exponents. A Fortran program for de- termining the largest exponent is contained in appendix B. Our algorithm requires input parame- ters whose selection is discussed in section 6. Sec- tion 7 concerns sources of error in the calculations and the quality and quantity of data required for accurate exponent estimation. Our method is ap- plied to model systems and experimental data in section 8, and the conclusions are given in section 9. 2. The Lyapunov spectrum defined We now define [8, 9] the spectrum of Lyapunov exponents in the manner most relevant to spectral calculations. Given a continuous dynamical sys- tem in an n-dimensional phase space, we monitor the long-term evolution of an infinitesimal n-sphere of initial conditions; the sphere will become an n-ellipsoid due to the locally deforming nature of the flow. The ith one-dimensional Lyapunov expo- nent is then defined in terms of the length of the ellipsoidal principal axis pi(t) : h~ = lim 1 log 2 pc( t ) t--,oo t pc(O)' (1) where the )h are ordered from largest to smallestt. Thus the Lyapunov exponents are related to the expanding or contracting nature of different direc- tions in phase space. Since the orientation of the ellipsoid changes continuously as it evolves, the directions associated with a given exponent vary in a complicated way through the attractor. One can- not, therefore, speak of a well-defined direction associated with a given exponent. Notice that the linear extent of the ellipsoid grows as 2 htt, the area defined by the first two principal axes grows as 2 (x~*x2)t, the volume de- fined by the first three principal axes grows as 2 (x'+x2+x~)t, and so on. This property yields another definition of the spectrum of exponents: tWhile the existence of this limit has been questioned [8, 9, 22], the fact is that the orbital divergence of any data set may be quantified. Even if the limit does not exist for the underlying system, or cannot be approached due to having finite amounts of noisy data, Lyapunov exponent estimates could still provide a useful characterization of a given data set. (See section 7.1.) A. Wolf et aL / Determining Lyapunov exponents from a time series 287 the sum of the first j exponents is defined by the long term exponential growth rate of a j-volume element. This alternate definition will provide the basis of our spectral technique for experimental data. Any continuous time-dependent dynamical sys- tem without a fixed point will have at least one zero exponent [22], corresponding to the slowly changing magnitude of a principal axis tangent to the flow. Axes that are on the average expanding (contracting) correspond to positive (negative) ex- ponents. The sum of the Lyapunov exponents is the time-averaged divergence of the phase space velocity; hence any dissipative dynamical system will have at least one negative exponent, the sum of all of the exponents is negative, and the post- transient motion of trajectories will occur on a zero volume limit set, an attractor. The exponential expansion indicated by a posi- tive Lyapunov exponent is incompatible with mo- tion on a bounded attractor unless some sort of folding process merges widely separated trajecto- ries. Each positive exponent reflects a "direction" in which the system experiences the repeated stretching and folding that decorrelates nearby states on the attractor. Therefore, the long-term behavior of an initial condition that is specified with any uncertainty cannot be predicted; this is chaos. An attractor for a dissipatiVe system with one or more positive Lyapunov exponents is said to be "strange" or "chaotic". The signs of the Lyapunov exponents provide a qualitative picture of a system's dynamics. One- dimensional maps are characterized by a single Lyapunov exponent which is positive for chaos, zero for a marginally stable orbit, and negative for a periodic orbit. In a three-dimensional continuous dissipative dynamical system the only possible spectra, and the attractors they describe, are as follows: (+ ,0 , - ) , a strange attractor; (0,0,-) , a two-toms; (0, - , - ) , a limit cycle; and ( - , - , - ) , a fixed point. Fig. 1 illustrates the expanding, "slower than exponential," and contracting char- acter of the flow for a three,dimensional system, the Lorenz model [23]. (All of the model systems that we will discuss are defined in table I.) Since Lyapunov exponents involve long-time averaged behavior, the short segments of the trajectories shown in the figure cannot be expected to accu- rately characterize the positive, zero, and negative exponents; nevertheless, the three distinct types of behavior are clear. In a continuous four-dimen- sional dissipative system there are three possible types of strange attractors: their Lyapunov spectra are (+, + ,0 , - ) , (+ ,0 ,0 , - ) , and (+,0 , - , - ) . An example of the first type is Rossler's hyper- chaos attractor [24] (see table I). For a given system a change in parameters will generally change the Lyapunov spectrum and may also change both the type of spectrum and type of attractor. The magnitudes of the Lyapunov exponents quantify an attractor's dynamics in information theoretic terms. The exponents measure the rate at which system processes create or destroy informa- tion [10]; thus the exponents are expressed in bits of information/s or bits/orbit for a continuous system and bits/iteration for a discrete system. For example, in the Lorenz attractor the positive exponent has a magnitude of 2.16 bits/s (for the parameter values shown in table I). Hence if an initial point were specified with an accuracy of one part per million (20 bits), the future behavior could not be predicted after about 9 s [20 bits/(2.16 bits/s)], corresponding to about 20 orbits. After this time the small initial uncertainty will essen- tially cover the entire attractor, reflecting 20 bits of new information that can be gained from an ad: ditional measurement of the system. This new information arises from scales smaller than our initial uncertainty and results in an inability to specify the state of the system except to say that it is somewhere on the attractor. This process is sometimes called an information gain- reflecting new information from the heat bath, and some- times is called an information loss-bits shifted out of a phase space variable "register" when bits from the heat bath are shifted in. The average rate at which information con- tained in transients is lost can be determined from 288 A. Wolf et al. / Determining Lyapunov exponents from a time series • ° • o • ° ° ° . ° • ° • ° ° • • " . . t "." : . . . . ".• . • . • . . - . . . . . . . "• . . . . . . , : .~ "'.•'. •... . . . . , - - . : - : . : : - : . . . . . . . . . • . • .,:'..~..--..:~::.-.:..:'..:..:.. .. .. ,~ . , . : . . : ' . . " • - . . . • . . s tart •.'~"." .- • " " " : ' " : " : : " ' " N~ " ~" "" ~'"" """ • ".':':"'( ( ,~ ' . "~,m ~" " " : ' " ' " ' " " . • - ' . : . . " "%V4" ; ' : ' " " ' - " . ' . . . " . . . . ~ " . • - . . , . . : /~ . • : • . . . . " . . . . . . . . . . . . . i IIIIIl[lll t ime-~ • • • • • . . ° • . ° . ° ° • • - . . o , - : . . . -. . • . • • . . . . . . . , : , , . : : . " .• . • . (b ) . . . ' , .<'."~::.:. ' : . . "... "-.. ,...~.'~:'-:,~:..~ "r.. : - : • . ~ . . . . . .:'~.?'-,',;~ "x -~ i l~ I - " . . . " "~,~ ' . :~ '~" , -~,~ '~xxx ~.~' - .'.• .: • " :::,.!...":"""-...- • " • ". " : : : , . f ,~ ,_ , , ,~ .~ - , . - , , . , . - . . . : . . . . . . • $ :L~. . . . . ~- "~. ' . . . _ _ '•• .~. . . . : - . : . . ." . . . . : : . : . . ; : • • . - ' .~ . . .~ . ; . ; . , . . . . . . . . . . . • . . . . . . . . . . . . , " " . : , v ' : . :~ . : " "~. ' . " ' " . . . . . " " : ~star t [ l , ,m,, , , l , , , , , l l l i , , l l l l i , ,d l l t ime - -~ • . . . - . - . . . . . :; ..- ---~..- : .. • . . . ;;.~.....~.::........... s,,,, ~.~,.~;,..':..::.~:"~-:.v • : ' " " , ' , . . . :~ . ' . -2~'W~'" .~- . . . ".".'.. " ' : : . . . ' : " ' . : - - " ". "- : " • . , , • " •~, '•••• .• •" • . . . . . . . . . t ime Fig. 1. The short term evolution of the separation vector between three carefully chosen pairs of nearby points is shown for the Lorenz attractor, a) An expanding direction (~1 > 0); b) a "slower than exponential" direction (~'2 = 0); C) a contracting direction (X3 < 0) . the negative exponents• The asymptotic decay of a perturbation to the attractor is governed by the least negative exponent, which should therefore be the easiest of the negative exponents to estimatet. tWe have been quite successful with an algorithm for de- termiuing the dominant (smallest magnitude) negative expo- nent from pseudo-experimental data (a single time series ex- tracted from the solution of a model system and treated as an experimental observable) for systems that are nearly integer- dimensional. Unfortunately, our approach, which involves mea- suring the mean decay rate of many induced perturbations of the dynamical system, is unlikely to work on many experimen- tal systems. There are several fundamental problems with the calculation of negative exponents from experimental data, but For the Lorenz attractor the negative exponent is so large that a perturbed orbit typically becomes indistinguishable from the attractor, by "eye", in less than one mean orbital period (see fig. 1). of greatest importance is that post-transient data may not contain resolvable negative exponent information and per- turbed data must refl~t properties of the unperturbed system, that is, perturbations must only change the state of the system (current values of the dynamical variables). The response of a physical system to a non-delta function perturbation is difficult to interpret, as an orbit separating from the attractor may reflect either a locally repelling region of the attractor (a positive contribution to the negative exponent) or the finite rise time of the perturbation. A. Wolf et al. / Determining Lyapunov exponents from a time series 289 Table I The model systems considered in this paper and their Lyapunov spectra and dimensions as computed from the equations of motion Lyapunov Lyapunov System Parameter spectrum dimension* values (bits/s)t H~non: [25] ~1 = 0.603 X. +1 = 1 - aX;. + Yn { b = 1.4 h 2 = - 2.34 Y. + 1 = bX. = 0.3 (bits/iter.) Rossler-chaos: [26] ) ( = - (Y + Z) [ a = 0 .15 )k 1 = 0.13 ) '= X+ aY I b = 0.20 ~2 =0.00 = b + Z(X- c) c = 10.0 h 3 = - 14.1 Lorenz: [23] ) (= o (Y - X) [ o = 16.0 h 1 = 2.16 ~'= X( R - Z ) - Y I R=45.92 X 2 =0.00 = XY - bZ b = 4.0 ;k 3 = - 32.4 Rossler-hyperchaos: [24] Jr'= - (Y+ Z) ( a = 0.25 A t = 0.16 ) '= X+ aY+ W [ b= 3.0 X 2 =0.03 = b + XZ | c = 0.05 h 3 = 0.00 if" = cW - dZ k d = 0.5 h4 = - 39.0 Mackey-Glass: [27] ( a = 0.2 h t = 6.30E-3 j ( = aX( t + s ) - bX(t) / b = 0.1 )~2 = 2.62E-3 1 + [ X ( t + s) ] c ) c = 10.0 IX31 < 8.0E-6 s = 31.8 )'4 = - 1.39E-2 1.26 2.01 2.07 3.005 3.64 tA mean orbital period is well defined for Rossler chaos (6.07 seconds) and for hyperchaos (5.16 seconds) for the parameter values used here. For the Lorenz attractor a characteristic time (see footnote- section 3) is about 0.5 seconds. Spectra were computed for each system with the code in appendix A. ~As defined in eq. (2). The Lyapunov spectrum is closely related to the fractional dimension of the associated strange at- tractor. There are a number [19] of different frac- tional-dimension-like quantities, including the fractal dimension, information dimension, and the correlation exponent; the difference between them is often small. It has been conjectured by Kaplan and Yorke [28, 29] that the information dimension d r is related to the Lyapunov spectrum by the equation Ei-- 1~i df=J+ I?~j+il ' (2) where j is defined by the condition that j j+ l E)~i> 0 and EX,