Computing Lyapunov Exponents

Lyapunov exponents describe the average rate at which nearby trajectories of a model separate or converge. run_lyapunov runs COPASI’s Lyapunov Exponents task on the current model and returns three values: the exponents, their sum, and the average divergence of the model equations.

The sign of an exponent has a direct meaning:

• positive: nearby trajectories diverge exponentially, the signature of chaos

• about zero: motion is neutral in that direction, which is typical along a limit cycle

• negative: nearby trajectories converge, as when a cycle or steady state is attracting

This example uses the Brusselator, a chemical oscillator. Only X and Y are dynamic species, so the state space is two-dimensional and two exponents are the complete spectrum. With the parameters in the example, the model settles onto a stable oscillation.

from basico import *

Load the example. The name is matched against the models shipped with basico, so bruss selects the Brusselator.

model = load_example('bruss')

num_exponents sets how many exponents are computed. The default is 3; this model only has two independent directions, so we ask for 2. The average divergence is included by default.

exponents, sum_of_exponents, divergence = run_lyapunov(num_exponents=2)
print('exponents:        ', exponents)
print('sum of exponents: ', sum_of_exponents)
print('average divergence:', divergence)

exponents:         [-0.0007025776472087997, -2.5817449635993714]
sum of exponents:  -2.58244754124658
average divergence: -2.6300694998749843

The first exponent is about zero and the second is about -2.6. That spectrum belongs to a stable limit cycle: one direction follows the orbit, and perturbations away from it decay. A positive exponent would indicate chaos, and a spectrum of strictly negative exponents would indicate a stable steady state.

When every exponent of the system has been computed, their sum equals the average divergence, the time average of the trace of the Jacobian. Both are about -2.6 here; the remaining difference is numerical.

The same numbers can be read again with get_lyapunov_exponents, without repeating the integration. On larger models, calculate_divergence=False skips the divergence, and start_averaging_after sets the time at which averaging begins so that an initial transient is left out.