Numerical computation of Lyapunov exponent

Question

I'm trying to compute the Lyapunov exponent for a smooth continuous time dynamical system(say, $\dot{\bar{x}} = f(\bar x)$). I using the QR decomposition method. Here are the steps that I follow.

1. Choose some initial condition in the basin of the attractor. Call this $\bar v_0$. And have a blob (hypersphere, $U$) of unit radii around $\bar v_{0}$.
2. Iterate one time-step to get $\bar v_{1}$. The blob evolves as $U_{1} = D\bar f(\bar v_{0})\cdot U$. $D\bar{f}(\bar v_{0})$ being the Jacobian matrix.
3. Do a QR decomposition of $U_{1}$. $R$'s diag elements gives me the expansion/contraction rate. Store these lengths (in say J). And $Q$ is my new $U$.
4. Continue this process for say $n$ times.

The Lyapunov exponent(s) is given by, $\lambda = \frac1{k \cdot dt} \sum_{i=1}^{n} \log |J_{ii}|$. $dt$ is the time-step of the integrator.

I'm using this procedure to compute the Lyapunov exponent for the Lorenz system as a check. But unortunately this is not working. Is there anything wrong with the steps of the algorithm?

For example, for the canonical Lorenz system with the parameters, $\sigma = 10, r = 28, b = \frac83$, and starting with the initial condition $(1, 1, 1)$, I get all positive Lyapunov exponents! Them being, $25.6336, 20.1935, 16.76311$. Temporally they fluctuate for some time, and then settle on these values.

Here is the code that I'm using:

import numpy as np
from scipy.integrate import solve_ivp

#ODE system
def func(t, v, sigma, r, b):
    x, y, z = v #unpack the variables
    return [ sigma * (y - x), r * x - y - x * z, x * y - b * z ]

#Jacobian matrix
def JM(v, sigma, r, b):
    x, y, z = [k for k in v]
    return np.array([[-sigma, sigma, 0], [r - z, -1, -x], [y, x, -b]])

#initial parameters
sigma = 10
r = 28
b = 8/3

U = np.eye(3) #unit blob
v0 = np.ones(3) #initial condition
lyap = [] #empty list to store the lengths of the orthogonal axes

iters=10**3
dt=0.1
tf=iters * dt

#integrate the ODE system -- hopefully falls into an attractor
sol = solve_ivp(func, [0, tf], v0, t_eval=np.linspace(0, tf, iters), args=(sigma, r, b))
v_n = sol.y.T #transpose the solution

#do this for each iteration
for k in range(0, iters):
    v0 = v_n[k] #new v0 after iteration
    U_n = np.matmul(np.eye(3) + JM(v0, sigma, r, b) * dt, U)

    #do a Gram-Schmidt Orthogonalisation (GSO)
    Q, R = np.linalg.qr(U_n)
    lyap.append(np.log(abs(R.diagonal())))

    U = Q #new axes after iteration

[sum([lyap[k][j] for k in range(iters)]) / (dt * iters) for j in range(3)]

Note: The funution $\text{JM}$ gives the Jacobian. $\text{funs}$ gives the system of ODE (Lorenz system here).

Update: Updated code. As per Lutz Lehmann's answer.

Answer 1

It is not the Jacobian of $f$ that you need to use to propagate the local basis, but the Jacobian of the propagator of the numerical method. That is, if $$ v_1=v_0+\Phi_f(v_0,dt)dt $$ then $$ U_1=(I+D_v\Phi_f(v_0,dt)dt)U_0 $$ In a pinch you can of course replace the derivative $D_v\Phi_f(v_0,dt)=\frac{\partial \Phi_v}{\partial v}(v_0,dt)$ of the method step with the derivative $f'(v_0)=\frac{\partial f}{\partial v}(v_0)$ of the Euler step, as the Lyapunov exponents computed in this fashion will only be first order correct anyway, for whatever their exact definition may be. So then $$ U_1=(I+f'(v_0)dt)U_0 $$ The change in the code is just in this line

    U_n = np.matmul(I3+dt*JM(v0, sigma, r, b), U)

where I3 is the 3-dimensional identity matrix. The computed exponents are

dt = 0.1   : [ 4.44417446,  0.9529651 , -8.87227238]
dt = 0.05  : [  2.94558777,   1.16792384, -33.18463954]
dt = 0.02  : [  1.31725716,   1.05858356, -17.4156334 ]
dt = 0.01  : [  1.15509043,   0.45619985, -15.74654749]
dt = 0.001 : [  0.93408526,  -0.05527013, -14.57888732]

The above was computed by setting tf=100 and increasing the step number. As one can see, the result is highly unstable until the step size of the "Euler method" becomes small enough at $dt=0.02$. Choosing a different initial point does not change this much. The last two lines show the expected behavior of the result.

That there is a first-order error in the middle entry ios due to the fact that the solutions do not move in equal speed in a straight line. For a quasi-periodic dynamic that should even out, but only in the long term

A differential approach

Using the Euler method as propagator with extremely small step size $h\ll dt$, the procedure to compute the Lyapunov exponents is

