ComponentArrays

https://github.com/SciML/ComponentArrays.jl are array blocks that can be accessed through a named index.

You can compose differential equations without a DSL like ModelingToolkit.

using ComponentArrays
using OrdinaryDiffEq
using SimpleUnPack: @unpack

Composing Lorenz and Lotka-Volterra systems¶

Lorenz system¶

function lorenz!(D, u, p, t; f=0.0)
    @unpack σ, ρ, β = p
    @unpack x, y, z = u

    D.x = σ*(y - x)
    D.y = x*(ρ - z) - y - f
    D.z = x*y - β*z
    return nothing
end

lorenz! (generic function with 1 method)

tspan = (0.0, 20.0)
lorenz_p = (σ=10.0, ρ=28.0, β=8/3)
lorenz_ic = ComponentArray(x=1.0, y=0.0, z=0.0)
lorenz_prob = ODEProblem(lorenz!, lorenz_ic, tspan, lorenz_p)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(x = 1, y = 2, z = 3)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 20.0)
u0: ComponentVector{Float64}(x = 1.0, y = 0.0, z = 0.0)

Lotka-Volterra system¶

function lotka!(D, u, p, t; f=0.0)
    @unpack α, β, γ, δ = p
    @unpack x, y = u

    D.x =  α*x - β*x*y + f
    D.y = -γ*y + δ*x*y
    return nothing
end

lotka! (generic function with 1 method)

lotka_p = (α=2/3, β=4/3, γ=1.0, δ=1.0)
lotka_ic = ComponentArray(x=1.0, y=1.0)
lotka_prob = ODEProblem(lotka!, lotka_ic, tspan, lotka_p)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(x = 1, y = 2)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 20.0)
u0: ComponentVector{Float64}(x = 1.0, y = 1.0)

Composed system¶

function composed!(D, u, p, t)
    c = p.c
    @unpack lorenz, lotka = u

    lorenz!(D.lorenz, lorenz, p.lorenz, t, f=c*lotka.x)
    lotka!(D.lotka, lotka, p.lotka, t, f=c*lorenz.x)
    return nothing
end

composed! (generic function with 1 method)

comp_p = (lorenz=lorenz_p, lotka=lotka_p, c=0.01)
comp_ic = ComponentArray(lorenz=lorenz_ic, lotka=lotka_ic)
comp_prob = ODEProblem(composed!, comp_ic, tspan, comp_p)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(lorenz = ViewAxis(1:3, Axis(x = 1, y = 2, z = 3)), lotka = ViewAxis(4:5, Axis(x = 1, y = 2)))}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 20.0)
u0: ComponentVector{Float64}(lorenz = (x = 1.0, y = 0.0, z = 0.0), lotka = (x = 1.0, y = 1.0))

Solve problem¶

We can solve the composed system...

@time comp_sol = solve(comp_prob);

 11.345849 seconds (20.22 M allocations: 1.307 GiB, 3.09% gc time, 99.99% compilation time)

Extract a certain state variable.

map(comp_sol(0:0.1:20).u) do u
    u.lotka.x
end

201-element Vector{Float64}:
 1.0
 0.9365676092645384
 0.8792888960773919
 0.8317135260245978
 0.798159657171654
 0.7657768717743701
 0.7244406470706198
 0.6818769363821192
 0.6418734059212958
 0.6050209536516953
 ⋮
 3.4580438620363654
 3.666085134693197
 3.875947221099235
 4.0803447007279
 4.268914746962764
 4.434969037310843
 4.56792255836667
 4.638349109679531
 4.607328225450531

or we can unit test one of the component systems

lotka_sol = solve(lotka_prob)
lotka_sol(1.0).x

0.5742850840860039

Symbolic mapping in `ODEFunction`¶

using syms=....
Only valid in non-nested ComponentArrays.

lotka_func = ODEFunction(lotka!; syms=keys(lotka_ic))

SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, typeof(Main.var"##227".lotka!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Dict{Symbol, Int64}, Nothing, Nothing, Nothing, Dict{Any, Any}}, Nothing, Nothing}(Main.var"##227".lotka!, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, SymbolicIndexingInterface.SymbolCache{Dict{Symbol, Int64}, Nothing, Nothing, Nothing, Dict{Any, Any}}(Dict(:y => 2, :x => 1), nothing, nothing, nothing, Dict{Any, Any}()), nothing, nothing)

@time lab_sol = solve(ODEProblem(lotka_func, lotka_ic, tspan, lotka_p));

  6.688113 seconds (9.91 M allocations: 652.088 MiB, 1.87% gc time, 99.99% compilation time)

lab_sol(0.0:0.1:20, idxs=:x)

t: 0.0:0.1:20.0
u: 201-element Vector{Float64}:
 1.0
 0.9356430867230899
 0.8761715473663337
 0.8218112641466815
 0.7726016170719231
 0.7284340850413794
 0.6890986518315521
 0.6543285164232394
 0.6238163035211391
 0.597240645280455
 ⋮
 0.6035774598802119
 0.6219652417879952
 0.6417886677925122
 0.6630691711464052
 0.6858289176558621
 0.7100908056806172
 0.7358784661339564
 0.763216262482708
 0.792129290747252

This notebook was generated using Literate.jl.