How to specify nonhomogeneous system with ModelingToolkit

Specific Domains Modelling & Simulations
1

DifferentialEquations.jl can deal with nonhomogeneous systems specified as an ODEProblem, where some parameter also depends on time.

How do I specify the example as an ODESystem using ModelingToolkit?

So far I tried:

using ModelingToolkit, DifferentialEquations
modelingtoolkitize(prob)

using ModelingToolkit, DifferentialEquations
@parameters t
@variables θ(t) ω(t)
D = Differential(t)
@parameters l m g M(t)
eqs = [
    D(θ) ~ ω,
    #D(ω) ~ -3g/(2l)*sin(θ) + 3/(m*l^2)*M(t),
    D(ω) ~ -3g/(2l)*sin(θ) + 3/(m*l^2)*M,
    ]
p = [
    l => 1.0                             # length [m]
    m => 1.0                             # mass[m]
    g => 9.81                            # gravitational acceleration [m/s²]
    M => (t->0.1sin(t))                    # external torque [Nm]
        ]     
@named sys = ODESystem(eqs)    
sys = structural_simplify(sys)
states(sys)
equations(sys) # looks good
u0 = [
    θ => 0.01,                           # initial angular deflection [rad]
    ω => 0.0    
    ]
tspan = (0.0,10.0)    
prob = ODEProblem(sys,u0, tspan, p)
#prob = ODEProblem(sesam1,u0, tspan, p, jac=true)
sol = solve(prob)

but get error (Julia 1.7.1):
LoadError: MethodError: no method matching *(::Int64, ::var"#40#41")

2

Just make M a state and it will simplify out.

3

thanks, this worked.
This way one can use t at the left hand side:

using ModelingToolkit, DifferentialEquations
@parameters t
@variables θ(t) ω(t) M(t)
D = Differential(t)
@parameters l m g 
Mf(t) = 0.1sin(t)
eqs = [
    D(θ) ~ ω,
    #D(ω) ~ -3g/(2l)*sin(θ) + 3/(m*l^2)*M(t),
    D(ω) ~ -3g/(2l)*sin(θ) + 3/(m*l^2)*M,
    M ~ Mf(t)
    ]
p = [
    l => 1.0                             # length [m]
    m => 1.0                             # mass[m]
    g => 9.81                            # gravitational acceleration [m/s²]
        ]     
@named sys = ODESystem(eqs)    
sys = structural_simplify(sys)
states(sys)
equations(sys)
u0 = [
    θ => 0.01,                           # initial angular deflection [rad]
    ω => 0.0,    
    ]
tspan = (0.0,10.0)    
prob = ODEProblem(sys,u0, tspan, p)
sol = solve(prob)
2 Likes