@@ -23,11 +23,8 @@ D = Differential(t)
23
23
@parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
24
24
M1 = ODESystem([
25
25
D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
26
- ], t; defaults = [
27
- y => 0.0 # start on the ground
28
- ], initialization_eqs = [
29
- #x ~ 0.0, # TODO: handle? # hide
30
- D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
26
+ ], t; initialization_eqs = [
27
+ D(x) ~ D(y) # equal initial horizontal and vertical velocity (45°)
31
28
], name = :M)
32
29
M1s = structural_simplify(M1)
33
30
```
@@ -43,15 +40,26 @@ We will demonstrate the last method by changing the independent variable from $t
43
40
This transformation is well-defined for any non-zero horizontal velocity $v$.
44
41
``` @example changeivar
45
42
M2 = change_independent_variable(M1, x; dummies = true)
46
- @assert M2.x isa Num # hide
47
43
M2s = structural_simplify(M2; allow_symbolic = true)
44
+ # a sanity test on the 10 x/y variables that are accessible to the user # hide
45
+ @assert allequal([x, M1s.x]) # hide
46
+ @assert allequal([M2.x, M2s.x]) # hide
47
+ @assert allequal([y, M1s.y]) # hide
48
+ @assert allunique([M1.x, M1.y, M2.y, M2s.y]) # hide
49
+ M2s # display this # hide
48
50
```
49
51
The derivatives are now with respect to the new independent variable $x$, which can be accessed with ` M2.x ` .
52
+
53
+ !!! warn
54
+ At this point ` x ` , ` M1.x ` , ` M1s.x ` , ` M2.x ` , ` M2s.x ` are * three* different variables.
55
+ Meanwhile ` y ` , ` M1.y ` , ` M1s.y ` , ` M2.y ` and ` M2s.y ` are * four* different variables.
56
+ It can be instructive to inspect these yourself to see their subtle differences.
57
+
50
58
Notice how the number of equations has also decreased from three to two, as $\mathrm{d}x/\mathrm{d}t$ has been turned into an observed equation.
51
59
It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
52
60
``` @example changeivar
53
61
using OrdinaryDiffEq, Plots
54
- prob = ODEProblem(M2s, [], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
62
+ prob = ODEProblem(M2s, [M2s.y => 0.0 ], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
55
63
sol = solve(prob, Tsit5())
56
64
plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
57
65
```
0 commit comments