# Calculating Lyapunov Exponents with ModelingToolkit

**URL:** <https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641>\
**Category:** Modelling & Simulations\
**Created:** [July 22, 2019, 5:05pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641 "2019-07-22T17:05:29Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 22, 2019, 5:05pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/1 "2019-07-22T17:05:29Z")

</div>

Hi all,  
suggested in this thread here [Chaostools tangent\_integrator stiff solvers - #26 by jamblejoe](https://discourse.julialang.org/t/chaostools-tangent-integrator-stiff-solvers/26464/26), 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:

```julia
@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:

```julia
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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2019, 2:52am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/2 "2019-07-23T02:52:03Z")

</div>

> [@jamblejoe](#):
>
> Following the tutorial on the github page of `ModelingToolkit` we now can create an `ODEFunction` using the following code:

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 23, 2019, 11:30am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/3 "2019-07-23T11:30:20Z")

</div>

Of course! I changed it to

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2019, 7:25pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/4 "2019-07-23T19:25:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 23, 2019, 7:28pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/5 "2019-07-23T19:28:10Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2019, 7:30pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/6 "2019-07-23T19:30:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 23, 2019, 7:30pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/7 "2019-07-23T19:30:45Z")

</div>

OK, I will. Thanks Chris!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 28, 2019, 4:06pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/8 "2019-07-28T16:06:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 29, 2019, 9:38am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/9 "2019-07-29T09:38:46Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 29, 2019, 9:57am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/10 "2019-07-29T09:57:20Z")

</div>

Yes, the tag system is taking awhile to merge.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 29, 2019, 1:42pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/11 "2019-07-29T13:42:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 30, 2019, 6:02pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/12 "2019-07-30T18:02:22Z")

</div>

I just continue in this post here:

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

```julia
@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 ? 🙂

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 1, 2019, 9:36pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/13 "2019-08-01T21:36:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [August 9, 2019, 12:11pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/14 "2019-08-09T12:11:49Z")

</div>

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

```julia
@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

```julia
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]  
> [8] initialize!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}},true,Array{Float64,2},Float64,Nothing,Float64,Float64,Float64,Array{Array{Float64,2},1},OrdinaryDiffEq.ODECompositeSolution{Float64,3,Array{Array{Float64,2},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,2},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{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},Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{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},Array{Float64,2},Nothing},DiffEqDiffTools.UJacobianWrapper{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},Float64,Nothing},DefaultLinSolve,DiffEqDiffTools.JacobianCache{Array{Float64,2},Array{Float64,2},Array{Float64,2},UnitRange{Int64},Nothing,Val{:forward},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,2},Array{Float64,2},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}}},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},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{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},Array{Float64,2},Nothing},DiffEqDiffTools.UJacobianWrapper{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},Float64,Nothing},DefaultLinSolve,DiffEqDiffTools.JacobianCache{Array{Float64,2},Array{Float64,2},Array{Float64,2},UnitRange{Int64},Nothing,Val{:forward},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,2},Array{Float64,2},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE\_DEFAULT\_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE\_DEFAULT\_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE\_DEFAULT\_PROG\_MESSAGE),typeof(DiffEqBase.ODE\_DEFAULT\_UNSTABLE\_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,2},Float64,Nothing}, ::OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}) at /home/goran/.julia/packages/OrdinaryDiffEq/BhP0W/src/perform\_step/low\_order\_rk\_perform\_step.jl:623  
> [9] initialize!(::OrdinaryDiffEq.ODEIntegrator{CompositeAlgorithm{Tuple{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}},true,Array{Float64,2},Float64,Nothing,Float64,Float64,Float64,Array{Array{Float64,2},1},OrdinaryDiffEq.ODECompositeSolution{Float64,3,Array{Array{Float64,2},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,2},1},1},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}},OrdinaryDiffEq.CompositeInterpolationData{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},Array{Array{Float64,2},1},Array{Float64,1},Array{Array{Array{Float64,2},1},1},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{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},Array{Float64,2},Nothing},DiffEqDiffTools.UJacobianWrapper{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},Float64,Nothing},DefaultLinSolve,DiffEqDiffTools.JacobianCache{Array{Float64,2},Array{Float64,2},Array{Float64,2},UnitRange{Int64},Nothing,Val{:forward},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,2},Array{Float64,2},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}}},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},OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{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},Array{Float64,2},Nothing},DiffEqDiffTools.UJacobianWrapper{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},Float64,Nothing},DefaultLinSolve,DiffEqDiffTools.JacobianCache{Array{Float64,2},Array{Float64,2},Array{Float64,2},UnitRange{Int64},Nothing,Val{:forward},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,2},Array{Float64,2},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE\_DEFAULT\_NORM),typeof(LinearAlgebra.opnorm),CallbackSet{Tuple{},Tuple{}},typeof(DiffEqBase.ODE\_DEFAULT\_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE\_DEFAULT\_PROG\_MESSAGE),typeof(DiffEqBase.ODE\_DEFAULT\_UNSTABLE\_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1}},Array{Float64,2},Float64,Nothing}, ::OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.Rosenbrock23Cache{Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.Rosenbrock23ConstantCache{Float64,typeof(identity),typeof(identity)},DiffEqDiffTools.TimeGradientWrapper{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},Array{Float64,2},Nothing},DiffEqDiffTools.UJacobianWrapper{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},Float64,Nothing},DefaultLinSolve,DiffEqDiffTools.JacobianCache{Array{Float64,2},Array{Float64,2},Array{Float64,2},UnitRange{Int64},Nothing,Val{:forward},Float64,Val{true}},DiffEqDiffTools.GradientCache{Nothing,Array{Float64,2},Array{Float64,2},Val{:forward},Float64,Val{true}}}},AutoSwitch{Tsit5,Rosenbrock23{0,false,DefaultLinSolve,DataType},Rational{Int64},Float64}}) at /home/goran/.julia/packages/OrdinaryDiffEq/BhP0W/src/perform\_step/composite\_perform\_step.jl:38  
> [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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 9, 2019, 12:59pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/15 "2019-08-09T12:59:55Z")

</div>

File an issue. It’ll get lost here.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [August 9, 2019, 2:37pm UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/16 "2019-08-09T14:37:12Z")

</div>

ok 🙂

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 11, 2019, 10:38am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/17 "2019-08-11T10:38:20Z")

</div>

For those following, parameters weren’t passed:

```julia
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())

```

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [August 14, 2019, 10:49am UTC](https://discourse.julialang.org/t/calculating-lyapunov-exponents-with-modelingtoolkit/26641/18 "2019-08-14T10:49:52Z")

</div>

Sorry, that was my bad!

Your idea to make the error messages of ModelingToolkit a bit more non expert user friendly is great! Thanks @ChrisRackauckas for all your effort and time!