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

Question

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 &gt; max_step:
        h_abs = max_step
    elif self.h_abs &lt; 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 &lt; 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) &gt; 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 &lt; 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. — Lutz Lehmann, Dec 07 '23 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. — Axel Wang, Dec 07 '23 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. — Lutz Lehmann, Dec 07 '23 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. — Lutz Lehmann, Dec 07 '23 at 23:16
@LutzLehmann I think you are referring to this for the time loop method? I think you can only do constant time steps (which is fine for Lorenz) or construct your own time step size controls. My main reason of modifying SciPy is so that I can take advantage of the step size control already built in there and look at some relaxation oscillation problems that can only be solved by adaptive time methods. — Axel Wang, Dec 08 '23 at 02:09
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. — Axel Wang, Dec 08 '23 at 02:10

score 2 · Accepted Answer · answered Dec 08 '23 at 06:15

2

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.

answered Dec 08 '23 at 06:15

Axel Wang

197
7

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. – Laurent90 Dec 08 '23 at 07:05
Yes that's right. But now the stacked system is error controlled and can step adaptively. – Axel Wang Dec 09 '23 at 01:30

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

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

Now I switch out Scipy's original function

Now I define my Lorenz eqs and variational eqs:

Lorenz system by itself

Stacked solve

Now I specify the parameters and initial conditions

Now I call solve_ivp to integrate

After the integration is done, I calculate Lyapunov exponents from the local ones

Calculate Lyapunov exponent

1 Answers1