Skip to content

Commit 2b9ec60

Browse files
Merge pull request #3437 from hersle/change_ivar
Change independent variable of ODE systems
2 parents 87872dd + 7f821e2 commit 2b9ec60

File tree

7 files changed

+513
-39
lines changed

7 files changed

+513
-39
lines changed

docs/pages.jl

+1
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ pages = [
1010
"tutorials/stochastic_diffeq.md",
1111
"tutorials/discrete_system.md",
1212
"tutorials/parameter_identifiability.md",
13+
"tutorials/change_independent_variable.md",
1314
"tutorials/bifurcation_diagram_computation.md",
1415
"tutorials/SampledData.md",
1516
"tutorials/domain_connections.md",

docs/src/systems/ODESystem.md

+1
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@ ODESystem
3030
structural_simplify
3131
ode_order_lowering
3232
dae_index_lowering
33+
change_independent_variable
3334
liouville_transform
3435
alias_elimination
3536
tearing
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,146 @@
1+
# Changing the independent variable of ODEs
2+
3+
Ordinary differential equations describe the rate of change of some dependent variables with respect to one independent variable.
4+
For the modeler it is often most natural to write down the equations with a particular independent variable, say time $t$.
5+
However, in many cases there are good reasons for changing the independent variable:
6+
7+
1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
8+
2. Some differential equations vary more nicely (e.g. less stiff) with respect to one independent variable than another.
9+
3. It can reduce the number of equations that must be solved (e.g. $y(x)$ is one equation, while $(x(t), y(t))$ are two).
10+
11+
To manually change the independent variable of an ODE, one must rewrite all equations in terms of a new variable and transform differentials with the chain rule.
12+
This is mechanical and error-prone.
13+
ModelingToolkit provides the utility function [`change_independent_variable`](@ref) that automates this process.
14+
15+
## 1. Get one dependent variable as a function of another
16+
17+
Consider a projectile shot with some initial velocity in a vertical gravitational field with constant horizontal velocity.
18+
19+
```@example changeivar
20+
using ModelingToolkit
21+
using ModelingToolkit: t_nounits as t, D_nounits as D
22+
@variables x(t) y(t)
23+
@parameters g=9.81 v # gravitational acceleration and horizontal velocity
24+
eqs = [D(D(y)) ~ -g, D(x) ~ v]
25+
initialization_eqs = [D(x) ~ D(y)] # 45° initial angle
26+
M1 = ODESystem(eqs, t; initialization_eqs, name = :M)
27+
M1s = structural_simplify(M1)
28+
@assert length(equations(M1s)) == 3 # hide
29+
M1s # hide
30+
```
31+
32+
This is the standard parametrization that arises naturally from kinematics and Newton's laws.
33+
It expresses the position $(x(t), y(t))$ as a function of time $t$.
34+
But suppose we want to determine whether the projectile hits a target 10 meters away.
35+
There are at least three ways of answering this:
36+
37+
- Solve the ODE for $(x(t), y(t))$ and use a callback to terminate when $x$ reaches 10 meters, and evaluate $y$ at the final time.
38+
- Solve the ODE for $(x(t), y(t))$ and use root finding to find the time when $x$ reaches 10 meters, and evaluate $y$ at that time.
39+
- Solve the ODE for $y(x)$ and evaluate it at 10 meters.
40+
41+
We will demonstrate the last method by changing the independent variable from $t$ to $x$.
42+
This transformation is well-defined for any non-zero horizontal velocity $v$, so $x$ and $t$ are one-to-one.
43+
44+
```@example changeivar
45+
M2 = change_independent_variable(M1, x)
46+
M2s = structural_simplify(M2; allow_symbolic = true)
47+
# a sanity test on the 10 x/y variables that are accessible to the user # hide
48+
@assert allequal([x, M1s.x]) # hide
49+
@assert allequal([M2.x, M2s.x]) # hide
50+
@assert allequal([y, M1s.y]) # hide
51+
@assert allunique([M1.x, M1.y, M2.y, M2s.y]) # hide
52+
@assert length(equations(M2s)) == 2 # hide
53+
M2s # display this # hide
54+
```
55+
56+
The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`.
57+
58+
!!! warn
59+
60+
At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables.
61+
Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables.
62+
It can be instructive to inspect these yourself to see their subtle differences.
63+
64+
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.
65+
It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
66+
67+
```@example changeivar
68+
using OrdinaryDiffEq, Plots
69+
prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
70+
sol = solve(prob, Tsit5())
71+
plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
72+
```
73+
74+
!!! tip "Usage tips"
75+
76+
Look up the documentation of [`change_independent_variable`](@ref) for tips on how to use it.
77+
78+
For example, if you also need $t(x)$, you can tell it to add a differential equation for the old independent variable in terms of the new one using the [inverse function rule](https://en.wikipedia.org/wiki/Inverse_function_rule) (here $\mathrm{d}t/\mathrm{d}x = 1 / (\mathrm{d}x/\mathrm{d}t)$). If you know an analytical expression between the independent variables (here $t = x/v$), you can also pass it directly to the function to avoid the extra differential equation.
79+
80+
## 2. Alleviating stiffness by changing to logarithmic time
81+
82+
In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe.
83+
In terms of conformal time $t$, they can be written
84+
85+
```@example changeivar
86+
@variables a(t) Ω(t)
87+
a = GlobalScope(a) # global var needed by all species
88+
function species(w; kw...)
89+
eqs = [D(Ω) ~ -3(1 + w) * D(a) / a * Ω]
90+
return ODESystem(eqs, t, [Ω], []; kw...)
91+
end
92+
@named r = species(1 // 3) # radiation
93+
@named m = species(0) # matter
94+
@named Λ = species(-1) # dark energy / cosmological constant
95+
eqs = [Ω ~ r.Ω + m.Ω + Λ.Ω, D(a) ~ √(Ω) * a^2]
96+
initialization_eqs = [Λ.Ω + r.Ω + m.Ω ~ 1]
97+
M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M)
98+
M1 = compose(M1, r, m, Λ)
99+
M1s = structural_simplify(M1)
100+
```
101+
102+
Of course, we can solve this ODE as it is:
103+
104+
```@example changeivar
105+
prob = ODEProblem(M1s, [M1s.a => 1.0, M1s.r.Ω => 5e-5, M1s.m.Ω => 0.3], (0.0, -3.3), [])
106+
sol = solve(prob, Tsit5(); reltol = 1e-5)
107+
@assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide
108+
plot(sol, idxs = [M1.a, M1.r.Ω / M1.Ω, M1.m.Ω / M1.Ω, M1.Λ.Ω / M1.Ω])
109+
```
110+
111+
But the solver becomes unstable due to stiffness.
112+
Also notice the interesting dynamics taking place towards the end of the integration (in the early universe), which gets compressed into a very small time interval.
113+
These ODEs would benefit from being defined with respect to a logarithmic "time" that better captures the evolution of the universe through *orders of magnitude* of time, rather than linear time.
114+
115+
It is therefore common to write these ODEs in terms of $b = \ln a$.
116+
To do this, we will change the independent variable in two stages; first from $t$ to $a$, and then from $a$ to $b$.
117+
Notice that $\mathrm{d}a/\mathrm{d}t > 0$ provided that $\Omega > 0$, and $\mathrm{d}b/\mathrm{d}a > 0$, so the transformation is well-defined since $t \leftrightarrow a \leftrightarrow b$ are one-to-one.
118+
First, we transform from $t$ to $a$:
119+
120+
```@example changeivar
121+
M2 = change_independent_variable(M1, M1.a)
122+
@assert !ModelingToolkit.isautonomous(M2) # hide
123+
M2 # hide
124+
```
125+
126+
Unlike the original, notice that this system is *non-autonomous* because the independent variable $a$ appears explicitly in the equations!
127+
This means that to change the independent variable from $a$ to $b$, we must provide not only the rate of change relation $db(a)/da = \exp(-b)$, but *also* the equation $a(b) = \exp(b)$ so $a$ can be eliminated in favor of $b$:
128+
129+
```@example changeivar
130+
a = M2.a # get independent variable of M2
131+
Da = Differential(a)
132+
@variables b(a)
133+
M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)])
134+
```
135+
136+
We can now solve and plot the ODE in terms of $b$:
137+
138+
```@example changeivar
139+
M3s = structural_simplify(M3; allow_symbolic = true)
140+
prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), [])
141+
sol = solve(prob, Tsit5())
142+
@assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide
143+
plot(sol, idxs = [M3.r.Ω / M3.Ω, M3.m.Ω / M3.Ω, M3.Λ.Ω / M3.Ω])
144+
```
145+
146+
Notice that the variables vary "more nicely" with respect to $b$ than $t$, making the solver happier and avoiding numerical problems.

src/ModelingToolkit.jl

+2-1
Original file line numberDiff line numberDiff line change
@@ -259,7 +259,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
259259
tunable_parameters, isirreducible, getdescription, hasdescription,
260260
hasunit, getunit, hasconnect, getconnect,
261261
hasmisc, getmisc, state_priority
262-
export ode_order_lowering, dae_order_lowering, liouville_transform
262+
export ode_order_lowering, dae_order_lowering, liouville_transform,
263+
change_independent_variable
263264
export PDESystem
264265
export Differential, expand_derivatives, @derivatives
265266
export Equation, ConstrainedEquation

src/systems/diffeqs/basic_transformations.jl

+156-15
Original file line numberDiff line numberDiff line change
@@ -14,26 +14,20 @@ the final value of `trJ` is the probability of ``u(t)``.
1414
Example:
1515
1616
```julia
17-
using ModelingToolkit, OrdinaryDiffEq, Test
17+
using ModelingToolkit, OrdinaryDiffEq
1818
1919
@independent_variables t
2020
@parameters α β γ δ
2121
@variables x(t) y(t)
2222
D = Differential(t)
23+
eqs = [D(x) ~ α*x - β*x*y, D(y) ~ -δ*y + γ*x*y]
24+
@named sys = ODESystem(eqs, t)
2325
24-
eqs = [D(x) ~ α*x - β*x*y,
25-
D(y) ~ -δ*y + γ*x*y]
26-
27-
sys = ODESystem(eqs)
2826
sys2 = liouville_transform(sys)
29-
@variables trJ
30-
31-
u0 = [x => 1.0,
32-
y => 1.0,
33-
trJ => 1.0]
34-
35-
prob = ODEProblem(complete(sys2),u0,tspan,p)
36-
sol = solve(prob,Tsit5())
27+
sys2 = complete(sys2)
28+
u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0]
29+
prob = ODEProblem(sys2, u0, tspan, p)
30+
sol = solve(prob, Tsit5())
3731
```
3832
3933
Where `sol[3,:]` is the evolution of `trJ` over time.
@@ -46,12 +40,159 @@ Optimal Transport Approach
4640
Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya
4741
https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf
4842
"""
49-
function liouville_transform(sys::AbstractODESystem)
43+
function liouville_transform(sys::AbstractODESystem; kwargs...)
5044
t = get_iv(sys)
5145
@variables trJ
5246
D = ModelingToolkit.Differential(t)
5347
neweq = D(trJ) ~ trJ * -tr(calculate_jacobian(sys))
5448
neweqs = [equations(sys); neweq]
5549
vars = [unknowns(sys); trJ]
56-
ODESystem(neweqs, t, vars, parameters(sys), checks = false)
50+
ODESystem(
51+
neweqs, t, vars, parameters(sys);
52+
checks = false, name = nameof(sys), kwargs...
53+
)
54+
end
55+
56+
"""
57+
change_independent_variable(
58+
sys::AbstractODESystem, iv, eqs = [];
59+
add_old_diff = false, simplify = true, fold = false
60+
)
61+
62+
Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``u(t)``).
63+
The transformation is well-defined when the mapping between the new and old independent variables are one-to-one.
64+
This is satisfied if one is a strictly increasing function of the other (e.g. ``du(t)/dt > 0`` or ``du(t)/dt < 0``).
65+
66+
Any extra equations `eqs` involving the new and old independent variables will be taken into account in the transformation.
67+
68+
# Keyword arguments
69+
70+
- `add_old_diff`: Whether to add a differential equation for the old independent variable in terms of the new one using the inverse function rule ``dt/du = 1/(du/dt)``.
71+
- `simplify`: Whether expanded derivative expressions are simplified. This can give a tidier transformation.
72+
- `fold`: Whether internal substitutions will evaluate numerical expressions.
73+
74+
# Usage before structural simplification
75+
76+
The variable change must take place before structural simplification.
77+
In following calls to `structural_simplify`, consider passing `allow_symbolic = true` to avoid undesired constraint equations between between dummy variables.
78+
79+
# Usage with non-autonomous systems
80+
81+
If `sys` is non-autonomous (i.e. ``t`` appears explicitly in its equations), consider passing an algebraic equation relating the new and old independent variables (e.g. ``t = f(u(t))``).
82+
Otherwise the transformed system can be underdetermined.
83+
If an algebraic relation is not known, consider using `add_old_diff` instead.
84+
85+
# Usage with hierarchical systems
86+
87+
It is recommended that `iv` is a non-namespaced variable in `sys`.
88+
This means it can belong to the top-level system or be a variable in a subsystem declared with `GlobalScope`.
89+
90+
# Example
91+
92+
Consider a free fall with constant horizontal velocity.
93+
Physics naturally describes position as a function of time.
94+
By changing the independent variable, it can be reformulated for vertical position as a function of horizontal position:
95+
```julia
96+
julia> @variables x(t) y(t);
97+
98+
julia> @named M = ODESystem([D(D(y)) ~ -9.81, D(D(x)) ~ 0.0], t);
99+
100+
julia> M = change_independent_variable(M, x);
101+
102+
julia> M = structural_simplify(M; allow_symbolic = true);
103+
104+
julia> unknowns(M)
105+
3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
106+
xˍt(x)
107+
y(x)
108+
yˍx(x)
109+
```
110+
"""
111+
function change_independent_variable(
112+
sys::AbstractODESystem, iv, eqs = [];
113+
add_old_diff = false, simplify = true, fold = false
114+
)
115+
iv2_of_iv1 = unwrap(iv) # e.g. u(t)
116+
iv1 = get_iv(sys) # e.g. t
117+
118+
if is_dde(sys)
119+
error("System $(nameof(sys)) contains delay differential equations (DDEs). This is currently not supported!")
120+
elseif isscheduled(sys)
121+
error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!")
122+
elseif !iscall(iv2_of_iv1) || !isequal(only(arguments(iv2_of_iv1)), iv1)
123+
error("Variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))!")
124+
end
125+
126+
# Set up intermediate and final variables for the transformation
127+
iv1name = nameof(iv1) # e.g. :t
128+
iv2name = nameof(operation(iv2_of_iv1)) # e.g. :u
129+
iv2, = @independent_variables $iv2name # e.g. u
130+
iv1_of_iv2, = GlobalScope.(@variables $iv1name(iv2)) # inverse, e.g. t(u), global because iv1 has no namespacing in sys
131+
D1 = Differential(iv1) # e.g. d/d(t)
132+
div2_of_iv1 = GlobalScope(default_toterm(D1(iv2_of_iv1))) # e.g. uˍt(t)
133+
div2_of_iv2 = substitute(div2_of_iv1, iv1 => iv2) # e.g. uˍt(u)
134+
div2_of_iv2_of_iv1 = substitute(div2_of_iv2, iv2 => iv2_of_iv1) # e.g. uˍt(u(t))
135+
136+
# If requested, add a differential equation for the old independent variable as a function of the old one
137+
if add_old_diff
138+
eqs = [eqs; Differential(iv2)(iv1_of_iv2) ~ 1 / div2_of_iv2] # e.g. dt(u)/du ~ 1 / uˍt(u) (https://en.wikipedia.org/wiki/Inverse_function_rule)
139+
end
140+
@set! sys.eqs = [get_eqs(sys); eqs] # add extra equations we derived
141+
@set! sys.unknowns = [get_unknowns(sys); [iv1, div2_of_iv1]] # add new variables, will be transformed to e.g. t(u) and uˍt(u)
142+
143+
# Create a utility that performs the chain rule on an expression, followed by insertion of the new independent variable:
144+
# e.g. (d/dt)(f(t)) -> (d/dt)(f(u(t))) -> df(u(t))/du(t) * du(t)/dt -> df(u)/du * uˍt(u)
145+
function transform(ex)
146+
# 1) Replace the argument of every function; e.g. f(t) -> f(u(t))
147+
for var in vars(ex; op = Nothing) # loop over all variables in expression (op = Nothing prevents interpreting "D(f(t))" as one big variable)
148+
is_function_of_iv1 = iscall(var) && isequal(only(arguments(var)), iv1) # of the form f(t)?
149+
if is_function_of_iv1 && !isequal(var, iv2_of_iv1) # prevent e.g. u(t) -> u(u(t))
150+
var_of_iv1 = var # e.g. f(t)
151+
var_of_iv2_of_iv1 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(u(t))
152+
ex = substitute(ex, var_of_iv1 => var_of_iv2_of_iv1; fold)
153+
end
154+
end
155+
# 2) Repeatedly expand chain rule until nothing changes anymore
156+
orgex = nothing
157+
while !isequal(ex, orgex)
158+
orgex = ex # save original
159+
ex = expand_derivatives(ex, simplify) # expand chain rule, e.g. (d/dt)(f(u(t)))) -> df(u(t))/du(t) * du(t)/dt
160+
ex = substitute(ex, D1(iv2_of_iv1) => div2_of_iv2_of_iv1; fold) # e.g. du(t)/dt -> uˍt(u(t))
161+
end
162+
# 3) Set new independent variable
163+
ex = substitute(ex, iv2_of_iv1 => iv2; fold) # set e.g. u(t) -> u everywhere
164+
ex = substitute(ex, iv1 => iv1_of_iv2; fold) # set e.g. t -> t(u) everywhere
165+
return ex
166+
end
167+
168+
# Use the utility function to transform everything in the system!
169+
function transform(sys::AbstractODESystem)
170+
eqs = map(transform, get_eqs(sys))
171+
unknowns = map(transform, get_unknowns(sys))
172+
unknowns = filter(var -> !isequal(var, iv2), unknowns) # remove e.g. u
173+
ps = map(transform, get_ps(sys))
174+
ps = filter(!isinitial, ps) # remove Initial(...) # TODO: shouldn't have to touch this
175+
observed = map(transform, get_observed(sys))
176+
initialization_eqs = map(transform, get_initialization_eqs(sys))
177+
parameter_dependencies = map(transform, get_parameter_dependencies(sys))
178+
defaults = Dict(transform(var) => transform(val)
179+
for (var, val) in get_defaults(sys))
180+
guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
181+
assertions = Dict(transform(ass) => msg for (ass, msg) in get_assertions(sys))
182+
systems = get_systems(sys) # save before reconstructing system
183+
wascomplete = iscomplete(sys) # save before reconstructing system
184+
sys = typeof(sys)( # recreate system with transformed fields
185+
eqs, iv2, unknowns, ps; observed, initialization_eqs,
186+
parameter_dependencies, defaults, guesses,
187+
assertions, name = nameof(sys), description = description(sys)
188+
)
189+
systems = map(transform, systems) # recurse through subsystems
190+
sys = compose(sys, systems) # rebuild hierarchical system
191+
if wascomplete
192+
wasflat = isempty(systems)
193+
sys = complete(sys; flatten = wasflat) # complete output if input was complete
194+
end
195+
return sys
196+
end
197+
return transform(sys)
57198
end

0 commit comments

Comments
 (0)