Modeling Toolkit and Symbolics

Up to now, when we needed to solve a system of ODEs or non-linear equations, for example, we have programmed these and passed them to a solver. If you are comfortable with writing code, this is not a problem. There is however a few tools available in Julia that can make the task a lot easier, both for experienced and novice programmers - a modern Computer Algebra System, Symbolics.jl.

Symbolics.jl and the Symbolics Ecosystem

Symbolics allows you to simply specify the equations for your system and then do symbolic manipulations and simplifications to the system. Since Symbolics is written in pure Julia, it interacts seamlessly with normal Julia code and you can do anything with a symbolic variable you would do with any other Julia variable. This gives great flexibility. It doesn’t stop there however - it will also generate a Julia function for you - optimised and parallelised, if specified.

In addition, there are several other packages that build on Symbolics:

  • ModelingToolkit: Provides the tools to symbolically specify common numerical systems like ODEs, PDEs, non-linear systems, control problems, causal and acausal modelling etc.
  • Catalyst: Provides tools for specifying reaction networks and kinetics
  • DataDrivenDiffEq: Automatic identifications of ODEs/DAEs from data
  • SymbolicRegression: Genetic programming to find equations from data
Note

ModelingToolkit is developing rapidly, which is great in terms of new features being added all the time. Unfortunately this also means syntax improvements are happening all the time, so you may need to check the documentation from time to time to catch up with new developments and new syntax.

Using Symbolics

This is a very brief introduction to Symbolics. For more detail, refer to the manual.

At the simplest level, Symbolics works with symbolic variables and equations. You specify equations with the @variables macro and then simply build equations with these symbolic variables. You can simplify equations, using simplify(), define and calculate derivatives and Jacobians etc. If you are using Symbolics in a Jupyter notebook or in Pluto, the equations will be rendered as \(LaTeX\). In the REPL or a script, there will be text output. You can also generate \(LaTeX\) expressions with the Latexify package.

using Symbolics

@variables x, y
# 2-element Vector{Num}:
#  x
#  y

z = (x + y)*(x - y)

Dx = Differential(x)
# (::Differential) (generic function with 2 methods)

Dy = Differential(y)
# (::Differential) (generic function with 2 methods)

expand_derivatives(Dx(z))
# 2x

expand_derivatives(Dy(z))
# -2y

simplify(2x + x^2 - y - x -2x^2)
# x - y - (x^2)

A = [x x*y y^2;
     y 2*x*y x^2*y]

# 2×3 Matrix{Num}:
#  x   x*y      y^2
#  y  2x*y  y*(x^2)

latexify(A)
# L"\begin{equation}
# \left[
# \begin{array}{ccc}
# x & x y & y^{2} \\
# y & 2 x y & x^{2} y \\
# \end{array}
# \right]
# \end{equation}
# "

The last result, when rendered in \(LaTeX\), or in a Jupyter notebook, for example, looks like this:

\[ \begin{equation} \left[ \begin{array}{ccc} x & x y & y^{2} \\ y & 2 x y & x^{2} y \\ \end{array} \right] \end{equation} \]

So, we can do algebra and calculus and make pretty equations for reports. What else can we do? Well, a whole lot, as it turns out. Symbolics is a full-featured CAS (Computer Algebra System). Whats more, since it is fully written in Julia, it interacts with the rest of Julia in a pretty seamless manner. Anything you would do with a Julia variable, you can also do with a symbolic (type Num) variable.

For example, we can call a normal, generic function:

f(u) = 2*u[1] - u[2]
# f (generic function with 1 method)

f([x, y])
# 2x - y

We can also define more complicated systems, where the variables are dependent on each other:

using Symbolics, Latexify

@variables t x(t) y(t) # Declare an unknown dependency on t for x and y
# 3-element Vector{Num}:
#     t
#  x(t)
#  y(t)

z = x*t + y*t^2
# t*x(t) + (t^2)*y(t)

Dt = Differential(t)
# (::Differential) (generic function with 2 methods)

expand_derivatives(Dt(z))
# t*Differential(t)(x(t)) + (t^2)*Differential(t)(y(t)) + 2t*y(t) + x(t)

latexify(expand_derivatives(Dt(z)))
# L"\begin{equation}
# t \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} + t^{2} \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} + 2 t y\left( t \right) + x\left( t \right)
# \end{equation}
# "

Where the result from the last line is: \[ \begin{equation} t \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} + t^{2} \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} + 2 t y\left( t \right) + x\left( t \right) \end{equation} \]

In other words, since we indicated that x and y depend on t, this information was included in the derivative calculation.

The final, but most important part of Symbolics is its ability to convert symbolic equations into Julia code.


using Symbolics

@variables x y
# 2-element Vector{Num}:
#  x
#  y

eqs = [x^2 + y, y^2 + x]
# 2-element Vector{Num}:
#  y + x^2
#  x + y^2)

