|
| 1 | +# [Building ODEs and DAEs with ModelingToolkit.jl](@id getting_started_ode) |
| 2 | + |
| 3 | +This is an introductory tutorial for ModelingToolkit.jl (MTK). We will demonstrate the |
| 4 | +basics of the package by demontrating how to build systems of Ordinary Differential |
| 5 | +Equations (ODEs) and Differential-Algebraic Equations (DAEs). |
| 6 | + |
| 7 | +## Installing ModelingToolkit |
| 8 | + |
| 9 | +To install ModelingToolkit, use the Julia package manager. This can be done as follows: |
| 10 | + |
| 11 | +```julia |
| 12 | +using Pkg |
| 13 | +Pkg.add("ModelingToolkit") |
| 14 | +``` |
| 15 | + |
| 16 | +## The end goal |
| 17 | + |
| 18 | +TODO |
| 19 | + |
| 20 | +## Basics of MTK |
| 21 | + |
| 22 | +ModelingToolkit.jl is a symbolic-numeric system. This means it allows specifying a model |
| 23 | +(such as an ODE) in a similar way to how it would be written on paper. Let's start with a |
| 24 | +simple example. The system to be modeled is a first-order lag element: |
| 25 | + |
| 26 | +```math |
| 27 | +\dot{x} = \frac{f(t) - x(t)}{\tau} |
| 28 | +``` |
| 29 | + |
| 30 | +Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) unknown |
| 31 | +variable, ``f(t)`` is an external forcing function, and ``\tau`` is a |
| 32 | +parameter. |
| 33 | + |
| 34 | +For simplicity, we will start off by setting the forcing function to a constant `1`. Every |
| 35 | +ODE has a single independent variable. MTK has a common definition for time `t` and the |
| 36 | +derivative with respect to it. |
| 37 | + |
| 38 | +```@example ode2 |
| 39 | +using ModelingToolkit |
| 40 | +using ModelingToolkit: t_nounits as t, D_nounits as D |
| 41 | +``` |
| 42 | + |
| 43 | +Next, we declare the (dependent) variables and the parameters of our model: |
| 44 | + |
| 45 | +```@example ode2 |
| 46 | +@variables x(t) |
| 47 | +@parameters τ |
| 48 | +``` |
| 49 | + |
| 50 | +Note the syntax `x(t)`. We must declare that the variable `x` is a function of the independent |
| 51 | +variable `t`. Next, we define the equations of the system: |
| 52 | + |
| 53 | +```@example ode2 |
| 54 | +eqs = [D(x) ~ (1 - x) / τ] |
| 55 | +``` |
| 56 | + |
| 57 | +Since `=` is reserved as the assignment operator, MTK uses `~` to denote equality between |
| 58 | +expressions. Now we must consolidate all of this information about our system of ODEs into |
| 59 | +ModelingToolkit's `System` type. |
| 60 | + |
| 61 | +```@example ode2 |
| 62 | +sys = System(eqs, t, [x], [τ]; name = :sys) |
| 63 | +``` |
| 64 | + |
| 65 | +The `System` constructor accepts a `Vector{Equation}` as the first argument, followed by the |
| 66 | +independent variable, a list of dependent variables, and a list of parameters. Every system |
| 67 | +must be given a name via the `name` keyword argument. Most of the time, we want to name our |
| 68 | +system the same as the variable it is assigned to. The `@named` macro helps with this: |
| 69 | + |
| 70 | +```@example ode2 |
| 71 | +@named sys = System(eqs, t, [x], [τ]) |
| 72 | +``` |
| 73 | + |
| 74 | +Additionally, it may become inconvenient to specify all variables and parameters every time |
| 75 | +a system is created. MTK allows omitting these arguments, and will automatically infer them |
| 76 | +from the equations. |
| 77 | + |
| 78 | +```@example ode2 |
| 79 | +@named sys = System(eqs, t) |
| 80 | +``` |
| 81 | + |
| 82 | +Our system is not quite ready for simulation yet. First, we must use the `mtkcompile` |
| 83 | +function which transforms the system into a form that MTK can handle. For our trivial |
| 84 | +system, this does not do much. |
| 85 | + |
| 86 | +```@example ode2 |
| 87 | +simp_sys = mtkcompile(sys) |
| 88 | +``` |
| 89 | + |
| 90 | +Since building and simplifying a system is a common workflow, MTK provides the `@mtkcompile` |
| 91 | +macro for convenience. |
| 92 | + |
| 93 | +```@example ode2 |
| 94 | +@mtkcompile sys = System(eqs, t) |
| 95 | +``` |
| 96 | + |
| 97 | +We can now build an `ODEProblem` from the system. ModelingToolkit generates the necessary |
| 98 | +code for numerical ODE solvers to solve this system. We need to provide an initial |
| 99 | +value for the variable `x` and a value for the parameter `p`, as well as the time span |
| 100 | +for which to simulate the system. |
| 101 | + |
| 102 | +```@example ode2 |
| 103 | +prob = ODEProblem(sys, [x => 0.0, τ => 3.0], (0.0, 10.0)) |
| 104 | +``` |
| 105 | + |
| 106 | +Here, we are saying that `x` should start at `0.0`, `τ` should be `3.0` and the system |
| 107 | +should be simulated from `t = 0.0` to `t = 10.0`. To solve the system, we must import a |
| 108 | +solver. |
| 109 | + |
| 110 | +```@example ode2 |
| 111 | +using OrdinaryDiffEq |
| 112 | +
|
| 113 | +sol = solve(prob) |
| 114 | +``` |
| 115 | + |
| 116 | +[OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/) contains a large number of |
| 117 | +numerical solvers. It also comes with a default solver which is used when calling |
| 118 | +`solve(prob)` and is capable of handling a large variety of systems. |
| 119 | + |
| 120 | +We can obtain the timeseries of `x` by indexing the solution with the symbolic variable: |
| 121 | + |
| 122 | +```@example ode2 |
| 123 | +sol[x] |
| 124 | +``` |
| 125 | + |
| 126 | +We can even obtain timeseries of complicated expressions involving the symbolic variables in |
| 127 | +the model |
| 128 | + |
| 129 | +```@example ode2 |
| 130 | +sol[(1 - x) / τ] |
| 131 | +``` |
| 132 | + |
| 133 | +Perhaps more interesting is a plot of the solution. This can easily be achieved using Plots.jl. |
| 134 | + |
| 135 | +```@example ode2 |
| 136 | +using Plots |
| 137 | +
|
| 138 | +plot(sol) |
| 139 | +``` |
| 140 | + |
| 141 | +Similarly, we can plot different expressions: |
| 142 | + |
| 143 | +```@example ode2 |
| 144 | +plot(sol; idxs = (1 - x) / τ) |
| 145 | +``` |
0 commit comments