Build and run your first simulation with Julia's SciML

In this tutorial, we will build and run our first simulation with SciML!

To build our simulation, we will use the ModelingToolkit system for modeling and simulation. ModelingToolkit is a bit higher level than directly defining code for a differential equation system: it's a symbolic system that will automatically simplify our models, optimize our code, and generate compelling visualizations. Sounds neat? Let's dig in.

Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

Our Problem: Simulate the Lotka-Volterra Predator-Prey Dynamics

The dynamics of our system are given by the Lotka-Volterra dynamical system: Let $x(t)$ be the number of rabbits in the environment and $y(t)$ be the number of wolves. The equation that defines the evolution of the species is given as follows:

\[\begin{align} \frac{dx}{dt} &= \alpha x - \beta x y\\ \frac{dy}{dt} &= -\gamma y + \delta x y \end{align}\]

where $\alpha, \beta, \gamma, \delta$ are parameters. Starting from equal numbers of rabbits and wolves, $x(0) = 1$ and $y(0) = 1$, we want to simulate this system from time $t_0 = 0$ to $t_f = 10$. Luckily, a local guide provided us with some parameters that seem to match the system! These are $\alpha = 1.5$, $\beta = 1.0$, $\gamma = 3.0$, $\delta = 1.0$. How many rabbits and wolves will there be 10 months from now? And if z = x + y, i.e. the total number of animals at a given time, can we visualize this total number of animals at each time?

Solution as Copy-Pastable Code

using ModelingToolkit, DifferentialEquations, Plots

# Define our state variables: state(t) = initial condition
@variables t x(t)=1 y(t)=1 z(t)=2

# Define our parameters
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0

# Define our differential: takes the derivative with respect to `t`
D = Differential(t)

# Define the differential equations
eqs = [D(x) ~ α * x - β * x * y
       D(y) ~ -γ * y + δ * x * y
       z ~ x + y]

# Bring these pieces together into an ODESystem with independent variable t
@mtkbuild sys = ODESystem(eqs, t)

# Convert from a symbolic to a numerical problem to simulate
tspan = (0.0, 10.0)
prob = ODEProblem(sys, [], tspan)

# Solve the ODE
sol = solve(prob)

# Plot the solution
p1 = plot(sol, title = "Rabbits vs Wolves")
p2 = plot(sol, idxs = z, title = "Total Animals")

plot(p1, p2, layout = (2, 1))

Step-by-Step Solution

Step 1: Install and Import the Required Packages

To do this tutorial, we will need a few components:

• ModelingToolkit.jl, our modeling environment
• DifferentialEquations.jl, the differential equation solvers
• Plots.jl, our visualization tool

To start, let's add these packages as demonstrated in the installation tutorial:

using Pkg
Pkg.add(["ModelingToolkit", "DifferentialEquations", "Plots"])

Now we're ready. Let's load in these packages:

using ModelingToolkit, DifferentialEquations, Plots

Step 2: Define our ODE Equations

Now let's define our ODEs. We use the ModelingToolkit.@variabes statement to declare our variables. We have the independent variable time t, and then define our 3 state variables:

# Define our state variables: state(t) = initial condition
@variables t x(t)=1 y(t)=1 z(t)=2

\[ \begin{equation} \left[ \begin{array}{c} t \\ x\left( t \right) \\ y\left( t \right) \\ z\left( t \right) \\ \end{array} \right] \end{equation} \]

Notice here that we use the form state = default, where on the right-hand side the default value of a state is interpreted to be its initial condition. This is then done similarly for parameters, where the default value is now the parameter value:

# Define our parameters
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0

\[ \begin{equation} \left[ \begin{array}{c} \alpha \\ \beta \\ \gamma \\ \delta \\ \end{array} \right] \end{equation} \]

Next, we define our set of differential equations. To define the Differential operator D, we need to first tell it what to differentiate with respect to, here the independent variable t, Then, once we have the operator, we apply that into the equations.

