Recently a new method for the computation of Lyapunov exponents that does not require rescaling and reorthogonalization was proposed ͓Rangarajan, Habib, and Ryne, Phys. Rev. Lett. 80, 3747 ͑1998͔͒. In this paper we make a detailed numerical comparison of the new method and a standard algorithm, as regards accuracy and efficiency, by applying them to some typical two-, three-, and four-dimensional systems. We find that in most cases there is reasonable agreement between the Lyapunov spectra obtained using the two algorithms. The CPU times required for computation are also comparable. However, in certain strongly chaotic cases, the new method was found to be either inefficient ͑taking a lot of CPU time for computation͒ or inaccurate. ͓S1063-651X͑99͒50907-0

K Ramasubramanian

M S Sriram

CiteSeerX

RAPID COMMUNICATIONS
PHYSICAL REVIEW E AUGUST 1999VOLUME 60, NUMBER 2Alternative algorithm for the computation of Lyapunov spectra of dynamical systems
K. Ramasubramanian and M. S. Sriram*
Department of Theoretical Physics, University of Madras, Guindy Campus, Chennai 600 025, India
~Received 25 November 1998!
Recently a new method for the computation of Lyapunov exponents that does not require rescaling and
reorthogonalization was proposed@Rangarajan, Habib, and Ryne, Phys. Rev. Lett.80, 3747 ~1998!#. In this
paper we make a detailed numerical comparison of the new method and a standard algorithm, as regards
accuracy and efficiency, by applying them to some typical two-, three-, and four-dimensional systems. We find
that in most cases there is reasonable agreement between the Lyapunov spectra obtained using the two
algorithms. The CPU times required for computation are also comparable. However, in certain strongly chaotic
cases, the new method was found to be either inefficient~taking a lot of CPU time for computation! or
inaccurate.@S1063-651X~99!50907-0#
PACS number~s!: 05.45.2a, 02.20.Qsal
e
t
te
he
fu
on
n
s
ne
W
d
n
ov
and
nov
it
-
d
se
-
ta-
his
ges
blesI. INTRODUCTION
Consider ann-dimensional continuous-time dynamic
system
dZ
dt
5F~Z,t !, ~1!
whereZ andF aren-dimensional vector fields. To determin
the n Lyapunov exponents of the system, corresponding
some initial conditionZ~0!, we have to find the long term
evolution of the axes of an infinitesimal sphere of sta
aroundZ~0!. For this, consider the tangent map given by t
set of equations
ddZ
dt
5DF•dZ, ~2!
whereDF is then3n Jacobian matrix with
DFi j 5
]Fi
]Zj
. ~3!
One of the standard methods used to determine the
Lyapunov spectrum due to Benettinet al. and Shimada and
Nagashima@1# uses the Gram-Schmidt reorthonormalizati
~GSR! procedure. An explicit source code for computatio
based on this procedure is given by Wolfet al. @2#. In this
method we have to integraten(n11) coupled equations, a
there aren equations for the fiducial trajectory in Eq.~1! and
n copies of then tangent map equations in Eq.~2!. We refer
to this method as the standard method in the following.
Recently, Rangarajan, Habib, and Ryne proposed a
algorithm for the computation of Lyapunov exponents@3#
based on the QR method@1# for the decomposition of the
tangent map. This does not require the GSR procedure.
summarize the essentials of this method below. The rea
can refer to@3# for details.
A solution of Eq.~2! can be formally written as
*Electronic address: tphysmu@imsc.ernet.inPRE 601063-651X/99/60~2!/1126~4!/$15.00o
s
ll
s
w
e
er
dZ~ t !5M „Z~ t !,t…dZ~0!, ~4!
whereM „Z(t),t… is the tangent map matrix whose evolutio
equation is easily seen to be
dM
dt
5DF•M . ~5!
The idea of the new method is to evaluate the Lyapun
exponents without using the vectorsd Z directly and conse-
quently without using the associated reorthogonalization
rescaling. For this one uses the fact thatM can be written as
M5QR, a product of an orthogonaln3n matrix Q and an
upper triangular matrixR with positive diagonal entries@4#.
Then it can be easily shown that
Q̃Q̇1ṘR215Q̃DFQ[S, ~6!
where the overdot denotes a time derivative. The Lyapu
exponentsl i are equal tos i /t in the limit t˜` wheres i
5 ln(Rii) @5#. ṘR
21 is also an upper triangular matrix and
is easily shown that the evolution equations fors i are con-
trolled by the diagonal elements ofS:
ṡ i5Sii , i 51,...,n. ~7!
Now Q, which is ann3n orthogonal matrix, is essen
tially the diagonalizing matrix for the tangent map flow an
is parametrized byn(n21)/2 angles~u’s!. Q̃Q̇ is an anti-
symmetric matrix and the evolution equations for the
angles can be obtained from the subdiagonal elements ofS in
Eq. ~6!. For n,4, we can work with any explicit represen
tation for Q. For n54, we employ a representation forQ
based on the well known fact thatSO(4);SO(3)3SO(3)
@6#. This simplifies the calculations and numerical compu
tions considerably. Hence, we have to solven(n13)/2
coupled equations to find the Lyapunov exponents in t
method, as there aren equations for the fiducial trajectory in
Eq. ~1!, n equations for the exponents in Eq.~7!, and n(n
21)/2 equations for the angles.
It might be thought that the new method has advanta
over the standard methods, as a minimal number of variaR1126 © 1999 The American Physical Society
lim
fo
ea
ta
cc
u
m
ex
a
a
u
u
a
.
s
-
n
a
r
te
of
e
.
e
the
nal
en
ta-
size
-
ng
ase
by
in
lly
em,
-
s,
en
m,
died
nts
f
e
r a
wo
po-
rs,
re-
ies
RAPID COMMUNICATIONS
PRE 60 R1127ALTERNATIVE ALGORITHM FOR THE COMPUTATION . . .is used and rescaling and reorthogonalization are also e
nated. However, in this method the evolution equations
the angles and Lyapunov exponents are highly nonlin
involving sines and cosines of the angles, whereas the s
dard method uses the linearized equations ford Z directly.
Hence, there is a need to compare the efficiency and a
racy of this method with a standard method. That is the s
ject of the present investigation. Here, we consider so
typical nonlinear systems of physical interest withn52, 3,
and 4. The driven Van der Pol oscillator is taken as an
ample of a two-dimensional system, whereas the stand
Lorenz system is chosen forn53. Forn54, we consider the
coupled quartic oscillators and anisotropic Kepler problem
examples of conservative Hamiltonian systems and Ro¨ssler
hyperchaos system as an example of a dissipative system
all these cases, the full Lyapunov spectrum is computed
ing both methods. The time of integration is chosen to ens
reasonable convergence of the Lyapunov exponents.
II. COMPARISON OF THE TWO METHODS
In this section, we take up the following systems for
detailed comparison of the two methods.
~i! Driven Van der Pol oscillator (n52).
d
dt S z1z2D5S z22d~12z12!z22z11b cosvt,D , ~8!
whereb andd are parameters andv is the driven frequency
In our numerical work we have chosend525.0, b55.0,
andv52.47 as the parameter values.
~ii ! Lorenz system (n53).
d
dt S z1z2
z3
D 5S s~z22z1!z1~r2z3!2z2
z1z22bz3
D . ~9!
This system is too well known to require any further discu
sion. For computations we sets510.0, r528.0, andb
5 8/3.
~iii ! Coupled quartic oscillators (n54). This is a con-
servative system and the Hamiltonian is given by
H5
z3
2
2
1
z4
2
2
1z1
41z2
41az1
2z2
2 , ~10!
wherez1 andz2 are the canonical coordinates,z3 andz4 are
the corresponding momenta, anda is a parameter. The equa
tions of motion are readily obtained from the Hamiltonia
This system is known to be integrable fora50, 2 and 6.
~iv! Anisotropic Kepler problem (n54). The Hamil-
tonian of this system is given by
H5
pr
2
2
1g
pz
2
2
2
e2
Ar21z2
, ~11!
whereg is a number.
The Hamiltonian given above describes the motion of
electron in the Coloumb field in an anisotropic crystal, whe
its effective mass along thex-y plane andz direction are
different.g51 corresponds to the isotropic case and is ini-
r
r,
n-
u-
b-
e
-
rd
s
. In
s-
re
-
.
n
e
-
grable. WhengÞ1, the system is nonintegrable. Because
the singularity atr5z50, the Hamiltonian in the above
form is hardly suitable for numerical integration. For this w
choosez15Ar1z andz25Ar2z as the canonical variables
We can find the corresponding canonical momentaz3 ndz4
in terms ofpr and pz . We also use a reparametrized tim
variablet defined bydt5dt(z1
21z2
2).
The original Hamiltonian with old variables and energyE
corresponds to the following Hamiltonian withH852 in
terms of the new variables@7#:
H8525
1
2
~z3
21z4
2!2E~z1
21z2
2!1~g21!
~z1z32z2z4!
2
2~z1
21z2
2!
.
~12!
The equations of motion can be easily obtained from
above Hamiltonian. We have choseng50.61 for computa-
tional purposes.
~v! Rössler hyperchaos system (n54). This is a dissi-
pative system and an extension of the three-dimensio
Rössler attractor@8,9#. It is described by the equations
d
dt S z1z2z3
z4
D 5S 2~z21z3!z11az21z4b1z1z3
cz42dz3
D , ~13!
wherea, b, c, andd are parameters whose values are tak
to be 0.25, 3.0, 0.05, and 0.5, respectively, for our compu
tions.
For all the systems, we have used a variable step-
Runge-Kutta routine~RKQC! for integration with an error
tolerancee;1026– 1028. All the computations were per
formed on a DEC Alpha based workstation runni
OpenVMS. We also noted the CPU time taken for each c
with either of the algorithms. This is the actual time taken
the CPU to accomplish a specific process~independent of the
other processes running in the system!. The details of the
comparison between the two methods are summarized
Table I.
It may be noticed that the two methods yield essentia
the same Lyapunov spectrum. For any dynamical syst
one of the Lyapunov exponents has to be zero~correspond-
ing to the difference vectordz lying along the trajectory
itself!. For the Lorenz system, the Ro¨ssler hyperchaos sys
tem ~both dissipative!, and the coupled quartic oscillator
this condition is satisfied by both algorithms. For the driv
Van der Pol oscillator and the anisotropic Kepler proble
both methods fail the test. This aspect needs to be stu
further. For the coupled quartic oscillators, all the expone
should be zero corresponding to the integrable case oa
56. This is indeed satisfied by both algorithms. In Fig. 1 w
give plots of Lyapunov exponents as functions of time fo
typical case. Again, there is little difference between the t
algorithms as far as the convergence of the Lyapunov ex
nents is concerned.
However, for the system of coupled quartic oscillato
the CPU time is abnormally high for the new method, cor
sponding to the nonintegrable case ofa58. This is true for
both small and large energies. For large energ
(;25 000), since the energy varied by;15 when we used
th
ne
ow
o
io
-
a
ce
us
le
e
en
77
nt
ua-
in-
the
ese
are
ity
on-
ven
of
ated
nal
uire
ser
eral
we
cal
nd
ng a
son-
two
lso
ns,
he
ms
d
lgo-
RAPID COMMUNICATIONS
R1128 PRE 60K. RAMASUBRAMANIAN AND M. S. SRIRAMthe RKQC routine, we also used a symplectic procedure
eliminates secular variations in the energy. With this routi
the CPU times were nearly the same for both methods. H
ever, the new method yields poor results for the Lyapun
spectrum. For instance, corresponding to the initial condit
z157.0, z257.0, z355.0, andz454.0, the Lyapunov spec
tra computed using the new and the standard methods
~1.5506, 0.3254,20.3261,21.5499) and~1.5205, 0.0001,
20.0001,21.5205), respectively.
We also compared the new method with another pro
dure for computing Lyapunov spectra with continuo
Gram-Schmidt orthonormalization@10#. Here the number of
equations that need to be integrated to obtain the comp
spectrum isn(n12), as compared ton(n11) equations in
the standard method andn(n13)/2 equations in the new
method, wheren is the order of the system. The CPU tim
for this method, corresponding to the initial conditions giv
in Table I for a58, is 7658.57 s, as compared to 39 012.
TABLE I. Comparison of the two methods for some syste
with n52, 3, and 4. The values given in parentheses correspon
the standard method.
System
Initial
condition
Lyapunov spectrum,
sum (s), and
CPU time (t) in sec
Driven 0.0981~0.0987!
van der Pol z1521.0 26.8400 (26.8411)
Oscillator z251.0 s526.7419 (26.7424)
(n52) t5519.22~825.56!
0.9056~0.9051!
Lorenz z150.0 0.0000~0.0000!
system z251.0 214.5723 (214.5718)
(n53) z350.0 s5213.6667 (213.6667)
t52394.30~1668.68!
0.1360~0.1332!
Anisotropic z151.0 0.0831~0.0832!
Kepler z252.0 20.0833 (20.0833)
problem z351.0 20.1357 (20.1331)
(n54) z450.5 s520.0000 (20.0000)
t5350.18~201.04!
0.1128~0.1121!
z15220.0 0.0214~0.0196!
Rössler z250.0 20.0000 (20.0000)
hyperchaos z350.0 224.7527 (225.1886)
(n54) z4515.0 s5224.6185 (225.0568)
t51527.58~5594.99!
0.0001~0.0001!
z150.8 0.0001~0.0001!
Coupled z250.5 20.0001 (20.0001)
quartic oscr. z351.0 20.0001 (20.0001)
(n54, a56) z451.3 s50.0000~0.0000!
t5803.49~492.09!
0.1806~0.1738!
z150.8 0.0001~0.0001!
Coupled z250.5 20.0001 (20.0001)
quartic oscr. z351.0 20.1806 (20.1738)
(n54, a58) z451.3 s50.0000~0.0000!
t539012.77~855.64!at
,
-
v
n
re
-
te
with the new method~with hardly any difference in the
Lyapunov spectrum!.
In the standard method, as well as in Ref.@10#, after solv-
ing for the fiducial trajectory, the equations for the tange
flow are linearized equations. In the new method, these eq
tions are replaced by the equations for the angles determ
ing the principal axes or the bases associated with
Lyapunov spectrum and the Lyapunov exponents. Th
equations involving sines and cosines of the angles
highly nonlinear. For dissipative systems this nonlinear
does not pose a problem. However, in many cases, this n
linearity renders the new method less efficient, and can e
lead to inaccuracies in strongly chaotic situations.
III. CONCLUSIONS
In a recently proposed new method for the computation
Lyapunov exponents, the Lyapunov exponents are calcul
directly, so to say, by utilizing representations of orthogo
matrices applied to the tangent map. Since it does not req
renormalization or reorthogonalization and requires a les
number of equations, it has been claimed that it has sev
advantages over existing methods. To test this claim,
have computed the full Lyapunov spectrum of some typi
nonlinear systems with two, three, and four variables a
made a detailed comparison with the results obtained usi
standard algorithm. For dissipative systems, there is rea
able agreement between the spectra obtained using the
algorithms. The CPU time taken for the computation is a
comparable. However, in certain strongly chaotic situatio
the new algorithm could lead to inaccuracies in t
to
FIG. 1. Plots of the Lyapunov exponent for the Ro¨ssler hyper-
chaos system.~a! Highest exponentl1 , ~b! lowest exponentl4 .
The thick and thin lines correspond to the new and standard a
rithms, respectively.
a
b
ab
fa
-
a
a
e
lgo-
-
RAPID COMMUNICATIONS
PRE 60 R1129ALTERNATIVE ALGORITHM FOR THE COMPUTATION . . .Lyapunov spectrum, when one uses fixed step-size integr
routines, however small the step size may be. This could
remedied by a variable step-size routine with a reason
value of error tolerance. But this makes the new method
less efficient~taking abnormally long times for computa
tion!, compared to the standard algorithm. This could be
tributed to the highly nonlinear nature of the evolution equ
tions for the tangent map inherent in this method. Howev-
g.
,
tor
e
le
r
t-
-
r,
the proposed new method is still useful as an alternate a
rithm for the computation of Lyapunov spectra.
ACKNOWLEDGMENT
One of the authors~K.R.! thanks the Council for Scien
tific and Industrial Research, India for financial support.@1# G. Benettin, L. Galgani, A. Giorgilli, and J. M. Strelcyn, Mec
canica 15, 9 ~1980!; I. Shimada and T. Nagashima, Pro
Theor. Phys.61, 1605~1979!.
@2# A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano
Physica D16, 285 ~1985!.
@3# G. Rangarajan, S. Habib, and R. D. Ryne, Phys. Rev. Lett.80,
3747 ~1998!.
@4# G. H. Golub and C. F. Van Loan,Matrix Computations, 3rd
ed. ~John Hopkins University Press, Baltimore, 1996!.
@5# K. Geist, U. Parlitz, and W. Lauterborn, Prog. Theor. Phys.83,875 ~1990!.
@6# See, for instance, B. G. Wybourne,Classical Groups for
Physicists~Wiley, New York, 1974!.
@7# M. Kuwata, A. Harada, and H. Hasegawa, J. Phys. A23, 3227
~1990!.
@8# O. E. Rössler, Phys. Lett.71, 155 ~1979!.
@9# H. Yoshida, Phys. Lett. A150, 262 ~1990!.
@10# I. Goldhirsch, P. L. Sulem, and S. A. Orszag, Physica D27,
311 ~1987!; F. Christiansen and H. H. Rugh, Nonlinearity10,
1063 ~1997!.


Alternative algorithm for the computation of Lyapunov spectra of dynamical systems

http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=87CA91796FE34CC6D4B4CFCFCC06C608?doi=10.1.1.1056.7951&rep=rep1&type=pdf

Alternative algorithm for the computation of Lyapunov spectra of dynamical systems

Abstract

Similar works

Full text

Available Versions

CiteSeerX