ModelingToolkit.jl

ModelingToolkit.jl is an intermediate representation (IR) of computational graphs for scientific computing problems. Its purpose is to be a common target for modeling DSLs in order to allow for a common platform for model inspection and transformation. It uses a tagged variable IR in order to allow specification of complex models and allow for transformations of models. It has ways to plug into its function registration and derivative system so that way it can interact nicely with user-defined routines. Together, this is an abstract form of a scientific model that is easy for humans to generate but also easy for programs to manipulate.

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation which contains the un-released features.

High Level Examples

First let's define a second order riff on the Lorenz equations, symbolically lower it to a first order system, symbolically generate the Jacobian function for the numerical integrator, and solve it.

using ModelingToolkit, OrdinaryDiffEq

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(D(x)) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

sys = ODESystem(eqs)
sys = ode_order_lowering(sys)

u0 = [D(x) => 2.0,
      x => 1.0,
      y => 0.0,
      z => 0.0]

p  = [σ => 28.0,
      ρ => 10.0,
      β => 8/3]

tspan = (0.0,100.0)
prob = ODEProblem(sys,u0,tspan,p,jac=true)
sol = solve(prob,Tsit5())
using Plots; plot(sol,vars=(x,y))

This automatically will have generated fast Jacobian functions, making it more optimized than directly building a function. In addition, we can then use ModelingToolkit to compose multiple ODE subsystems. Now let's define two interacting Lorenz equations and simulate the resulting Differential-Algebriac Equation (DAE):

using ModelingToolkit, OrdinaryDiffEq

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs,name=:lorenz1)
lorenz2 = ODESystem(eqs,name=:lorenz2)

@variables a
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a*γ]
connected = ODESystem(connections,t,[a],[γ],systems=[lorenz1,lorenz2])

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0,
      a => 2.0]

p  = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3,
      γ => 2.0]

tspan = (0.0,100.0)
prob = ODEProblem(connected,u0,tspan,p)
sol = solve(prob,Rodas5())

using Plots; plot(sol,vars=(a,lorenz1.x,lorenz2.z))