2
$\begingroup$

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!

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Dec 7, 2023 at 11:10
  • $\begingroup$ 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. $\endgroup$
    – Axel Wang
    Dec 7, 2023 at 21:46
  • $\begingroup$ 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. $\endgroup$ Dec 7, 2023 at 23:13
  • $\begingroup$ 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. $\endgroup$ Dec 7, 2023 at 23:16
  • 1
    $\begingroup$ 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. $\endgroup$
    – Axel Wang
    Dec 8, 2023 at 2:10

1 Answer 1

2
$\begingroup$

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

The original code is corrected and should be working now:

enter image description here enter image description here

The final values I got for LEs are:

array([  0.67550971,   0.05611249, -13.99882156])

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

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Laurent90
    Dec 8, 2023 at 7:05
  • $\begingroup$ Yes that's right. But now the stacked system is error controlled and can step adaptively. $\endgroup$
    – Axel Wang
    Dec 9, 2023 at 1:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.