Calculating Lyapunov Exponents with ModelingToolkit

Hi all,
suggested in this thread here Chaostools tangent_integrator stiff solvers, one could use `ModelingToolkit` the approach the following scenario.

One has a dynamical system defined by n differential equations given by the following:

d/dt x = f(x)

Now we want to calculate Lyapunov exponents of this system. As a first step we need to solve the coupled system

d/dt x = f(x)
d/dt W = Df(x)*W,

where W is nxn-Matrix. If one does not want to provide the Jacobian-Matrix Df by hand one could calculate it with `ModelingToolkit` and then set up the coupled system. See this example with Lorenz system:

``````@parameters t sigma rho beta
@variables x(t) y(t) z(t)
@derivatives dt'~t

# The Lorenz system
f(x,y,z) = [sigma*(y-x),
x*(rho-z)-y,
x*y - beta*z]

# creating the tangent dynamics
@variables w11(t) w12(t) w13(t) w21(t) w22(t) w23(t) w31(t) w32(t) w33(t)

w = [w11 w12 w13; w21 w22 w23; w31 w32 w33]

jac(x,y,z) = expand_derivatives.(calculate_jacobian(f(x,y,z), [x, y, z]))
jac_w = jac(x,y,z)*w

eqs = []
push!(eqs, dt(x) ~ f(x,y,z)[1])
push!(eqs, dt(y) ~ f(x,y,z)[2])
push!(eqs, dt(z) ~ f(x,y,z)[3])

for (index, val) in enumerate(w)
push!(eqs, dt(val) ~ jac_w[index])
end
eqs
``````

This works well and results in

12-element Array{Any,1}:
Equation((D’~t())(x(t())), sigma() * (y(t()) - x(t())))
Equation((D’~t())(y(t())), x(t()) * (rho() - z(t())) - y(t()))
Equation((D’~t())(z(t())), x(t()) * y(t()) - beta() * z(t()))
Equation((D’~t())(w11(t())), ((sigma() * -1) * w11(t()) + sigma() * w21(t())) + 0 * w31(t()))
Equation((D’~t())(w21(t())), ((rho() - z(t())) * w11(t()) + -1 * w21(t())) + (x(t()) * -1) * w31(t()))
Equation((D’~t())(w31(t())), (y(t()) * w11(t()) + x(t()) * w21(t())) + (-1 * beta()) * w31(t()))
Equation((D’~t())(w12(t())), ((sigma() * -1) * w12(t()) + sigma() * w22(t())) + 0 * w32(t()))
Equation((D’~t())(w22(t())), ((rho() - z(t())) * w12(t()) + -1 * w22(t())) + (x(t()) * -1) * w32(t()))
Equation((D’~t())(w32(t())), (y(t()) * w12(t()) + x(t()) * w22(t())) + (-1 * beta()) * w32(t()))
Equation((D’~t())(w13(t())), ((sigma() * -1) * w13(t()) + sigma() * w23(t())) + 0 * w33(t()))
Equation((D’~t())(w23(t())), ((rho() - z(t())) * w13(t()) + -1 * w23(t())) + (x(t()) * -1) * w33(t()))
Equation((D’~t())(w33(t())), (y(t()) * w13(t()) + x(t()) * w23(t())) + (-1 * beta()) * w33(t()))

Following the tutorial on the github page of `ModelingToolkit` we now can create an `ODEFunction` using the following code:

``````de = ODESystem(eqs)
ode_func = ODEFunction(de, [x,y,z], [sigma,rho,beta])
``````

The problem now is, that the last line executes for a very long time… I interrupted after 5minutes. The above approach is my first guess on how to use `ModelingToolkit` to create an ODEFunction which consists of the coupled systems - especially the way of creating `w` does not look good and I expect the count of variables to cause the last line not to finish. What can I do here? Maybe @ChrisRackauckas?
Thanks in advance for all replies!

You need to change it to include all of your dependent variables.

Of course! I changed it to

``````ode_func = ODEFunction(de, [x,y,z,w11,w12,w13,w21,w22,w23,w31,w32,w33], [sigma,rho,beta])
``````

But it still takes ages, meaning it does not finish in under 5 minutes, while using 100% processing power of one core.

EDIT: So what I found out by interrupting the function call after different times is that it spends all time in `calculate_factorized_W` and `generate_factorized_W`. I am not 100% sure but I guess this involves inverting a Matrix, in this case a 12x12 Matrix, which can be of course very time-consuming when done symbolically. I do not know, if there can be made further improvements to this functions? @ChrisRackauckas

As a quick fix, can I disable calculating W when creating the `ODEFunction`?

And, where do I find the information, which solver algorithms use the W-matrix?

Do we default to calculating W? That’s not intentional.

I did not look at the implementation details but running the above code, the last line does not finish in sufficient short time and when I interrupt, I get above results…

Remind me after my JuliaCon talk. It’s likely just flipping the kwarg

1 Like

OK, I will. Thanks Chris!

I don’t know what caused the difference but with v0.6.0 it takes less than a 0.01 seconds for building that function, so whatever it was is fixed and I setup the tagging process.

1 Like

You mean the `ModelingToolkit` v.0.6.0? I only see v0.5.0 on github.

Yes, the tag system is taking awhile to merge.

Do you know roughly when v0.6.0 will be available via the `update` function in julia?

