ModelingToolkit.jl, An IR and Compiler for Scientific Models

Chris Rackauckas

A lot of people are building modeling languages for their specific domains. However, while the syntax my vary greatly between these domain-specific languages (DSLs), the internals of modeling frameworks are surprisingly similar: building differential equations, calculating Jacobians, etc.

ModelingToolkit.jl is metamodeling systemitized

After building our third modeling interface, we realized that this problem can be better approached by having a reusable internal structure which DSLs can target. This internal is ModelingToolkit.jl: an Intermediate Representation (IR) with a well-defined interface for defining system transformations and compiling to Julia functions for use in numerical libraries. Now a DSL can easily be written by simply defining the translation to ModelingToolkit.jl's primatives and querying for the mathematical quantities one needs.

Basic usage: defining differential equation systems, with performance!

Let's explore the IR itself. ModelingToolkit.jl is friendly to use, and can used as a symbolic DSL in its own right. Let's define and solve the Lorenz differential equation system using ModelingToolkit to generate the functions:

using ModelingToolkit

### Define a differential equation system

@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]
de = ODESystem(eqs)
ode_f = ODEFunction(de, [x,y,z], [σ,ρ,β])

### Use in DifferentialEquations.jl

using OrdinaryDiffEq
u₀ = ones(3)
tspan = (0.0,100.0)
p = [10.0,28.0,10/3]
prob = ODEProblem(ode_f,u₀,tspan,p)
sol = solve(prob,Tsit5())

using Plots
plot(sol,vars=(1,2,3))

ModelingToolkit is a compiler for mathematical systems

At its core, ModelingToolkit is a compiler. It's IR is its type system, and its output are Julia functions (it's a compiler for Julia code to Julia code, written in Julia).

DifferentialEquations.jl wants a function f(du,u,p,t) for defining an ODE system, which is what ModelingToolkit.jl is building.

generate_function(de, [x,y,z], [σ,ρ,β])
:((##383, u, p, t)->begin
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:44 =#
          let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
              ##383[1] = (*)(σ, (-)(y, x))
              ##383[2] = (-)((*)(x, (-)(ρ, z)), y)
              ##383[3] = (-)((*)(x, y), (*)(β, z))
          end
      end)

A special syntax in DifferentialEquations.jl for small static ODE systems uses f(u,p,t), which can be generated as well:

generate_function(de, [x,y,z], [σ,ρ,β]; version=ModelingToolkit.SArrayFunction)
:((u, p, t)->begin
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:48 =#
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:49 =#
          X = let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
                  ((*)(σ, (-)(y, x)), (-)((*)(x, (-)(ρ, z)), y), (-)((*)(x, y), (*)(β, z)))
              end
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:50 =#
          T = StaticArrays.similar_type(typeof(u), eltype(X))
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:51 =#
          T(X)
      end)

ModelingToolkit.jl can be used to calculate the Jacobian of the differential equation system:

jac = calculate_jacobian(de)
3×3 Array{ModelingToolkit.Expression,2}:
     σ() * -1           σ()  Constant(0)
 ρ() - z(t())  Constant(-1)  x(t()) * -1
       y(t())        x(t())     -1 * β()

It will automatically generate functions for using this Jacobian within the stiff ODE solvers for faster solving:

jac_expr = generate_jacobian(de)
:((##384, u, p, t)->begin
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:44 =#
          let (x, y, z, ρ, σ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
              ##384[1] = (*)(σ, -1)
              ##384[2] = (-)(ρ, z)
              ##384[3] = y
              ##384[4] = σ
              ##384[5] = -1
              ##384[6] = x
              ##384[7] = 0
              ##384[8] = (*)(x, -1)
              ##384[9] = (*)(-1, β)
          end
      end)

It can even do fancy linear algebra. Stiff ODE solvers need to perform an LU-factorization which is their most expensive part. But ModelingToolkit.jl can skip this operation and instead generate the analytical solution to a matrix factorization, and build a Julia function for directly computing the factorization, which is then optimized in LLVM compiler passes.

ModelingToolkit.generate_factorized_W(de)[1]
:((##385, u, p, gam, t)->begin
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:44 =#
          let (x, y, z, ρ, σ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
              ##385[1] = (+)((*)(σ, gam), true)
              ##385[2] = (*)(gam, (-)(ρ, z), -1, (inv)((+)((*)(σ, gam), true)))
              ##385[3] = (*)(gam, y, -1, (inv)((+)((*)(σ, gam), true)))
              ##385[4] = (*)(gam, σ, -1)
              ##385[5] = (-)((+)(gam, true), (*)(gam, (-)(ρ, z), gam, σ, (inv)((+)((*)(σ, gam), true))))
              ##385[6] = (*)((-)((*)(gam, x, -1), (*)(gam, y, gam, σ, (inv)((+)((*)(σ, gam), true)))), (inv)((-)((+)(gam, true), (
*)(gam, (-)(ρ, z), gam, σ, (inv)((+)((*)(σ, gam), true))))))
              ##385[7] = 0
              ##385[8] = (-)((*)(x, gam), 0)
              ##385[9] = (-)((-)((+)((*)(β, gam), true), 0), (*)((-)((*)(gam, x, -1), (*)(gam, y, gam, σ, (inv)((+)((*)(σ, gam), t
rue)))), (inv)((-)((+)(gam, true), (*)(gam, (-)(ρ, z), gam, σ, (inv)((+)((*)(σ, gam), true))))), (-)((*)(x, gam), 0)))
          end
      end)

