Calculating Lyapunov Exponents with ModelingToolkit

Specific Domains Modelling & Simulations
1

Hi all,
suggested in this thread here Chaostools tangent_integrator stiff solvers - #26 by jamblejoe, 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!

2

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

3

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?

4

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

5

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…

6

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

1 Like
7

OK, I will. Thanks Chris!

8

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
9

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

10

Yes, the tag system is taking awhile to merge.

11

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.

12

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 ? :slight_smile:

13

w[1:3,1:3]?

1 Like
14

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]
[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

15

File an issue. It’ll get lost here.

2 Likes
16

ok :slight_smile:

17

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
18

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!