ModelingToolkit: troubles using ode_order_lowering (A differentiated state's operation must be a `Sym`)

Specific Domains Modelling & Simulations
1

I am trying to model a simple mechanical chain composed of a few masses interconnected with springs and exposed to friction. Such model comes in the form of a few coupled second-order differential equations. These I model below.

using ModelingToolkit
using OrdinaryDiffEq

@parameters a b               # Some physical parameters.
@variables t z[1:10](t)       # Positions of individual masses.
D = Differential(t)           # Derivative with respect to time.
n = 5                         # Number of masses.

F = [a*(z[i+1] - 2*z[i] + z[i-1]) for i=2:n-1]
pushfirst!(F,a*(z[2] - z[1]))
push!(F,a*(z[n-1] - z[n]))    # Models stiffness of the springs to the left and right.

eqs = [D(D(z[i])) ~ -b*D(z[i]) + F[i] for i in 1:n]
                              # Acceleration given by the friction and interaction with neighbors.

sys = ODESystem(eqs)

So far so good. Now I want to reduce the set of second-order ODEs to a set of first-order ODEs:

sys = ode_order_lowering(sys)
ERROR: ArgumentError: A differentiated state's operation must be a `Sym`, so states like `D(u + u)` are disallowed. Got `z[1]`.
Stacktrace:
 [1] diff2term(O::Term{Real, Nothing})
   @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/utils.jl:109
 [2] lower_varname(var::Term{Real, Nothing}, idv::Sym{Real, Nothing}, order::Int64)
   @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/utils.jl:160
 [3] ode_order_lowering(eqs::Vector{Equation}, iv::Sym{Real, Nothing}, states::Vector{Term{Real, Nothing}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MYBNE/src/systems/diffeqs/first_order_transform.jl:29
 [4] ode_order_lowering(sys::ODESystem)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MYBNE/src/systems/diffeqs/first_order_transform.jl:9
 [5] top-level scope
   @ none:1

I confess I am not able to decipher the error message and make the model work. Anyone can help?

2

I think this is because of the recent change to MTK that implements arrays of symbolic variables as symbolic arrays (in time it will allow more straightforward vector operations in MTK). The easiest fix is to generate an array of symbolic variables to mimic the old behaviour of MTK. The way I’d do it is

@parameters a b               # Some physical parameters.
@variables t
z = map(first, [@variables $sym(t) for sym in (Symbol(:z, i) for i in 1:10)]) # Positions of individual masses.
D = Differential(t)           # Derivative with respect to time.
n = 5                         # Number of masses.

which uses interpolation to create the desired symbolic variables. (There might be a cleaner way of doing this but I couldn’t work out what it is short of copying code out of the @variables macro, but that seems a bit brittle.) This gives you

julia> z = map(first, [@variables $sym(t) for sym in (Symbol(:z, i) for i in 1:10)])
10-element Vector{Num}:
  z1(t)
  z2(t)
  z3(t)
  z4(t)
  z5(t)
  z6(t)
  z7(t)
  z8(t)
  z9(t)
 z10(t)

If you want the nice subscripts like before (e.g., z₁), you can use the map_subscripts function from Symbolics (available via ModelingToolkit.Symbolics.map_subscripts.

1 Like