Skip to content

Commit f15e3e5

Browse files
clean up line lengths and add a copy-paste section
1 parent dc14e6f commit f15e3e5

File tree

1 file changed

+100
-66
lines changed

1 file changed

+100
-66
lines changed

docs/src/tutorials/ode_modeling.md

Lines changed: 100 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,37 @@
22

33
This is an introductory example for the usage of ModelingToolkit (MTK).
44
It illustrates the basic user-facing functionality by means of some
5-
65
examples of Ordinary Differential Equations (ODE). Some references to
76
more specific documentation are given at appropriate places.
87

8+
## Copy-Pastable Simplified Example
9+
10+
A much deeper tutorial with forcing functions and sparse Jacobians is all below.
11+
But if you want to just see some code and run, here's an example:
12+
13+
```julia
14+
using ModelingToolkit
15+
16+
@variables t, x(t) RHS(t) # independent and dependent variables
17+
@parameters τ # parameters
18+
D = Differential(t) # define an operator for the differentiation w.r.t. time
19+
20+
# your first ODE, consisting of a single equation, indicated by ~
21+
@named fol_separate = ODESystem([ RHS ~ (1 - x)/τ,
22+
D(x) ~ RHS ])
23+
24+
using DifferentialEquations: solve
25+
using Plots: plot
26+
27+
prob = ODEProblem(structural_simplify(fol_separate), [x => 0.0], (0.0,10.0), [τ => 3.0])
28+
sol = solve(prob)
29+
plot(sol, vars=[x,RHS])
30+
```
31+
32+
![Simulation result of first-order lag element, with right-hand side](https://user-images.githubusercontent.com/13935112/111958403-7e8d3e00-8aed-11eb-9d18-08b5180a59f9.png)
33+
34+
Now let's start digging into MTK!
35+
936
## Your very first ODE
1037

1138
Let us start with a minimal example. The system to be modelled is a
@@ -16,9 +43,10 @@ first-order lag element:
1643
\dot{x} = \frac{f(t) - x(t)}{\tau}
1744
```
1845

19-
Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state variable,
20-
``f(t)`` is an external forcing function, and ``\tau`` is a constant parameter. In MTK, this system can be modelled as follows. For simplicity, we first set the forcing function to
21-
a constant value.
46+
Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state
47+
variable, ``f(t)`` is an external forcing function, and ``\tau`` is a constant
48+
parameter. In MTK, this system can be modelled as follows. For simplicity, we
49+
first set the forcing function to a constant value.
2250

2351
```julia
2452
using ModelingToolkit
@@ -57,6 +85,7 @@ from the actual symbolic elements to their values, represented as an array
5785
of `Pair`s, which are constructed using the `=>` operator.
5886

5987
## Algebraic relations and structural simplification
88+
6089
You could separate the calculation of the right-hand side, by introducing an
6190
intermediate variable `RHS`:
6291

@@ -72,10 +101,11 @@ intermediate variable `RHS`:
72101
# τ
73102
```
74103

75-
To directly solve this system, you would have to create a Differential-Algebraic Equation (DAE)
76-
problem, since besides the differential equation, there is an additional algebraic equation now.
77-
However, this DAE system can obviously be transformed into the single ODE we used in the
78-
first example above. MTK achieves this by means of structural simplification:
104+
To directly solve this system, you would have to create a Differential-Algebraic
105+
Equation (DAE) problem, since besides the differential equation, there is an
106+
additional algebraic equation now. However, this DAE system can obviously be
107+
transformed into the single ODE we used in the first example above. MTK achieves
108+
this by means of structural simplification:
79109

80110
```julia
81111
fol_simplified = structural_simplify(fol_separate)
@@ -88,14 +118,16 @@ equations(fol_simplified) == equations(fol_model)
88118
# true
89119
```
90120

91-
You can extract the equations from a system using `equations` (and, in the same way,
92-
`variables` and `parameters`). The simplified equation is exactly the same as the
93-
original one, so the simulation performence will also be the same.
121+
You can extract the equations from a system using `equations` (and, in the same
122+
way, `variables` and `parameters`). The simplified equation is exactly the same
123+
as the original one, so the simulation performence will also be the same.
94124
However, there is one difference. MTK does keep track of the eliminated
95-
algebraic variables as "observables" (see [Observables and Variable Elimination](@ref)).
125+
algebraic variables as "observables" (see
126+
[Observables and Variable Elimination](@ref)).
96127
That means, MTK still knows how to calculate them out of the information available
97128
in a simulation result. The intermediate variable `RHS` therefore can be plotted
98-
along with the state variable. Note that this has to be requested explicitly, though:
129+
along with the state variable. Note that this has to be requested explicitly,
130+
though:
99131

100132
```julia
101133
prob = ODEProblem(fol_simplified, [x => 0.0], (0.0,10.0), [τ => 3.0])
@@ -111,6 +143,7 @@ values of `x` matching `sol.t`, etc. Note that this works even for variables
111143
which have been eliminated, and thus `sol[RHS]` retrieves the values of `RHS`.
112144

113145
## Specifying a time-variable forcing function
146+
114147
What if the forcing function (the "external input") ``f(t)`` is not constant?
115148
Obviously, one could use an explicit, symbolic function of time:
116149

@@ -120,9 +153,11 @@ Obviously, one could use an explicit, symbolic function of time:
120153
```
121154

122155
But often there is time-series data, such as measurement data from an experiment,
123-
we want to embedd as data in the simulation of a PDE, or as a forcing function on the right-hand side of an ODE -- is it is the case here. For this, MTK allows to "register" arbitrary
124-
Julia functions, which are excluded from symbolic transformations but are just used
125-
as-is. So, you could, for example, interpolate a given time series using
156+
we want to embed as data in the simulation of a PDE, or as a forcing function on
157+
the right-hand side of an ODE -- is it is the case here. For this, MTK allows to
158+
"register" arbitrary Julia functions, which are excluded from symbolic
159+
transformations but are just used as-is. So, you could, for example, interpolate
160+
a given time series using
126161
[DataInterpolations.jl](https://github.com/PumasAI/DataInterpolations.jl). Here,
127162
we illustrate this option by a simple lookup ("zero-order hold") of a vector
128163
of random values:
@@ -141,13 +176,11 @@ plot(sol, vars=[x,f])
141176

142177
![Simulation result of first-order lag element, step-wise forcing function](https://user-images.githubusercontent.com/13935112/111958424-83ea8880-8aed-11eb-8f42-489f4b44c3bc.png)
143178

144-
145179
## Building component-based, hierarchical models
146180

147-
Working with simple one-equation systems is already fun, but composing more complex systems from
148-
simple ones is even more fun. Best practice for such a "modeling framework" could be to use
149-
150-
factory functions for model components:
181+
Working with simple one-equation systems is already fun, but composing more
182+
complex systems from simple ones is even more fun. Best practice for such a
183+
"modeling framework" could be to use factory functions for model components:
151184

152185
```julia
153186
function fol_factory(separate=false;name)
@@ -162,17 +195,17 @@ function fol_factory(separate=false;name)
162195
end
163196
```
164197

165-
Such a factory can then used to instantiate the same component multiple times, but allows
166-
for customization:
198+
Such a factory can then used to instantiate the same component multiple times,
199+
but allows for customization:
167200

168201
```julia
169202
@named fol_1 = fol_factory()
170203
@named fol_2 = fol_factory(true) # has observable RHS
171204
```
172205

173-
Now, these two components can be used as subsystems of a parent system, i.e. one level
174-
higher in the model hierarchy. The connections between the components again are just
175-
algebraic relations:
206+
Now, these two components can be used as subsystems of a parent system, i.e.
207+
one level higher in the model hierarchy. The connections between the components
208+
again are just algebraic relations:
176209

177210
```julia
178211
connections = [ fol_1.f ~ 1.5,
@@ -193,11 +226,10 @@ connections = [ fol_1.f ~ 1.5,
193226

194227
All equations, variables and parameters are collected, but the structure of the
195228
hierarchical model is still preserved. That is, you can still get information about
196-
`fol_1` by addressing it by `connected.fol_1`, or its parameter by `connected.fol_1.τ`.
197-
Before simulation, we again eliminate the algebraic variables and connection equations
198-
from the system using structural simplification:
199-
200-
229+
`fol_1` by addressing it by `connected.fol_1`, or its parameter by
230+
`connected.fol_1.τ`. Before simulation, we again eliminate the algebraic
231+
variables and connection equations from the system using structural
232+
simplification:
201233

202234
```julia
203235
connected_simp = structural_simplify(connected)
@@ -242,8 +274,8 @@ solve(prob) |> plot
242274

243275
More on this topic may be found in [Composing Models and Building Reusable Components](@ref).
244276

245-
246277
## Defaults
278+
247279
Often it is a good idea to specify reasonable values for the initial state and the
248280
parameters of a model component. Then, these do not have to be explicitly specified when constructing the `ODEProblem`.
249281

@@ -257,20 +289,22 @@ end
257289
ODEProblem(unitstep_fol_factory(name=:fol),[],(0.0,5.0),[]) |> solve
258290
```
259291

260-
Note that the defaults can be functions of the other variables, which is then resolved
261-
at the time of the problem construction.
262-
Of course, the factory function could accept additional arguments to optionally
263-
specify the initial state or parameter values, etc.
292+
Note that the defaults can be functions of the other variables, which is then
293+
resolved at the time of the problem construction. Of course, the factory
294+
function could accept additional arguments to optionally specify the initial
295+
state or parameter values, etc.
264296

265297
## Symbolic and sparse derivatives
298+
266299
One advantage of a symbolic toolkit is that derivatives can be calculated
267-
explicitly, and that the incidence matrix of partial derivatives (the "sparsity pattern")
268-
can also be explicitly derived. These two facts lead to a substantial speedup of all
269-
model calculations, e.g. when simulating a model over time using an ODE solver.
300+
explicitly, and that the incidence matrix of partial derivatives (the
301+
"sparsity pattern") can also be explicitly derived. These two facts lead to a
302+
substantial speedup of all model calculations, e.g. when simulating a model
303+
over time using an ODE solver.
270304

271-
By default, analytical derivatives and sparse matrices, e.g. for the Jacobian, the matrix of first
272-
partial derivatives, are not used. Let's benchmark this (`prob` still is the problem using the
273-
`connected_simp` system above):
305+
By default, analytical derivatives and sparse matrices, e.g. for the Jacobian, the
306+
matrix of first partial derivatives, are not used. Let's benchmark this (`prob`
307+
still is the problem using the `connected_simp` system above):
274308

275309
```julia
276310
using BenchmarkTools
@@ -280,8 +314,8 @@ using DifferentialEquations: Rodas4
280314
# 251.300 μs (873 allocations: 31.18 KiB)
281315
```
282316

283-
Now have MTK provide sparse, analytical derivatives to the solver. This has to be
284-
specified during the construction of the `ODEProblem`:
317+
Now have MTK provide sparse, analytical derivatives to the solver. This has to
318+
be specified during the construction of the `ODEProblem`:
285319

286320
```julia
287321
prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true, sparse=true)
@@ -290,38 +324,38 @@ prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true, sparse=true)
290324
# 142.899 μs (1297 allocations: 83.96 KiB)
291325
```
292326

293-
The speedup is significant. For this small dense model (3 of 4 entries are populated), using sparse matrices is counterproductive in terms of required memory allocations. For large,
327+
The speedup is significant. For this small dense model (3 of 4 entries are
328+
populated), using sparse matrices is counterproductive in terms of required
329+
memory allocations. For large,
294330

295331
hierarchically built models, which tend to be sparse, speedup and the reduction of
296332
memory allocation can be expected to be substantial. In addition, these problem
297333
builders allow for automatic parallelism using the structural information. For
298334
more information, see the [ODESystem](@ref ODESystem) page.
299335

300336
## Notes and pointers how to go on
337+
301338
Here are some notes that may be helpful during your initial steps with MTK:
302339

303340
* Sometimes, the symbolic engine within MTK is not able to correctly identify the
304-
independent variable (e.g. time) out of all variables. In such a case, you usually
305-
get an error that some variable(s) is "missing from variable map". In most cases,
306-
it is then sufficient to specify the independent variable as second argument to
307-
`ODESystem`, e.g. `ODESystem(eqs, t)`.
308-
309-
* A completely macro-free usage of MTK is possible and is discussed in a separate tutorial. This is for package
310-
developers, since the macros are only essential for automatic symbolic naming for modelers.
311-
312-
313-
* Vector-valued parameters and variables are possible. A cleaner, more consistent treatment of
314-
these is work in progress, though. Once finished, this introductory tutorial will also
315-
cover this feature.
341+
independent variable (e.g. time) out of all variables. In such a case, you
342+
usually get an error that some variable(s) is "missing from variable map". In
343+
most cases, it is then sufficient to specify the independent variable as second
344+
argument to `ODESystem`, e.g. `ODESystem(eqs, t)`.
345+
* A completely macro-free usage of MTK is possible and is discussed in a
346+
separate tutorial. This is for package developers, since the macros are only
347+
essential for automatic symbolic naming for modelers.
348+
* Vector-valued parameters and variables are possible. A cleaner, more
349+
consistent treatment of these is work in progress, though. Once finished,
350+
this introductory tutorial will also cover this feature.
316351

317352
Where to go next?
318353

319-
* Not sure how MTK relates to similar tools and packages? Read [Comparison of ModelingToolkit vs Equation-Based Modeling Languages](@ref).
320-
321-
* Depending on what you want to do with MTK, have a look at some of the other **Symbolic Modeling Tutorials**.
322-
323-
* If you want to automatically convert an existing function to a symbolic representation, you might go through the **ModelingToolkitize Tutorials**.
324-
325-
* To learn more about the inner workings of MTK, consider the sections under **Basics** and **System Types**.
326-
```suggestion
327-
* To learn more about the details of using MTK, consider the sections under **Basics** and **System Types**.
354+
* Not sure how MTK relates to similar tools and packages? Read
355+
[Comparison of ModelingToolkit vs Equation-Based Modeling Languages](@ref).
356+
* Depending on what you want to do with MTK, have a look at some of the other
357+
**Symbolic Modeling Tutorials**.
358+
* If you want to automatically convert an existing function to a symbolic
359+
representation, you might go through the **ModelingToolkitize Tutorials**.
360+
* To learn more about the inner workings of MTK, consider the sections under
361+
**Basics** and **System Types**.

0 commit comments

Comments
 (0)