EDIT: Does finish under a second now for me as well with version 0.6.0.

I just continue in this post here:

For creating the tangent dynamics I have to set up dummy variables with

``````@variables w11(t) w12(t) w13(t) w21(t) w22(t) w23(t) w31(t) w32(t) w33(t)

w = [w11 w12 w13; w21 w22 w23; w31 w32 w33]
``````

See e.g. the initial post. This can be very tedious if we want to move to bigger systems. What would be the right approach to handle arbitrarily sized systems? Maybe @ChrisRackauckas ?

`w[1:3,1:3]`?

1 Like

Hi @ChrisRackauckas!
This works fine until I want to solve the ode problem. See the following code:

``````@parameters t sigma rho beta
@variables x(t) y(t) z(t)
@derivatives dt'~t

f(x,y,z) = [sigma*(y-x),
x*(rho-z)-y,
x*y - beta*z]

@variables w[1:3,1:3](t)
jac(x,y,z) = expand_derivatives.(calculate_jacobian(f(x,y,z), [x, y, z]))
jac_w = jac(x,y,z)*w

eqs = []
push!(eqs, dt(x) ~ f(x,y,z)[1])
push!(eqs, dt(y) ~ f(x,y,z)[2])
push!(eqs, dt(z) ~ f(x,y,z)[3])

jac_w = jac(x,y,z)*w

for (index, val) in enumerate(w)
push!(eqs, dt(val) ~ jac_w[index])
end

de = ODESystem(eqs)
ode_func = ODEFunction(de, [x,y,z,w[1:3,1:3]...], [sigma,rho,beta])

# create the tangent ode problem
u0 = [19.,20.,50.]
Random.seed!(42)
w0 = Matrix(LinearAlgebra.qr(Random.rand(3,3)).Q)
ode_prob = ODEProblem(ode_func, hcat(u0, w0), (0., 10.))
``````

Everything works well unitl here. Then invoking

``````solve(ode_prob)
``````

results in

BoundsError

Stacktrace:
[1] getindex at ./number.jl:78 [inlined]
[2] ##363(::Array{Float64,2}, ::Array{Float64,2}, ::Int64, ::Float64) at /home/goran/.julia/packages/ModelingToolkit/3e5pc/src/utils.jl:55
[3] macro expansion at /home/goran/.julia/packages/ModelingToolkit/3e5pc/src/utils.jl:102 [inlined]
[4] macro expansion at ./none:0 [inlined]
[5] fast_invokelatest at ./none:0 [inlined]
[6] out_f_safe at /home/goran/.julia/packages/ModelingToolkit/3e5pc/src/systems/diffeqs/diffeqsystem.jl:218 [inlined]
[7] ODEFunction at /home/goran/.julia/packages/DiffEqBase/aPwRz/src/diffeqfunction.jl:193 [inlined]
[10] #__init#346(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Float64, ::Float64, ::Bool, ::Bool, ::Rational{Int64}, ::Nothing, ::Nothing, ::Rational{Int64}, ::Int64, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(LinearAlgebra.opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:default_set,),Tuple{Bool}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}, ::Array{Array{Float64,2},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/goran/.julia/packages/OrdinaryDiffEq/BhP0W/src/solve.jl:338
[11] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:default_set,),Tuple{Bool}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}, ::Array{Array{Float64,2},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0 (repeats 5 times)
[12] #__solve#345(::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:default_set,),Tuple{Bool}}}, ::Function, ::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}) at /home/goran/.julia/packages/OrdinaryDiffEq/BhP0W/src/solve.jl:4
[13] (::getfield(DiffEqBase, Symbol("#kw##__solve")))(::NamedTuple{(:default_set,),Tuple{Bool}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}) at ./none:0
[14] #__solve#2(::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::Nothing) at /home/goran/.julia/packages/DifferentialEquations/4jSpr/src/default_solve.jl:15
[15] #__solve at ./none:0 [inlined]
[16] #__solve#1 at /home/goran/.julia/packages/DifferentialEquations/4jSpr/src/default_solve.jl:5 [inlined]
[17] __solve at /home/goran/.julia/packages/DifferentialEquations/4jSpr/src/default_solve.jl:2 [inlined]
[18] #solve#381 at /home/goran/.julia/packages/DiffEqBase/aPwRz/src/solve.jl:41 [inlined]
[19] solve(::ODEProblem{Array{Float64,2},Tuple{Float64,Float64},true,Nothing,ODEFunction{true,getfield(ModelingToolkit, Symbol("#out_f_safe#62")){getfield(ModelingToolkit, Symbol("###363"))},LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}) at /home/goran/.julia/packages/DiffEqBase/aPwRz/src/solve.jl:27
[20] top-level scope at In[38]:1

File an issue. It’ll get lost here.

2 Likes

ok

For those following, parameters weren’t passed:

``````using ModelingToolkit, OrdinaryDiffEq
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
D(y) ~ x*(ρ-z)-y,
D(z) ~ x*y - β*z]

de = ODESystem(eqs)
f = ODEFunction(de, [x,y,z], [σ,ρ,β])
u0 = [19.,20.,50.]
ode_prob = ODEProblem(f, u0, (0., 10.), (10.0,28.0,2.66))
sol = solve(ode_prob, Tsit5())
``````
1 Like