# Define our differential: takes the derivative with respect to `t`
D = Differential(t)

# Define the differential equations
eqs = [D(x) ~ α * x - β * x * y
       D(y) ~ -γ * y + δ * x * y
       z ~ x + y]

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

Notice that in the display, it will automatically generate LaTeX. If one is interested in generating this LaTeX locally, one can simply do:

using Latexify # add the package first
latexify(eqs)

Step 3: Define the ODEProblem

Now we bring these pieces together. In ModelingToolkit, we can bring these pieces together to represent an ODESystem with the following:

# Bring these pieces together into an ODESystem with independent variable t
@mtkbuild sys = ODESystem(eqs, t)

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= x\left( t \right) \alpha - x\left( t \right) y\left( t \right) \beta \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= - y\left( t \right) \gamma + x\left( t \right) y\left( t \right) \delta \end{align} \]

Notice that in our equations we have an algebraic equation z ~ x + y. This is not a differential equation but an algebraic equation, and thus we call this set of equations a Differential-Algebraic Equation (DAE). The symbolic system of ModelingToolkit can eliminate such equations to return simpler forms to numerically approximate.

Notice that what is returned is an ODESystem, but now with the simplified set of equations. z has been turned into an “observable”, i.e. a state that is not computed but can be constructed on-demand. This is one of the ways that SciML reaches its speed: you can have 100,000 equations, but solve only 1,000 to then automatically reconstruct the full set. Here, it's just 3 equations to 2, but as models get more complex, the symbolic system will find ever more clever interactions!

Now that we have simplified our system, let's turn it into a numerical problem to approximate. This is done with the ODEProblem constructor, that transforms it from a symbolic ModelingToolkit representation to a numerical DifferentialEquations representation. We need to tell it the numerical details now:

1. Whether to override any of the default values for the initial conditions and parameters.
2. What is the initial time point.
3. How long to integrate it for.

In this case, we will use the default values for all our variables, so we will pass a blank override []. If for example we did want to change the initial condition of x to 2.0 and α to 4.0, we would do [x => 2.0, α => 4.0]. Then secondly, we pass a tuple for the time span, (0.0,10.0) meaning start at 0.0 and end at 10.0. This looks like:

# Convert from a symbolic to a numerical problem to simulate
tspan = (0.0, 10.0)
prob = ODEProblem(sys, [], tspan)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

Step 4: Solve the ODE System

Now we solve the ODE system. Julia's SciML solvers have a defaulting system that can automatically determine an appropriate solver for a given system, so we can just tell it to solve:

# Solve the ODE
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 34-element Vector{Float64}:
  0.0
  0.0776084743154256
  0.2326451370670694
  0.42911851563726466
  0.679082199936808
  0.9444046279774128
  1.2674601918628516
  1.61929140093895
  1.9869755481702074
  2.2640903679981617
  ⋮
  7.5848624442719235
  7.978067891667038
  8.483164641366145
  8.719247691882519
  8.949206449510513
  9.200184762926114
  9.438028551201125
  9.711807820573478
 10.0
u: 34-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0454942346944578, 0.8576684823217128]
 [1.1758715885890039, 0.6394595702308831]
 [1.419680958026516, 0.45699626144050703]
 [1.8767193976262215, 0.3247334288460738]
 [2.5882501035146133, 0.26336255403957304]
 [3.860709084797009, 0.27944581878759106]
 [5.750813064347339, 0.5220073551361045]
 [6.814978696356636, 1.917783405671627]
 [4.392997771045279, 4.194671543390719]
 ⋮
 [2.6142510825026886, 0.2641695435004172]
 [4.241070648057757, 0.30512326533052475]
 [6.79112182569163, 1.13452538354883]
 [6.265374940295053, 2.7416885955953294]
 [3.7807688120520893, 4.431164521488331]
 [1.8164214705302744, 4.064057991958618]
 [1.146502825635759, 2.791173034823897]
 [0.955798652853089, 1.623563316340748]
 [1.0337581330572414, 0.9063703732075853]

