Using ModelingToolkit.jl to model a matrix ODE

Specific Domains Modelling & Simulations
, ,
1

Hi,

I’m trying to create a matrix ODE using ModelingToolkit,jl (to simplify autodiff and code optimization), but I’m receiving an error which I did not find in any of the issues of the library. Any idea what I’m doing wrong?

MWE:

using ModelingToolkit
using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)
@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)
@named sys = ODESystem([D(y) ~ g * y + f * mat * y])

prob = ODEProblem(sys, [y => [1.0 0.0; 0.0 0.0]], [g => 1.0, f => 2.0])

solve(prob)

Stacktrace:

Stacktrace:
  [1] copyto!(dest::Vector{Any}, src::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base ./abstractarray.jl:890
  [2] _collect(cont::UnitRange{Int64}, itr::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:655
  [3] collect(itr::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base ./array.jl:649
  [4] broadcastable(x::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing})
    @ Base.Broadcast ./broadcast.jl:704
  [5] broadcasted
    @ ./broadcast.jl:1300 [inlined]
  [6] broadcast(::typeof(*), ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, ::Symbolics.ArrayOp{AbstractMatrix{Real}})
    @ Base.Broadcast ./broadcast.jl:798
  [7] materialize(bc::Base.Broadcast.Broadcasted{Symbolics.SymWrapBroadcast, Nothing, typeof(*), Tuple{Num, Symbolics.Arr{Num, 2}}})
    @ Symbolics ~/.julia/packages/Symbolics/H8dtg/src/array-lib.jl:150
  [8] broadcast_preserving_zero_d
    @ ./broadcast.jl:849 [inlined]
  [9] *(A::Num, B::Symbolics.Arr{Num, 2})
    @ Base ./arraymath.jl:52
 [10] *(::Num, ::Num, ::Symbolics.Arr{Num, 2})
    @ Base ./operators.jl:655
 [11] top-level scope
    @ ~/Code/julia/control/src/mwe.jl:10
 [12] include(fname::String)
    @ Base.MainInclude ./client.jl:451
 [13] top-level scope
    @ REPL[3]:1
in expression starting at ./mwe.jl:10
2

I asked a similar question in this issue "Axes of broadcast not known" with differential matrix equation · Issue #1740 · SciML/ModelingToolkit.jl · GitHub

If your remove the [ ] around the equation and use D.() you can create the ODESystem.

using ModelingToolkit

using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)

@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)

@named sys = ODESystem(D.(y) ~ (g * y + f * mat * y))

prob = ODEProblem(sys, [y => [1.0 0.0; 0.0 0.0]], [g => 1.0, f => 2.0])

But then your ODEProblem fails and I don’t know how to fix it.

ERROR: ArgumentError: SymbolicUtils.Symbolic[y[1, 1](x), y[2, 1](x), y[1, 2](x), y[2, 2](x)] are missing from the variable map.
3

You also need to specify the timespan of your ODEProblem

This works, but its not exactly what you want. (I hard coded the values of parameters g and f in the equation)

using ModelingToolkit

using DifferentialEquations: solve

@variables x y[1:2, 1:2](x)

@parameters f g

mat = [1.0 2.0; 0.0 0.5]

D = Differential(x)

@named sys = ODESystem(D.(y) ~ (1.0 * y + 2.0 * mat * y))

t0 = 0.0

tf = 1.0

prob = ODEProblem(sys,[1.0 0.0; 0.0 0.0],(t0,tf))

solve(prob)
4 
eqs = D.(y) ~ (g * y + f * mat * y)
@named sys = ODESystem(eqs)

y0 = collect(y .=> [1.0 0.0; 0.0 0.0])
p = [g => 1.0, f => 2.0]
tspan = (0.0, 0.5)

prob = ODEProblem(sys, y0, tspan, p)
3 Likes