ChaosTools lyapunovs with stiff solvers

Specific Domains Modelling & Simulations
, , ,
1

I’ve come across a problem regarding using lyapunovs() with stiff diffeq solvers. It is very similar to this one:
Problem with tangent_integrator: Chaostools tangent_integrator stiff solvers
tangent_integrator not compatible with stiff solvers · Issue #87 · JuliaDynamics/DynamicalSystemsBase.jl · GitHub
Since lyapunovs() also uses a tangent_integrator, the same problem arises when one wants to calculate the lyapunov exponents with a stiff solver. These solvers seem to be incompatible with the tangent_integrator. Turning autodiff off switches to numerical differentiation, which is super slow. Can anyone help me with a workaround (the link to the workaround in the issue above is no longer available), maybe leaving the variables untyped, so it doesn’t run into problems when converting Float64 into Dual, idk ,or tips to make numerical differentiation faster (Jacobi matrix optimalization etc)?

Minimal Working Example:

function lorenz(du,u,p,t)
    sigma = p[1]
    rho = p[2]
    beta = p[3]

    # eom
    du[1] = sigma*(u[2]-u[1])
    du[2] = u[1]*(rho-u[3])-u[2]
    du[3] = u[1]*u[2]-beta*u[3]

end

function lorenz_jac(J,u,p,t)

    sigma = p[1]
    rho = p[2]
    beta = p[3]

    J[1,1] = -sigma
    J[1,2] = sigma
    J[1,3] = 0
    J[2,1] = rho-u[3]
    J[2,2] = -1
    J[2,3] = -u[1]
    J[3,1] = u[2]
    J[3,2] = u[1]
    J[3,3] = -beta
end

state = rand(3)
param = [16.,45.92,4.]
ds = ContinuousDynamicalSystem(lorenz,state,param,lorenz_jac)
lyapunovs(ds,2000.;Ttr = 2000.,alg = Rosenbrock23())

In this case, it is possible to calculate the lyapunov exponents in a matter of minutes using numerical differentiation (autodiff = false). However in my case, the Jacobian looks something like this:

function njacobi(J,u,p)
    a = -p[1]*p[2]/(p[4]-1)
    c = -p[1]*(1-p[2])
    n = convert(Int32,p[4])
    x = u[1:2:end]
    y = u[2:2:end]
    for i in 1:(2*n)
        for j in 1:(2*n)
            if i<(n+1) && j<(n+1)
                if i!=j
                    J[i,j] = a*x[i]
                else
                    J[i,j] = p[1]*(1 - (1 - p[2])*y[i] - p[2]/(n-1)*skipsum_jacobi(u,i))
                end
            elseif i<(n+1) && j>n
                if (j-i) == n
                    J[i,j] = c*x[i]
                else
                    J[i,j] = 0
                end
            elseif i>n && j<(n+1)
                if (i-j) == n
                    J[i,j] = p[3]
                else 
                    J[i,j] = 0
                end
            else
                if i==j
                    J[i,j] = -p[3]
                else 
                    J[i,j] = 0
                end
            end
        end
    end
end

which is supposed to look like this:

4×4 Array{Float64,2}:
  4.92501  -2.27912  -5.7178   0.0    
 -1.36462   4.99815   0.0     -3.42351
  0.1       0.0      -0.1      0.0    
  0.0       0.1       0.0     -0.1

for a 4D version of my system. It can be treated as 4 2X2 block matrices, and the elements are given depending on which block I’m in, hence the conditionals in my Jacobian. The two bottom blocks are always constant and diagonal. The loops are there because I want to change the dimensions from time to time using a dimension parameter in the parameter container. The problem is, in this case, calculating lyapunovs via numerical differentiation runs for ages, especially when I try to increase the dimensions (50,100+). How can this be done properly?

2

Hi rusandris!
Sorry for the late response. If you already implemented your jacobian function in julia then there is no need to calculate that via auto-differentiation. Your stiff solvers should work.
Unfortunately I can not run your mwe from Issues · JuliaDynamics/DynamicalSystemsBase.jl · GitHub