Step 5: Visualize the Solution

Now let's visualize the solution! Notice that our solution only has two states. If we recall, the simplified system only has two states: z was symbolically eliminated. We can access any of the values, even the eliminated values, using the symbolic variable as the index. For example:

sol[z]

34-element Vector{Float64}:
 2.0
 1.9031627170161705
 1.815331158819887
 1.876677219467023
 2.2014528264722952
 2.8516126575541865
 4.1401549035846
 6.272820419483444
 8.732762102028264
 8.587669314435999
 ⋮
 2.8784206260031056
 4.546193913388282
 7.92564720924046
 9.007063535890381
 8.21193333354042
 5.880479462488893
 3.937675860459656
 2.5793619691938368
 1.9401285062648266

returns the time series of the observable z at time points corresponding to sol.t. We can use this with the automated plotting functionality. First let's create a plot of x and y over time using plot(sol) which will plot all of the states. Then next, we will explicitly tell it to make a plot with the index being z, i.e. idxs=z.

# Plot the solution
p1 = plot(sol, title = "Rabbits vs Wolves")

p2 = plot(sol, idxs = z, title = "Total Animals")

Finally, let's make a plot where we merge these two plot elements. To do so, we can take our two plot objects, p1 and p2, and make a plot with both of them. Then we tell Plots to do a layout of (2,1), or 2 rows and 1 columns. Let's see what happens when we bring these together:

plot(p1, p2, layout = (2, 1))

And tada, we have a full analysis of our ecosystem!

Bonus Step: Emoji Variables

If you made it this far, then congrats, you get to learn a fun fact! Since Julia code can use Unicode, emojis work for variable names. Here's the simulation using emojis of rabbits and wolves to define the system:

using ModelingToolkit, DifferentialEquations
@parameters α=1.5 β=1.0 γ=3.0 δ=1.0
@variables t 🐰(t)=1.0 🐺(t)=1.0
D = Differential(t)
eqs = [D(🐰) ~ α * 🐰 - β * 🐰 * 🐺,
    D(🐺) ~ -γ * 🐺 + δ * 🐰 * 🐺]

@mtkbuild sys = ODESystem(eqs, t)
prob = ODEProblem(sys, [], (0.0, 10.0))
sol = solve(prob)

retcode: Success
Interpolation: 3rd order Hermite
t: 34-element Vector{Float64}:
  0.0
  0.0776084743154256
  0.2326451370670694
  0.42911851563726466
  0.679082199936808
  0.9444046279774128
  1.2674601918628516
  1.61929140093895
  1.9869755481702074
  2.2640903679981617
  ⋮
  7.5848624442719235
  7.978067891667038
  8.483164641366145
  8.719247691882519
  8.949206449510513
  9.200184762926114
  9.438028551201125
  9.711807820573478
 10.0
u: 34-element Vector{Vector{Float64}}:
 [1.0, 1.0]
 [1.0454942346944578, 0.8576684823217128]
 [1.1758715885890039, 0.6394595702308831]
 [1.419680958026516, 0.45699626144050703]
 [1.8767193976262215, 0.3247334288460738]
 [2.5882501035146133, 0.26336255403957304]
 [3.860709084797009, 0.27944581878759106]
 [5.750813064347339, 0.5220073551361045]
 [6.814978696356636, 1.917783405671627]
 [4.392997771045279, 4.194671543390719]
 ⋮
 [2.6142510825026886, 0.2641695435004172]
 [4.241070648057757, 0.30512326533052475]
 [6.79112182569163, 1.13452538354883]
 [6.265374940295053, 2.7416885955953294]
 [3.7807688120520893, 4.431164521488331]
 [1.8164214705302744, 4.064057991958618]
 [1.146502825635759, 2.791173034823897]
 [0.955798652853089, 1.623563316340748]
 [1.0337581330572414, 0.9063703732075853]

Now go make your professor mad that they have to grade a fully emojified code. I'll vouch for you: the documentation told you to do this.