• $v_{n+1}=v_n+hf(v_n)$,
• $\tilde U_{n+1}=(I+hf'(v_n))U_n$, with the general idea that $U_n$ is nearly an eigenbasis or right singular basis of $f'(v_n)$,
• $QR = \tilde U_{n+1}$, $R$ with positive diagonal,
• $U_{n+1} = Q$,
• $L^{yap}_{n+1}=L^{yap}_n+\ln(diag(R))$.

Now one can try to make all incremental updates proportional to $h$ to extract first the difference quotients and then in the limit differential equations. That is, consider only terms that are first order in $h$.

• $A = U_n^Tf'(v_n)U_n$ instead of $\tilde U_{n+1}=U_n(I+hA)$
• $QR = I+hA$ has the same $R$ as the QR decomposition of $\tilde U_{n+1}$, with $U_{n+1}=U_nQ$. We expect both factors to be close to the identity, $Q=I+hX+O(h^2)$, $R=I+hB$. In first order, this reduces $$(I+hX)(I+hB)=I+hA$$ to $$X+B=A=L_A+D_A+R_A.$$ One can read off $X=L_A-L_A^T$ as it is skew-symmetric and $B=D_A+R_A+L_A^T$ is upper triangular, with $L_A^T$ and $R_A$ strictly upper triangular.
• In the end, $U_{n+1}=U_n+hU_nX+O(h^2)$ and $L^{yap}_{n+1}=L^{yap}_n+\ln(diag(I+hB))=L^{yap}_n+h\,diag(B)+O(h^2)$ leads to the differential equations $$\dot U = UX, ~~\dot L^{yap} = diag(B)=diag(A)$$

This extended system can be implemented as

def diff_Lorenz(u):
    x,y,z = u
    f = [sigma * (y - x), r * x - y - x * z, x * y - b * z]
    Df = [[-sigma, sigma, 0], [r - z, -1, -x], [y, x, -b]]
    return np.array(f), np.array(Df)

def LEC_system(u):
    #x,y,z = u[:3]             # n=3
    U = u[3:12].reshape([3,3]) # size n square matrix, sub-array from n to n+n*n=n*(n+1)
    L = u[12:15]               # vector, sub-array from n*(n+1) to n*(n+1)+n=n*(n+2)
    f,Df = diff_Lorenz(u[:3])
    A = U.T.dot(Df.dot(U))
    dL = np.diag(A).copy();
    for i in range(3):
        A[i,i] = 0
        for j in range(i+1,3): A[i,j] = -A[j,i]
    dU = U.dot(A)
    return np.concatenate([f,dU.flatten(),dL])

u0 = np.ones(3)
U0 = np.identity(3)
L0 = np.zeros(3)
u0 = np.concatenate([u0, U0.flatten(), L0])
t = np.linspace(0,200,501)
u = odeint(lambda u,t:LEC_system(u),u0,t, hmax=0.05)
L = u[5:,12:15].T/t[5:]

plt.plot(t[5:],L.T)

and gives for $\lambda=L/t$ the plots

with the final values

[0.8501182905895157, -0.0028087381205551165, -14.5139762191356]

Answer 2

Just for completeness’ sake, here is a code that computes the evolution of the tangent vectors ($U$) and the regular system ($v$) in parallel, i.e., uses Benettin’s method. Instead of performing an Euler step (U_n = np.matmul(np.eye(3) + JM(v0, sigma, r, b) * dt, U)), the differential equation is expanded accordingly. This avoids having two integration schemes (Euler and whatever Solve IVP uses) in parallel.

import numpy as np
from scipy.integrate import solve_ivp

# ODE system
def func(t, v, sigma, r, b):
    x, y, z = v
    return [ sigma * (y - x), r * x - y - x * z, x * y - b * z ]

# Jacobian matrix
def JM(t, v, sigma, r, b):
    x, y, z = v
    return np.array([[-sigma, sigma, 0], [r - z, -1, -x], [y, x, -b]])

# parameters
sigma = 10
r = 28
b = 8/3

# dimension of the system (to keep things general)
n = 3

# number of Lyapunov exponents
n_lyap = 3

# (dis)assemble state and tangent vectors (or their derivatives)
# into one vector for integrator:
assemble = lambda v,U: np.hstack((v,U.flatten()))
disassemble = lambda state: (state[:n], state[n:].reshape(n,n_lyap))

def func_with_lyaps(t, state, *pars):
    v, U = disassemble(state)
    dv = func(t, v, *pars)
    dU = JM(t, v, *pars) @ U
    return assemble(dv, dU)

# initial states:
v = np.random.random(n)
U = np.random.random((n_lyap,n))

lyaps = [] # local Lyapunov exponents

dt = 1
iters = 1000

for _ in range(iters):
    sol = solve_ivp(
            func_with_lyaps,
            [0, dt],
            assemble(v,U),
            t_eval=[dt],
            args=(sigma, r, b),
            max_step=dt,
        )
    v,U = disassemble(sol.y.flatten())
    U,R = np.linalg.qr(U)
    lyaps.append(np.log(abs(R.diagonal()))/dt)

transient_steps = 100
# result:
print(*np.average(lyaps[transient_steps:],axis=0))

It should be straightforward to adapt this code to other systems.