Solving Nonlinear systems

ModelingToolkit.jl is not just for differential equations. It can be used for any mathematical target that is representable by its IR. For example, let's solve a rootfinding problem F(x)=0. What we do is define a nonlinear system and generate a function for use in NLsolve.jl

@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z])
nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
:((##387, u, p)->begin
          #= /home/alex/.julia/packages/ModelingToolkit/S0mks/src/utils.jl:44 =#
          let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
              ##387[1] = (*)(σ, (-)(y, x))
              ##387[2] = (-)((*)(x, (-)(ρ, z)), y)
              ##387[3] = (-)((*)(x, y), (*)(β, z))
          end
      end)

We can then tell ModelingToolkit.jl to compile this function for use in NLsolve.jl, and then numerically solve the rootfinding problem:

nl_f = @eval eval(nlsys_func)
# Make a closure over the parameters for for NLsolve.jl
f2 = (du,u) -> nl_f(du,u,(10.0,26.0,2.33))

using NLsolve
nlsolve(f2,ones(3))
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.0, 1.0, 1.0]
 * Zero: [2.2228e-10, 2.2228e-10, -9.99034e-11]
 * Inf-norm of residuals: 0.000000
 * Iterations: 3
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 4
 * Jacobian Calls (df/dx): 4

Library of transformations on mathematical systems

The reason for using ModelingToolkit is not just for defining performant Julia functions for solving systems, but also for performing mathematical transformations which may be required in order to numerically solve the system. For example, let's solve a third order ODE. The way this is done is by transforming the third order ODE into a first order ODE, and then solving the resulting ODE. This transformation is given by the ode_order_lowering function.

