Notes on How to Numerically Cal-
culate the Maximum Lyapunov
Exponent
These notes, primarily for my own reference,
briefly describe how my programs calculate
the maximum Lyapunov exponent. There is
nothing extraordinary in how I do it the
method is described in ample detail by well-
known references, e.g. Wolf et al. (1985,
Physica D 16, 285) and Benettin et al.
(1976, Phys. Rev. A 14, 2338).
Consider two orbits, a "reference" orbit and
a "test" orbit, separated at time t
0
by a small
phase space distance d
0
. We will use the test
orbit as a means of calculating the value of
the maximum Lyapunov exponent. Under
evolution of the equations of motion, the
two orbits may (or may not) separate. If the
motion is chaotic, the orbits will, by defini-
tion, separate at an exponential rate. The
maximum Lyapunov exponent
is a measure
of this rate of separation:
(1)
k
=
t
dº
lim
1
t
-
t
0
ln
d
(
t
)
d
0
Hence, in the limit of infinite time,
(2)
k
=
t
dº
lim
1
t
-
t
0
ln
d
(
t
)
d
0
In practice, we cannot afford the luxury of
infinitely long integrations, so we instead
calculate the instantaneous maximum Lyapu-
nov exponent
(3)
k
(
t
) =
1
t
-
t
0
ln
d
(
t
)
d
0
and, ideally, wait long enough for
(t) to
settle to approximately its asymptotic value,
if indeed it is non zero for the orbit of
interest. A simple method of calculating
(t)
is shown in section 1. Another practical
problem is that, for chaotic orbits, the dis-
tance between reference and test particles,
d(t), quickly saturates. Hence we must peri-
odically renormalize the orbit separation.
This is shown in section 2. Section 3 pre-
sents the derivation of eq. (7) in more detail.
1. Exponent Calculation
We will leave the reference orbit alone and
rescale the test orbit whenever the separa-
tion d(t) has passed beyond a threshold value
D. It is important that D be set small
enough that it is still in the linear regime
(i.e., the regime in which the linearized equa-
tions of motion are an accurate description).
Define a rescaling parameter:
(4)
a
1
h
d
(
t
1
)
d
(
t
0
)
where t
1
is the time at which
. Then
d
(
t
)
m D
we can write
(5)
k
1
=
1
t
1
-
t
0
ln
d
1
d
0
=
1
t
1
-
t
0
ln a
1
where
and
. At this
k
i
h k
(
t
i
)
d
i
h d
(
t
i
)
point, the test orbit is then rescaled, as
shown in section 2. Similarly, for successive
threshold crossings and subsequent rescal-
ings, we have
f:\dynamics\text\lyapcalc.lwp
Maximum Lyapunov Exponent Calculation
7:07pm September 30, 1995
1