f_expr = build_function(eqs, [x, y])
Base.remove_linenums!.(f_expr)
# (:(function (ˍ₋arg1,)
#       begin
#           begin
#               (SymbolicUtils.Code.create_array)(typeof(ˍ₋arg1), nothing, Val{1}(), Val{(2,)}(), (+)(ˍ₋arg1[2], (^)(ˍ₋arg1[1], 2)), (+)(ˍ₋arg1[1], (^)(ˍ₋arg1[2], 2)))
#           end
#       end
#   end), :(function (ˍ₋out, ˍ₋arg1)
#       begin
#           begin
#               @inbounds begin
#                   ˍ₋out[1] = (+)(ˍ₋arg1[2], (^)(ˍ₋arg1[1], 2))
#                   ˍ₋out[2] = (+)(ˍ₋arg1[1], (^)(ˍ₋arg1[2], 2))
#                   nothing
#               end
#           end
#       end
#   end))
Note

The Base.remove_linenums!() call just removes a lot of additional line number information form the result to make it more readable. You would not normally bother with this.

Don’t worry if this looks unlike anything related to Julia code you have ever seen. What is returned is a Julia expression string - this is how Julia code is stored internally. It is also what you would generate if you wrote a Julia macro (meta-programming). A macro in Julia is a piece of code that, instead of returning a value, returns another piece of code that is then evaluated. Macros in Julia are normally identified by the @ symbol, e.g. @variables. This is a topic for another chapter.

So, what is going on here? The first thing to note is that the expression contains a tuple that holds two functions. The first starts with function (ˍ₋arg1,), while the second starts with function (ˍ₋out, ˍ₋arg1). If you look more closely at the code, you will see that the functions (we’ll call then func1 and func2! here) are equivalent to:

function func1(u)
    return u[2] + u[1]^2, u[1] + u[2]^2
end 

function func2!(z, u)
    z[1] = u[2] + u[1]^2
    z[2] = u[1] + u[2]^2
end

Both functions will do the same calculations, but func1 returns the values, which may require allocating an array/tuple to store it in, while func2! mutates the first parameter to return the value. These are simply two version of the same function. The first is called an out-of-place function and the second, an in-place function. Sound familiar? We saw these in the chapter on differential equations.

How do we use them? We need to tell Julia to evaluate the expression to generate the code:

myf1 = eval(f_expr[1])
# #7 (generic function with 1 method)

myf1([1, 2])
# 2-element Vector{Int64}:
#  3
#  5

myf2! = eval(f_expr[2])
# #9 (generic function with 1 method)

z = zeros(Int, 2)
# 2-element Vector{Int64}:
#  0
#  0

myf2!(z, [1, 2])

z
# 2-element Vector{Int64}:
#  3
#  5

This is the key functionality from Symbolics we want to use. But we won’t use it directly. Rather, we shall use the ModelingToolkit package, that builds on Symbolics and makes it easy to specify systems of equations, differential or non-linear, and then solve them with numerical solvers.

Modeling Toolkit

Non-linear system

Let’s start with an example of a simple set of non-linear equations:

\[y = 2x^{2} - 7\] \[y = 10\sqrt{x}\]

using ModelingToolkit, NonlinearSolve

@variables x y
@parameters a b c

eqs = [y ~ a*x^2 - b, y ~ c*√x]

@named nlsys = NonlinearSystem(eqs, [x, y], [a, b, c])

initial_guess = [x => 3.5, y => 20]
params = [a => 2, b => 7, c => 10]

prob = NonlinearProblem(nlsys, initial_guess, params)
sol = solve(prob, NewtonRaphson())

# u: 2-element Vector{Float64}:
#   3.6045577280701777
#  18.98567282997952

Here we load two packages: ModelingToolkit, to set up the problem, and NonlinearSolve, to solve it.

Using the Symbolics package (this is re-exported by ModelingToolkit - you don’t need to explicitly load Symbolics), we specify the variables and parameters for the problem. We could hard code the values of a, b and c, but specifying them as parameters allows to change their values easily later.

We specify the system of equations as an array, eqs. Note that in each symbolic equation, we use ~ instead of =. This is simply because Julia has rather strong opinions about you overloading the assignment operator, =, in a macro.

Now, there is one new step: we create a NonlinearSystem. This is effectively the step were we ask ModelingToolkit to create the code for our system of non-linear equations. We pass the array of equations, an array with the variables and an array with the parameters. The latter can also be an empty array, if you hard-coded the parameter values.

You may wonder why there is @named in front of this line…

We need a variable to refer to the object created by NonlinearSystem. Internally, however ModelingToolkit needs to also know what the name of this variable is. You could handle this by adding a keyword argument name with the same name as the variable (expressed as a Symbol), or you can just use the convenient macro, @named, that will handle it for you.