@derivatives D3'''~t
@derivatives D2''~t
@variables u(t), x(t)
eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
       D2(x) ~ D(x) + 2]
de = ODESystem(eqs)
de1 = ode_order_lowering(de)
ModelingToolkit.ODESystem(ModelingToolkit.DiffEq[DiffEq(u_tt, 1, ((2 * u_tt(t()) + u_t(t())) + x_t(t())) + 1), DiffEq(x_t, 1, x_t(
t()) + 2), DiffEq(u_t, 1, u_tt(t())), DiffEq(u, 1, u_t(t())), DiffEq(x, 1, x_t(t()))], t, ModelingToolkit.Variable[u, x, u_tt, u_t
, x_t], ModelingToolkit.Variable[], Base.RefValue{Array{ModelingToolkit.Expression,2}}(Array{Expression}(0,0)))
de1.eqs
5-element Array{ModelingToolkit.DiffEq,1}:
 ModelingToolkit.DiffEq(u_tt, 1, ((2 * u_tt(t()) + u_t(t())) + x_t(t())) + 1)
 ModelingToolkit.DiffEq(x_t, 1, x_t(t()) + 2)                                
 ModelingToolkit.DiffEq(u_t, 1, u_tt(t()))                                   
 ModelingToolkit.DiffEq(u, 1, u_t(t()))                                      
 ModelingToolkit.DiffEq(x, 1, x_t(t()))

This has generated a system of 5 first order ODE systems which can now be used in the ODE solvers.

Linear Algebra... for free?

Let's take a look at how to extend ModelingToolkit.jl in new directions. Let's define a Jacobian just by using the derivative primatives by hand:

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t Dx'~x Dy'~y Dz'~z
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
J = [Dx(eqs[1].rhs) Dy(eqs[1].rhs) Dz(eqs[1].rhs)
 Dx(eqs[2].rhs) Dy(eqs[2].rhs) Dz(eqs[2].rhs)
 Dx(eqs[3].rhs) Dy(eqs[3].rhs) Dz(eqs[3].rhs)]
3×3 Array{ModelingToolkit.Operation,2}:
          (D'~x(t()))(σ() * (y(t()) - x(t())))  …           (D'~z(t()))(σ() * (y(t()) - x(t())))
 (D'~x(t()))(x(t()) * (ρ() - z(t())) - y(t()))     (D'~z(t()))(x(t()) * (ρ() - z(t())) - y(t()))
   (D'~x(t()))(x(t()) * y(t()) - β() * z(t()))       (D'~z(t()))(x(t()) * y(t()) - β() * z(t()))

Notice that this writes the derivatives in a "lazy" manner. If we want to actually compute the derivatives, we can expand out those expressions:

J = expand_derivatives.(J)
3×3 Array{ModelingToolkit.Expression,2}:
     σ() * -1           σ()  Constant(0)
 ρ() - z(t())  Constant(-1)  x(t()) * -1
       y(t())        x(t())     -1 * β()

Here's the magic of ModelingToolkit.jl: Julia treats ModelingToolkit expressions like a Number, and so generic numerical functions are directly usable on ModelingToolkit expressions! Let's compute the LU-factorization of this Jacobian we defined using Julia's Base linear algebra library.

using LinearAlgebra
luJ = lu(J)
LinearAlgebra.LU{ModelingToolkit.Expression,Array{ModelingToolkit.Expression,2}}
L factor:
3×3 Array{ModelingToolkit.Expression,2}:
                    Constant(1)  …  Constant(0)
 (ρ() - z(t())) * inv(σ() * -1)     identity(0)
         y(t()) * inv(σ() * -1)     Constant(1)
U factor:
3×3 Array{ModelingToolkit.Expression,2}:
    σ() * -1  …                                                                                                                   
                                                                     Constant(0)
 identity(0)                                                                                                                      
                              x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * 0
 identity(0)     (-1 * β() - (y(t()) * inv(σ() * -1)) * 0) - ((x(t()) - (y(t()) * inv(σ() * -1)) * σ()) * inv(-1 - ((ρ() - z(t()))
 * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * 0)
luJ.L
3×3 Array{ModelingToolkit.Expression,2}:
                    Constant(1)  …  Constant(0)
 (ρ() - z(t())) * inv(σ() * -1)     identity(0)
         y(t()) * inv(σ() * -1)     Constant(1)

and the inverse?

invJ = inv(J)
3×3 Array{ModelingToolkit.Operation,2}:
 (σ() * -1) \ ((identity(true) - identity(0) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t()) * inv(σ(
) * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * 
identity(0))) \ ((identity(0) - (y(t()) * inv(σ() * -1)) * identity(true)) - ((x(t()) - (y(t()) * inv(σ() * -1)) * σ()) * inv(iden
tity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(true))))) - σ() *
 ((identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ()) \ ((identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(true)) - (
x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0)) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) -
 (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())) * 
inv(σ() * -1)) * identity(0))) \ ((identity(0) - (y(t()) * inv(σ() * -1)) * identity(true)) - ((x(t()) - (y(t()) * inv(σ() * -1)) 
* σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(
true)))))))  …  (σ() * -1) \ ((identity(0) - identity(0) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t
()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())) * inv(σ
() * -1)) * identity(0))) \ ((identity(true) - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t()) * inv(σ() * -1)) * σ()
) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0))))
) - σ() * ((identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ()) \ ((identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0
)) - (x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0)) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(
t()) - (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t(
))) * inv(σ() * -1)) * identity(0))) \ ((identity(true) - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t()) * inv(σ() *
 -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * ide
ntity(0)))))))
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
    (identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ()) \ ((identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity(true)) -
 (x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0)) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t())
 - (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())) 
* inv(σ() * -1)) * identity(0))) \ ((identity(0) - (y(t()) * inv(σ() * -1)) * identity(true)) - ((x(t()) - (y(t()) * inv(σ() * -1)
) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identit
y(true)))))                                                                                                                       
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
             (identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ()) \ ((identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * identity
(0)) - (x(t()) * -1 - ((ρ() - z(t())) * inv(σ() * -1)) * identity(0)) * (((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((
x(t()) - (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(
t())) * inv(σ() * -1)) * identity(0))) \ ((identity(true) - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t()) * inv(σ()
 * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * i
dentity(0)))))
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                     ((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t(
)) - (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - z(t())
) * inv(σ() * -1)) * identity(0))) \ ((identity(0) - (y(t()) * inv(σ() * -1)) * identity(true)) - ((x(t()) - (y(t()) * inv(σ() * -
1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) * ident
ity(true)))                                                                                                                       
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                                                                                  
                                                                           ((-1 * β() - (y(t()) * inv(σ() * -1)) * identity(0)) - 
((x(t()) - (y(t()) * inv(σ() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (x(t()) * -1 - ((ρ() - 
z(t())) * inv(σ() * -1)) * identity(0))) \ ((identity(true) - (y(t()) * inv(σ() * -1)) * identity(0)) - ((x(t()) - (y(t()) * inv(σ
() * -1)) * σ()) * inv(identity(-1) - ((ρ() - z(t())) * inv(σ() * -1)) * σ())) * (identity(0) - ((ρ() - z(t())) * inv(σ() * -1)) *
 identity(0)))

Thus ModelingToolkit.jl can utilize existing numerical code on symbolic codes

Let's follow this thread a little deeper.

Automatically convert numerical codes to symbolic

Let's take someone's code written to numerically solve the Lorenz equation:

function lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end
lorenz (generic function with 1 method)

Since ModelingToolkit can trace generic numerical functions in Julia, let's trace it with Operations. When we do this, it'll spit out a symbolic representation of their numerical code:

u = [x,y,z]
du = similar(u)
p = [σ,ρ,β]
lorenz(du,u,p,t)
du
3-element Array{ModelingToolkit.Operation,1}:
          σ() * (y(t()) - x(t()))
 x(t()) * (ρ() - z(t())) - y(t())
   x(t()) * y(t()) - β() * z(t())

We can then perform symbolic manipulations on their numerical code, and build a new numerical code that optimizes/fixes their original function!

J = [Dx(du[1]) Dy(du[1]) Dz(du[1])
     Dx(du[2]) Dy(du[2]) Dz(du[2])
     Dx(du[3]) Dy(du[3]) Dz(du[3])]
J = expand_derivatives.(J)
3×3 Array{ModelingToolkit.Expression,2}:
     σ() * -1           σ()  Constant(0)
 ρ() - z(t())  Constant(-1)  x(t()) * -1
       y(t())        x(t())     -1 * β()

Automated Sparsity Detection

In many cases one has to speed up large modeling frameworks by taking into account sparsity. While ModelingToolkit.jl can be used to compute Jacobians, we can write a standard Julia function in order to get a spase matrix of expressions which automatically detects and utilizes the sparsity of their function.

using SparseArrays
function SparseArrays.SparseMatrixCSC(M::Matrix{T}) where {T<:ModelingToolkit.Expression}
    idxs = findall(!iszero, M)
    I = [i[1] for i in idxs]
    J = [i[2] for i in idxs]
    V = [M[i] for i in idxs]
    return SparseArrays.sparse_IJ_sorted!(I, J, V, size(M)...)
end
sJ = SparseMatrixCSC(J)
3×3 SparseArrays.SparseMatrixCSC{ModelingToolkit.Expression,Int64} with 8 stored entries:
  [1, 1]  =  σ() * -1
  [2, 1]  =  ρ() - z(t())
  [3, 1]  =  y(t())
  [1, 2]  =  σ()
  [2, 2]  =  Constant(-1)
  [3, 2]  =  x(t())
  [2, 3]  =  x(t()) * -1
  [3, 3]  =  -1 * β()

Dependent Variables, Functions, Chain Rule

"Variables" are overloaded. When you are solving a differential equation, the variable u(t) is actually a function of time. In order to handle these kinds of variables in a mathematically correct and extensible manner, the ModelingToolkit IR actually treats variables as functions, and constant variables are simply 0-ary functions (t()).

We can utilize this idea to have parameters that are also functions. For example, we can have a parameter σ which acts as a function of 1 argument, and then utilize this function within our differential equations:

@parameters σ(..)
eqs = [D(x) ~ σ(t-1)*(y-x),
       D(y) ~ x*(σ(t^2)-z)-y,
       D(z) ~ x*y - β*z]
3-element Array{ModelingToolkit.Equation,1}:
 ModelingToolkit.Equation((D'~t())(x(t())), σ(t() - 1) * (y(t()) - x(t())))         
 ModelingToolkit.Equation((D'~t())(y(t())), x(t()) * (σ(t() ^ 2) - z(t())) - y(t()))
 ModelingToolkit.Equation((D'~t())(z(t())), x(t()) * y(t()) - β() * z(t()))

Notice that when we calculate the derivative with respect to t, the chain rule is automatically handled:

@derivatives Dₜ'~t
Dₜ(x*(σ(t^2)-z)-y)
expand_derivatives(Dₜ(x*(σ(t^2)-z)-y))
(σ(t() ^ 2) - z(t())) * (D'~t())(x(t())) + x(t()) * ((D'~t())(σ(t() ^ 2)) + -1 * (D'~t())(z(t()))) + -1 * (D'~t())(y(t()))

Hackability: Extend directly from the language

ModelingToolkit.jl is written in Julia, and thus it can be directly extended from Julia itself. Let's define a normal Julia function and call it with a variable:

_f(x) = 2x + x^2
_f(x)
2 * x(t()) + x(t()) ^ 2

Recall that when we do that, it will automatically trace this function and then build a symbolic expression. But what if we wanted our function to be a primative in the symbolic framework? This can be done by registering the function.

f(x) = 2x + x^2
@register f(x)
f (generic function with 2 methods)

Now this function is a new primitive:

f(x)
Main.WeaveSandBox22.f(x(t()))

and we can now define derivatives of our function:

function ModelingToolkit.derivative(::typeof(f), args::NTuple{1,Any}, ::Val{1})
    2 + 2args[1]
end
expand_derivatives(Dx(f(x)))
2 + 2 * x(t())

Appendix

This tutorial is part of the DiffEqTutorials.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqTutorials.jl

To locally run this tutorial, do the following commands:

using DiffEqTutorials
DiffEqTutorials.weave_file("ode_extras","01-ModelingToolkit.jmd")

Computer Information:

Julia Version 1.1.1
Commit 55e36cc308 (2019-05-16 04:10 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)

Package Information:

Status `~/.julia/environments/v1.1/Project.toml`
[7e558dbc-694d-5a72-987c-6f4ebed21442] ArbNumerics 0.5.4
[6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf] BenchmarkTools 0.4.2
[be33ccc6-a3ff-5ff2-a52e-74243cff1e17] CUDAnative 2.2.0
[159f3aea-2a34-519c-b102-8c37f9878175] Cairo 0.6.0
[3a865a2d-5b23-5a0f-bc46-62713ec82fae] CuArrays 1.0.2
[55939f99-70c6-5e9b-8bb0-5071ed7d61fd] DecFP 0.4.8
[abce61dc-4473-55a0-ba07-351d65e31d42] Decimals 0.4.0
[ebbdde9d-f333-5424-9be2-dbf1e9acfb5e] DiffEqBayes 1.1.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 3.8.2
[459566f4-90b8-5000-8ac3-15dfb0a30def] DiffEqCallbacks 2.5.2
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.9.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.6.0
[055956cb-9e8b-5191-98cc-73ae4a59e68a] DiffEqPhysics 3.1.0
[6d1b261a-3be8-11e9-3f2f-0b112a9a8436] DiffEqTutorials 0.1.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.4.0
[31c24e10-a181-5473-b8eb-7969acd0382f] Distributions 0.20.0
[497a8b3b-efae-58df-a0af-a86822472b78] DoubleFloats 0.9.1
[f6369f11-7733-5829-9624-2563aa707210] ForwardDiff 0.10.3
[c91e804a-d5a3-530f-b6f0-dfbca275c004] Gadfly 1.0.1
[7073ff75-c697-5162-941a-fcdaad2a7d2a] IJulia 1.18.1
[4138dd39-2aa7-5051-a626-17a0bb65d9c8] JLD 0.9.1
[23fbe1c1-3f47-55db-b15f-69d7ec21a316] Latexify 0.8.2
[eff96d63-e80a-5855-80a2-b1b0885c5ab7] Measurements 2.0.0
[961ee093-0014-501f-94e3-6117800e7a78] ModelingToolkit 0.2.0
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[2774e3e8-f4cf-5e23-947b-6d7e65073b56] NLsolve 4.0.0
[429524aa-4258-5aef-a3af-852621145aeb] Optim 0.18.1
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.8.1
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.1.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.25.2
[d330b81b-6aea-500a-939a-2ce795aea3ee] PyPlot 2.8.1
[731186ca-8d62-57ce-b412-fbd966d074cd] RecursiveArrayTools 0.20.0
[90137ffa-7385-5640-81b9-e52037218182] StaticArrays 0.11.0
[f3b207a7-027a-5e70-b257-86293d7955fd] StatsPlots 0.11.0
[c3572dad-4567-51f8-b174-8c6c989267f4] Sundials 3.6.1
[1986cc42-f94f-5a68-af5c-568840ba703d] Unitful 0.16.0
[44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9] Weave 0.9.1
[b77e0a4c-d291-57a0-90e8-8db25a27a240] InteractiveUtils
[37e2e46d-f89d-539d-b4ee-838fcccc9c8e] LinearAlgebra
[44cfe95a-1eb2-52ea-b672-e2afdf69b78f] Pkg