# Numerical calculation of Lyapunov exponents using SciPy's built-in solve_ivp

I have previously successfully implemented the QR decomposition method in MATLAB to calculate Lyapunov exponents for Lorenz equations. See here.

This method integrates the stacked system, i.e. the system equation + the variational equation and use one common integrator to integrate them at the same time steps.

What I am trying to do now is to use any of the integrators built in to Scipy's solve_ivp to do this. In order to do this, I would need to switch out the time step method in scipy.integrate._ivp.rk._step_impl_qr

Here is what I have:

import numpy as np
import scipy
import matplotlib.pyplot as plt

import scipy.integrate._ivp.rk as rk

### Define a new function for '_step_impl' this called '_step_impl_qr'
### The bits I changed from Scipy's source code are in between
### ###################################
###      This is what I changed
### ###################################

lyap = []
ndim = 3
def _step_impl_qr(self):
t = self.t
y = self.y

max_step = self.max_step
rtol = self.rtol
atol = self.atol

min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)

if self.h_abs > max_step:
h_abs = max_step
elif self.h_abs < min_step:
h_abs = min_step
else:
h_abs = self.h_abs

step_accepted = False
step_rejected = False

while not step_accepted:
if h_abs < min_step:
return False, self.TOO_SMALL_STEP

h = h_abs * self.direction
t_new = t + h

if self.direction * (t_new - self.t_bound) > 0:
t_new = self.t_bound

h = t_new - t
h_abs = np.abs(h)

y_new, f_new = rk.rk_step(self.fun, t, y, self.f, h, self.A,
self.B, self.C, self.K)
scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol
error_norm = self._estimate_error_norm(self.K, h, scale)

if error_norm < 1:
if error_norm == 0:
factor = rk.MAX_FACTOR
else:
factor = min(rk.MAX_FACTOR,
rk.SAFETY * error_norm ** self.error_exponent)

if step_rejected:
factor = min(1, factor)

h_abs *= factor

step_accepted = True
###################################
## QR algorithm
## 1. local lyapunov exponents
Mn = np.reshape(y_new[ndim:],(ndim,ndim)).transpose()
Q,R = np.linalg.qr(Mn)
lyap.append(np.log(abs(R.diagonal())))
## 2. change solution for the perturbation vectors
y_new[ndim:] = Q.transpose().reshape(-1)
###################################

else:
h_abs *= max(rk.MIN_FACTOR,
rk.SAFETY * error_norm ** self.error_exponent)
step_rejected = True

self.h_previous = h
self.y_old = y

self.t = t_new
self.y = y_new

self.h_abs = h_abs
self.f = f_new

return True, None

### Now I switch out Scipy's original function
rk.RungeKutta._step_impl = _step_impl_qr

### Now I define my Lorenz eqs and variational eqs:
# Lorenz system by itself
def lorenz_jac(x,y,z):
dxxdt = -sigma
dxydt = sigma
dxzdt = 0
dyxdt = rho-z
dyydt = -1
dyzdt = -x
dzxdt = y
dzydt = x
dzzdt = -beta

return np.array([[dxxdt,dxydt,dxzdt],
[dyxdt,dyydt,dyzdt],
[dzxdt,dzydt,dzzdt]])

# Stacked solve
def lorenz_system_stacked(t, s):
x,y,z = s[0:3]

pertb_vecs = s[3:].reshape((3,3)).transpose()

sysdot = np.array([sigma * (y - x),
x * (rho - z) - y,
x * y - beta * z])

Mdot = np.matmul(lorenz_jac(x,y,z),pertb_vecs).transpose().reshape(-1)

return np.concatenate((sysdot,Mdot))

### Now I specify the parameters and initial conditions
sigma=10
beta=8/3
rho=28
t = np.arange(0, 50, 0.01)  # Time points
initial_state = np.array([0,1,1.05])               # Initial conditions
pertb0= np.eye(3).reshape(-1,1)
pertb0=pertb0[:,0]
s0 = np.concatenate((initial_state,pertb0))

### Now I call solve_ivp to integrate

sol = scipy.integrate.solve_ivp(lorenz_system_stacked, (t[0], t[-1]), s0, method='RK23')

### After the integration is done, I calculate Lyapunov exponents from the local ones
# Calculate Lyapunov exponent
lyap = np.array(lyap).transpose()
le = []
for i in range(np.shape(lyap)[1]):

if i == 0:
le.append(lyap[:,:i+1].transpose())
else:
le.append(np.sum(lyap[:,:i+1],axis=1)/sol.t[i]))



The cool thing here is that I can use any of SciPy's built-in integrators (and their error controls) to calculate Lyapunov exponents. However, I am not getting the correct Lyapunov exponents with this. The system (Lorenz) equations are solved correctly. I am really scratching my head on this right now as the algorithm is exactly the same as the old MATLAB implementation.

If anyone can give this a fresh look and catch my (I believe) very subtle error, it will be greatly appreciated!

• How do you extract the Lyapunov exponents from the Mn list? In the version you commented out, you got the local increments of the exponents, proportional to the time step. If you cumulatively sum them up and divide by the time (cumulative time steps), you should get the standard LE. Commented Dec 7, 2023 at 11:10
• Yes I am doing that after the integration is done. However, the results I am getting are incorrect. I will edit my original post to reflect this step. Commented Dec 7, 2023 at 21:46
• You could also use np.cumsum and divide by t[:, None]. The result should not change, just the execution should be much faster. // It should not have been necessary to go this deep into the stepper class, just using the stepper class and building a time loop around it should be sufficient. You would of course lose the solve_ivp interface. Commented Dec 7, 2023 at 23:13
• It might be redundant and easy to reconstruct, but please make the question complete an add the data of the observe outcome and the expected outcome, an why their difference is unexpected. Commented Dec 7, 2023 at 23:16
• For this classic setup we know the Lyapunov exponents will be something like [0.8, 0,-14], but I am getting like [2.4, 4.5, -5], which is clearly wrong. The Lorenz system is solved correctly though, as can be verified with a quick plot. Commented Dec 8, 2023 at 2:10

Problem is resolved! I missed a couple of transposes...

The original code is corrected and should be working now:

The final values I got for LEs are:

array([  0.67550971,   0.05611249, -13.99882156])


This uses SciPy's built-in RK23 adaptive solver.

• Still, your current solution is only working for explicit RK methods in solve_ivp. The implicit methods use separate stepping functions. You could rather do the Lyapunov exponent computation as a postprocessing, using a cumulative sum on all the substeps that have been performed. Commented Dec 8, 2023 at 7:05
• Yes that's right. But now the stacked system is error controlled and can step adaptively. Commented Dec 9, 2023 at 1:30