function lorenz(du,u,p,t)
    sigma = p[1]
    rho = p[2]
    beta = p[3]

    # eom
    du[1] = sigma*(u[2]-u[1])
    du[2] = u[1]*(rho-u[3])-u[2]
    du[3] = u[1]*u[2]-beta*u[3]

end

function lorenz_jac(J,u,p,t)

    sigma = p[1]
    rho = p[2]
    beta = p[3]

    J[1,1] = -sigma
    J[1,2] = sigma
    J[1,3] = 0
    J[2,1] = rho-u[3]
    J[2,2] = -1
    J[2,3] = -u[1]
    J[3,1] = u[2]
    J[3,2] = u[1]
    J[3,3] = -beta
end

state = rand(3)
param = [16.,45.92,4.]
ds = ContinuousDynamicalSystem(lorenz,state,param,lorenz_jac)
lyapunovs(ds,2000.;Ttr = 2000.,alg = Rosenbrock23())

I will try tomorrow. Does this error? And if yes, what error does it give?

EDIT: My bad, the variational equations are the dynamical system and the Jacobian of that is calculated. I will check tomorrow…

3

Hi rusandris!

The problem with tangent_integrator not being able to use dual numbers still persists. Auto-differentiation is not working right now. This probably should be fixed sooner than later. But there is a method on symbolically calculating all Jacobians, using ModelingToolkit. An advantage here is that there is no need to calculate the Jacobian of the original system. See the example code for the Lorentz system below:

import Random
import LinearAlgebra

using DynamicalSystems
using DifferentialEquations
using ModelingToolkit
# define parameters used in our DEs
@parameters t 
@parameters sigma rho beta

# define the time derivative
@derivatives dt'~t

# define the variables we use to define
# the Lorentz right hand side (rhs)
@variables x[1:3](t)

lorenz = [sigma*(x[2]-x[1]),
            x[1]*(rho-x[3])-x[2],
            x[1]*x[2] - beta*x[3]]

# calculate the jacobian of the Lorentz rhs
@variables w[1:3,1:3](t)
lorenz_jac = expand_derivatives.(calculate_jacobian(lorenz, x))

# define the tangent system, i.e. Leorentz system + its jacobian
# This is the rhs we want to evolve
y = vcat(x..., w...)
tangent = vcat(lorenz..., lorenz_jac*w...)

# Together with the dynamics, time-evolution, this
# gives us the Variational Equations
var_eqs = []
for i in 1:length(y)
    push!(var_eqs, dt(y[i]) ~ tangent[i])
end
var_eqs

# now we need to define an ODESystem and an ODEFunction
ode_system = ODESystem(var_eqs)
ode_function = ODEFunction(ode_system, y, [sigma, rho, beta];jac=true)

# define the ODEProblem with parameters
u0 = [19.,20.,50.]     # start values for the lorentz rhs
params = [16.,45.92,4] #[sigma, rho, beta]
# create a random orthogonal matrix
Random.seed!(42)
w0 = Matrix(LinearAlgebra.qr(Random.rand(3,3)).Q) # start values for the jacobian (of the lorentz rhs)

ode_problem = ODEProblem(ode_function, hcat(u0, w0), (0., Inf), params)

# create the integrator
integ = init(ode_problem, Rosenbrock23(); save_everystep = false)

# The system is evolved until system time N * dt.
# Note that dt is not necessarily the steps the
# integrator does.
N = 1000
dt = 1
Ttr = 0 # pre evolving time
lyapunovs(integ, N, dt, Ttr) # defined in ChaosTools

The output:

3-element Array{Float64,1}:
1.481185142306926
0.010637107238963809
-22.61235512981889ockquote

Why is this working? The stiff solvers need the Jacobian of tangent system. This is now calculated symbolically and optimized julia functions are generated.

This is very detailed and very likely can be shorten or written more elegantly.

Best,

jamble

3 Likes
4

Hi jamblejoe!
Thanks! Looks like an awesome workaround! Will check it out tomorrow!