33ModelingToolkit has a system for transformations of mathematical
44systems. These transformations allow for symbolically changing
55the representation of the model to problems that are easier to
6- numerically solve. One simple to demonstrate transformation is the
6+ numerically solve. One simple to demonstrate transformation, is
77` structural_simplify ` , which does a lot of tricks, one being the
88transformation that turns an Nth order ODE into N
99coupled 1st order ODEs.
@@ -15,16 +15,28 @@ We utilize the derivative operator twice here to define the second order:
1515using ModelingToolkit, OrdinaryDiffEq
1616using ModelingToolkit: t_nounits as t, D_nounits as D
1717
18- @parameters σ ρ β
19- @variables x(t) y(t) z(t)
20-
21- eqs = [D(D(x)) ~ σ * (y - x),
22- D(y) ~ x * (ρ - z) - y,
23- D(z) ~ x * y - β * z]
24-
25- @named sys = ODESystem(eqs, t)
18+ @mtkmodel SECOND_ORDER begin
19+ @parameters begin
20+ σ = 28.0
21+ ρ = 10.0
22+ β = 8 / 3
23+ end
24+ @variables begin
25+ x(t) = 1.0
26+ y(t) = 0.0
27+ z(t) = 0.0
28+ end
29+ @equations begin
30+ D(D(x)) ~ σ * (y - x)
31+ D(y) ~ x * (ρ - z) - y
32+ D(z) ~ x * y - β * z
33+ end
34+ end
35+ @mtkbuild sys = SECOND_ORDER()
2636```
2737
38+ The second order ODE has been automatically transformed to two first order ODEs.
39+
2840Note that we could've used an alternative syntax for 2nd order, i.e.
2941` D = Differential(t)^2 ` and then ` D(x) ` would be the second derivative,
3042and this syntax extends to ` N ` -th order. Also, we can use ` * ` or ` ∘ ` to compose
@@ -33,28 +45,17 @@ and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compos
3345Now let's transform this into the ` ODESystem ` of first order components.
3446We do this by calling ` structural_simplify ` :
3547
36- ``` @example orderlowering
37- sys = structural_simplify(sys)
38- ```
39-
4048Now we can directly numerically solve the lowered system. Note that,
4149following the original problem, the solution requires knowing the
42- initial condition for ` x' ` , and thus we include that in our input
43- specification:
50+ initial condition for both ` x ` and ` D(x) ` .
51+ The former already got assigned a default value in the ` @mtkmodel ` ,
52+ but we still have to provide a value for the latter.
4453
4554``` @example orderlowering
46- u0 = [D(x) => 2.0,
47- x => 1.0,
48- y => 0.0,
49- z => 0.0]
50-
51- p = [σ => 28.0,
52- ρ => 10.0,
53- β => 8 / 3]
54-
55+ u0 = [D(sys.x) => 2.0]
5556tspan = (0.0, 100.0)
56- prob = ODEProblem(sys, u0, tspan, p , jac = true)
57+ prob = ODEProblem(sys, u0, tspan, [] , jac = true)
5758sol = solve(prob, Tsit5())
58- using Plots;
59- plot(sol, idxs = (x, y));
59+ using Plots
60+ plot(sol, idxs = (sys. x, sys. y))
6061```
0 commit comments