Simple Demonstration of a Symbolic System Structure

In this tutorial we will show how to implement a system structure type for defining the symbolic indexing of a domain-specific language. This tutorial will show how the SymbolCache type is defined to take in arrays of symbols for its independent, dependent, and parameter variable names and uses that to define the symbolic indexing interface.

Defining the ODE

For this example, we will use the Robertson equations:

\[\begin{aligned} \frac{dy_1}{dt} &= -0.04y₁ + 10^4 y_2 y_3 \\ \frac{dy_2}{dt} &= 0.04 y_1 - 10^4 y_2 y_3 - 3*10^7 y_{2}^2 \\ \frac{dy_3}{dt} &= 3*10^7 y_{2}^2 \\ \end{aligned}\]

The in-place function for this ODE system can be defined as:

function rober!(du, u, p, t)
    y₁, y₂, y₃ = u
    k₁, k₂, k₃ = p
    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
    du[2] = k₁ * y₁ - k₂ * y₂^2 - k₃ * y₂ * y₃
    du[3] = k₂ * y₂^2
    nothing
end

rober! (generic function with 1 method)

To add symbolic names for the states in this example, a SymbolCache can be created and passed as the sys keyword argument to the ODEFunction constructor, as shown below:

using OrdinaryDiffEq, SymbolicIndexingInterface

sys = SymbolCache([:y₁, :y₂, :y₃])
odefun = ODEFunction(rober!; sys = sys)

This is then used to create and solve the ODEProblem

prob = ODEProblem(odefun, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4])
sol = solve(prob, Rosenbrock23())

retcode: Success
Interpolation: specialized 2nd order "free" stiffness-aware interpolation
t: 61-element Vector{Float64}:
      0.0
      3.196206628740808e-5
      0.00014400709669791316
      0.00025605212710841824
      0.00048593872520561496
      0.0007179482298405134
      0.0010819240431281
      0.0014801655696439716
      0.0020679567767069723
      0.0028435846274559506
      ⋮
  25371.934242574636
  30784.11997687391
  37217.42637183451
  44850.61378119757
  53893.69155452613
  64593.73799781129
  77241.71960460565
  92180.81944906656
 100000.0
u: 61-element Vector{Vector{Float64}}:
 [1.0, 0.0, 0.0]
 [0.9999987215181657, 1.2780900152625978e-6, 3.9181897521319503e-10]
 [0.9999942397327672, 5.718510591908378e-6, 4.175664083814196e-8]
 [0.9999897579685716, 9.992106860321925e-6, 2.499245680985423e-7]
 [0.99998056266788, 1.7833624284594367e-5, 1.6037078354141817e-6]
 [0.9999712826600025, 2.403488607694387e-5, 4.682453920643883e-6]
 [0.9999567250106862, 3.0390689558961382e-5, 1.2884299754832173e-5]
 [0.999940798607161, 3.388427373871301e-5, 2.531711910029394e-5]
 [0.9999172960308617, 3.583508668595173e-5, 4.686888245237194e-5]
 [0.9998862913719596, 3.641240165367128e-5, 7.729622638673522e-5]
 ⋮
 [0.05563507731000894, 2.3546320422346656e-7, 0.9443646872267863]
 [0.047925348764047165, 2.0121495043224827e-7, 0.952074450021001]
 [0.041233421490033735, 1.7192788166891112e-7, 0.9587664065820826]
 [0.035436998370307615, 1.4688361975260182e-7, 0.9645628547460712]
 [0.030425371816164216, 1.2546809318472803e-7, 0.9695745027157406]
 [0.026099132201150312, 1.0715608431718074e-7, 0.9739007606427644]
 [0.02236969169455032, 9.149844876161202e-8, 0.9776302168069995]
 [0.019158563313533397, 7.81109638013881e-8, 0.9808413585755011]
 [0.017827893850751935, 7.258919982166416e-8, 0.9821720335600467]

The solution can now be indexed symbolically:

sol[:y₁]

61-element Vector{Float64}:
 1.0
 0.9999987215181657
 0.9999942397327672
 0.9999897579685716
 0.99998056266788
 0.9999712826600025
 0.9999567250106862
 0.999940798607161
 0.9999172960308617
 0.9998862913719596
 ⋮
 0.05563507731000894
 0.047925348764047165
 0.041233421490033735
 0.035436998370307615
 0.030425371816164216
 0.026099132201150312
 0.02236969169455032
 0.019158563313533397
 0.017827893850751935

sol(1e3, idxs=:y₁)

0.3367546834697414

However, we did not give names to the parameters or the independent variables. They can be specified using SymbolCache as well:

sys = SymbolCache([:y₁, :y₂, :y₃], [:k₁, :k₂, :k₃], :t)
odefun = ODEFunction(rober!; sys = sys)
prob = ODEProblem(odefun, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4])
sol = solve(prob, Rosenbrock23())
getk1 = getp(sys, :k₁)

getk1(prob)

0.04

getk1(sol)

0.04

sol[:t]

61-element Vector{Float64}:
      0.0
      3.196206628740808e-5
      0.00014400709669791316
      0.00025605212710841824
      0.00048593872520561496
      0.0007179482298405134
      0.0010819240431281
      0.0014801655696439716
      0.0020679567767069723
      0.0028435846274559506
      ⋮
  25371.934242574636
  30784.11997687391
  37217.42637183451
  44850.61378119757
  53893.69155452613
  64593.73799781129
  77241.71960460565
  92180.81944906656
 100000.0