@@ -23,11 +23,8 @@ D = Differential(t)
2323@parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
2424M1 = ODESystem([
2525 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°)
3128], name = :M)
3229M1s = structural_simplify(M1)
3330```
@@ -43,15 +40,26 @@ We will demonstrate the last method by changing the independent variable from $t
4340This transformation is well-defined for any non-zero horizontal velocity $v$.
4441``` @example changeivar
4542M2 = change_independent_variable(M1, x; dummies = true)
46- @assert M2.x isa Num # hide
4743M2s = 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
4850```
4951The 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+
5058Notice 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.
5159It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
5260``` @example changeivar
5361using 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
5563sol = solve(prob, Tsit5())
5664plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
5765```
0 commit comments