Next, we create a NonlinearProblem, just like we did in the chapter on solving non-linear systems. Instead of passing a function that encodes the system of equations, we pass the object created by NonlinearSystem. Then we solve the problem, exactly like before.

Solving ODEs

The example above was a fairly straightforward problem. ModelingToolkit can however do a lot more for us. To illustrate some of the magic, we’ll look at an example from the ModelingToolkit documentation.

Here we have a system three of ODEs, the first of which is a second order ODE.

\[\begin{align} \frac{\mathrm{d}}{\mathrm{d}t} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \sigma \left( - x\left( t \right) + y\left( t \right) \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& - y\left( t \right) + \left( \rho - z\left( t \right) \right) x\left( t \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& x\left( t \right) y\left( t \right) - \beta z\left( t \right) \end{align}\]

using ModelingToolkit, DifferentialEquations, Plots

@parameters σ ρ β
@variables t x(t) y(t) z(t)
D = Differential(t)

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

@named sys = ODESystem(eqs)

sys = structural_simplify(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())

plot(sol, idxs = (x, y))

Just like before, we specify the variables and parameters. We then, again, specify an array of equations. As discussed above in the section on Symbolics, we use Differential to specify the derivatives. A second (or higher) derivative is simply a nested call of this function.

Like in our previous example, we now need to ask ModelingToolkit to build the code for our system of ODEs. This is done with a call to ODESystem. This time, the call to ODESystem looks a lot simpler. Where are the specifications for the variables and parameters?

As it happens, the full call would be:

@named sys = ODESystem(eqs, t, [x, y, z], [σ, ρ, β])

The parameters are the equations, the independent variables, the dependent variables and the parameters. ModelingToolkit is clever enough to figure out which is which for ODESystems. This doesn’t work for NonlinearSystems yet, but one would image this functionality is on its way.

The next part is where the magic happens:

sys = structural_simplify(sys)

We call the all-singing, all-dancing structural_simplify(), which is really the heart of ModelingToolkit. It will do a whole host of transformations of your system, using everything available in Symbolics. Let’s look at the intermediate results to understand what it is doing for us:

@named sys = ODESystem(eqs, t, [x, y, z], [σ, ρ, β])

# Model sys with 3 equations
# States (3):
#   x(t)
#   y(t)
#   z(t)
# Parameters (3):
#   σ
#   ρ
#   β

sys = structural_simplify(sys)

# Model sys with 4 equations
# States (4):
#   y(t)
#   z(t)
#   x(t)
#   xˍt(t)       
# Parameters (3):
#   σ
#   ρ
#   β
# Incidence matrix:4×9 SparseArrays.SparseMatrixCSC{Num, Int64} with 13 stored entries:
#  ×  ×  ×  ⋅  ⋅  ⋅  ×  ⋅  ⋅
#  ×  ×  ×  ⋅  ⋅  ⋅  ⋅  ×  ⋅
#  ⋅  ⋅  ⋅  ×  ⋅  ×  ⋅  ⋅  ⋅
#  ×  ⋅  ×  ⋅  ⋅  ⋅  ⋅  ⋅  ×

equations(sys)

# 4-element Vector{Equation}:
#  Differential(t)(y(t)) ~ (ρ - z(t))*x(t) - y(t)
#  Differential(t)(z(t)) ~ x(t)*y(t) - β*z(t)
#  Differential(t)(x(t)) ~ xˍt(t)
#  Differential(t)(xˍt(t)) ~ xˍtt(t)

From the results of each line, you can see that structural_simplify() has converted our system of two first- and one second-order ODEs into a system of four first order ODEs. We can inspect these equations with the equations() function, and get the pretty versions with latexify():

\[\begin{align} \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} =& - y\left( t \right) + \left( \rho - z\left( t \right) \right) x\left( t \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& x\left( t \right) y\left( t \right) - \beta z\left( t \right) \\ \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& xˍt\left( t \right) \\ \frac{\mathrm{d} xˍt\left( t \right)}{\mathrm{d}t} =& xˍtt\left( t \right) \end{align}\]

Wait a minute! Aren’t we missing some details on the last two equations? What is x_tt(t) defined as? For reasons that are as clear as mud, you need to also call observed(sys) to see the missing bits:

observed(sys)
# 1-element Vector{Equation}:
#  xˍtt(t) ~ σ*(y(t) - x(t))

The rest of the code should be familiar. We specify the initial values and the time span for the integration, create an ODEProblem and call the solver. The only part that will be new is in the call to ODEProblem:

prob = ODEProblem(sys, u0, tspan, p, jac = true)

We have added the keyword parameter jac = true. This tells the ODE solver that we have also specified the Jacobian matrix of the system, which was done automatically, and symbolically, for you by ModelingToolkit. This can greatly enhance the solving of difficult systems.

Then we solve the problem and finally, we plot the results:

This was a very short introduction to Symbolics and ModelingToolkit. There is a lot more to learn in the documentation, including several tutorials.

Back to top