Solving system of ode with Modeling Toolkit

Specific Domains Modelling & Simulations
1

I solved this ode system with DifferentialEquations and got the correct result.
I am trying with ModelingToolkit but I don’t get the correct result. What am I doing wrong?

p₀, q₀ = 0.0, 3.0

@variables t
@parameters P(t) Q(t)
∂t = Differential(t)

eqs = [∂t(Q) ~ 0.2P
       ∂t(P) ~ -80.0sin(Q)]
@named sys = ODESystem(eqs)

u₀ = [P=>p₀, Q=>q₀]

t_span = (0.0, 20.0)
prob = ODEProblem(sys, u₀, t_span)
sol = solve(prob, Tsit5())

The solution with DifferentialEquations that gives me the correct result is:

p₀, q₀ = 0.0, 3.0
function mySystem(du,u,p,t)
    du[1] = 0.2u[2]
    du[2] = -80.0sin(u[1])
end

u0 = [q₀, p₀]
tspan = (0.0,20.0)
prob = ODEProblem(mySystem, u0, tspan)
sol = solve(prob, Tsit5())
2

Interesting… the fix is:

using ModelingToolkit, OrdinaryDiffEq
p₀, q₀ = 0.0, 3.0

@variables t
@variables P(t) Q(t)
∂t = Differential(t)

eqs = [∂t(Q) ~ 0.2P
       ∂t(P) ~ -80.0sin(Q)]
@named sys = ODESystem(eqs)

u₀ = [P=>p₀, Q=>q₀]

t_span = (0.0, 20.0)
prob = ODEProblem(sys, u₀, t_span)
sol = solve(prob, Tsit5())

The only change was to make @variables P(t) Q(t) instead of parameters. It seems that by making them parameters, the automatic system detection put them as both variables and parameters, printing out:

julia> @named sys = ODESystem(eqs)
Model sys with 2 equations
States (2):
  Q(t)
  P(t)
Parameters (2):
  Q(t)
  P(t)

So it would be great if you could open an issue on ModelingToolkit.jl. We should add a sanity check that the intersection between states and parameters should be null. Somehow this one slipped by. Thanks for pointing this out!

3 Likes
3

Sorry, my mistake. Declaring P and Q as parameters is wrong. P and Q must be declared as variables. There are no problems in Modeling Toolkit. Thank you.

4

There is a problem with MTK. It should’ve thrown an error here to help you find this. We’ll add that.

1 Like