ComponentArrays

ComponentArrays#

SciML/ComponentArrays.jl are array blocks that can be accessed through a named index. You can compose differential equations without ModelingToolkit.

using ComponentArrays
using OrdinaryDiffEq
using SimpleUnPack: @unpack

tspan = (0.0, 20.0)

# 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_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)


# 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_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)

# Composed Lorenz and Lotka-Volterra system
function composed!(D, u, p, t)
    c = p.c #coupling parameter
    @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

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…

comp_sol = solve(comp_prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 241-element Vector{Float64}:
  0.0
  6.181923459793355e-5
  0.0006800115805772691
  0.004299130644888261
  0.010835274422622885
  0.01959555343271236
  0.0316252322166906
  0.047427991793856386
  0.06799287907875844
  0.09424926634485992
  ⋮
 19.470984724131988
 19.56149826542021
 19.638571826913314
 19.71202349285504
 19.776111547247073
 19.842847075547546
 19.911770857107726
 19.98556258790879
 20.0
u: 241-element Vector{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)))}}}}:
 ComponentVector{Float64}(lorenz = (x = 1.0, y = 0.0, z = 0.0), lotka = (x = 1.0, y = 1.0))
 ComponentVector{Float64}(lorenz = (x = 0.9993825332995702, y = 0.0017297323086565809, z = 5.3446465987803325e-8), lotka = (x = 0.9999594060148846, y = 0.9999999987452488))
 ComponentVector{Float64}(lorenz = (x = 0.9932873611408131, y = 0.018962886302111005, z = 6.422601353340659e-6), lotka = (x = 0.999553537391496, y = 0.9999998481912489))
 ComponentVector{Float64}(lorenz = (x = 0.9604308390386475, y = 0.11762827803691758, z = 0.000246932331946337), lotka = (x = 0.9971800835961485, y = 0.9999939361252497))
 ComponentVector{Float64}(lorenz = (x = 0.9125872373872632, y = 0.2874050074697573, z = 0.0014720096775698572), lotka = (x = 0.9929055676177644, y = 0.9999615274469091))
 ComponentVector{Float64}(lorenz = (x = 0.8693644514509589, y = 0.5019139490156515, z = 0.004480547354506518), lotka = (x = 0.987202270832691, y = 0.9998743818678912))
 ComponentVector{Float64}(lorenz = (x = 0.8439920748633932, y = 0.7815147428132007, z = 0.010834300546636073), lotka = (x = 0.9794224855372796, y = 0.9996736178757919))
 ComponentVector{Float64}(lorenz = (x = 0.8618249107449826, y = 1.1413508107615669, z = 0.023033153078366626), lotka = (x = 0.9693021639697955, y = 0.9992685312839769))
 ComponentVector{Float64}(lorenz = (x = 0.9603072729766949, y = 1.6328057689454554, z = 0.04697042251721498), lotka = (x = 0.9563171342036174, y = 0.9985041943402345))
 ComponentVector{Float64}(lorenz = (x = 1.2008321945955154, y = 2.3649414730469003, z = 0.09822934050055844), lotka = (x = 0.9400721183189727, y = 0.997146076484591))
 ⋮
 ComponentVector{Float64}(lorenz = (x = -6.456640310418434, y = -11.301591348479082, z = 13.579301255981097), lotka = (x = 4.048583254557733, y = 0.09036713029021927))
 ComponentVector{Float64}(lorenz = (x = -11.98914582689141, y = -18.65379996284021, z = 22.268644179842802), lotka = (x = 4.238187094020435, y = 0.12011880751447605))
 ComponentVector{Float64}(lorenz = (x = -15.310935211801239, y = -15.5440122666827, z = 36.08580331875077), lotka = (x = 4.388362143336394, y = 0.155077960118357))
 ComponentVector{Float64}(lorenz = (x = -11.968917785545145, y = -3.7661598961030274, z = 38.78443474987127), lotka = (x = 4.51909613045653, y = 0.19987430738288942))
 ComponentVector{Float64}(lorenz = (x = -6.516077535769612, y = 1.2554761404938108, z = 33.15893795214686), lotka = (x = 4.620776349063695, y = 0.2512766104607464))
 ComponentVector{Float64}(lorenz = (x = -2.4716487755263947, y = 1.8751117763255676, z = 27.307129257301032), lotka = (x = 4.707316550624519, y = 0.32092674024836887))
 ComponentVector{Float64}(lorenz = (x = -0.40465156571559013, y = 1.5281547362335937, z = 22.583051410704424), lotka = (x = 4.764760203116338, y = 0.41527390525354513))
 ComponentVector{Float64}(lorenz = (x = 0.5766169462153872, y = 1.520820221743716, z = 18.565160004371325), lotka = (x = 4.774793872051657, y = 0.5486699097009238))
 ComponentVector{Float64}(lorenz = (x = 0.7079514594147894, y = 1.5886271745693228, z = 17.87815297902187), lotka = (x = 4.769017950232028, y = 0.5793790495600354))

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

201-element Vector{Float64}:
 1.0
 0.9365676111577816
 0.8792888862135674
 0.8317131775607912
 0.7981623360172602
 0.7657794013009516
 0.7244491807705086
 0.6818931117860487
 0.6418920515039641
 0.6050346686175818
 ⋮
 3.463870472793916
 3.6773336686361615
 3.8947925428201797
 4.110486391957825
 4.314781179879635
 4.498597362739154
 4.654513152532488
 4.75781342700078
 4.769017950232028

…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 available in non-nested ComponentArrays

lotka_func = ODEFunction(lotka!; syms=keys(lotka_ic))
lab_sol = solve(ODEProblem(lotka_func, lotka_ic, tspan, lotka_p))
lab_sol(1.0, idxs=:x)

0.5742850840860039

This notebook was generated using Literate.jl.