Rewrite change_independent_variable to handle any derivative order and nonlinear equations

Author: hersle

docs/src/tutorials/change_independent_variable.md

@@ -39,7 +39,7 @@ There are at least three ways of answering this:
 We will demonstrate the last method by changing the independent variable from $t$ to $x$.
 This transformation is well-defined for any non-zero horizontal velocity $v$.
 ```@example changeivar
-M2 = change_independent_variable(M1, x; dummies = true)
+M2 = change_independent_variable(M1, x)
 M2s = structural_simplify(M2; allow_symbolic = true)
 # a sanity test on the 10 x/y variables that are accessible to the user # hide
 @assert allequal([x, M1s.x]) # hide
@@ -110,7 +110,7 @@ To do this, we will change the independent variable in two stages; from $t$ to $
 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.
 First, we transform from $t$ to $a$:
 ```@example changeivar
-M2 = change_independent_variable(M1, M1.a; dummies = true)
+M2 = change_independent_variable(M1, M1.a)
 ```
 Unlike the original, notice that this system is *non-autonomous* because the independent variable $a$ appears explicitly in the equations!
 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$:

src/systems/diffeqs/basic_transformations.jl

@@ -51,22 +51,19 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
 end
 
 """
-    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, add_old_diff = false, simplify = true, fold = false, kwargs...)
+    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_old_diff = false, simplify = true, fold = false, kwargs...)
 
 Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``u(t)``).
-An equation in `sys` must define the rate of change of the new independent variable (e.g. ``du(t)/dt``).
-This or other additional equations can also be specified through `eqs`.
-
 The transformation is well-defined when the mapping between the new and old independent variables are one-to-one.
 This is satisfied if one is a strictly increasing function of the other (e.g. ``du(t)/dt > 0`` or ``du(t)/dt < 0``).
 
+Any extra equations `eqs` involving the new and old independent variables will be taken into account in the transformation.
+
 # Keyword arguments
 
-- `dummies`: Whether derivatives of the new independent variable with respect to the old one are expressed through dummy equations or explicitly inserted into the equations.
-- `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.
+- `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)``.
 - `simplify`: Whether expanded derivative expressions are simplified. This can give a tidier transformation.
 - `fold`: Whether internal substitutions will evaluate numerical expressions.
-Additional keyword arguments `kwargs...` are forwarded to the constructor that rebuilds `sys`.
 
 # Usage before structural simplification
 
@@ -77,6 +74,7 @@ Subsequently, consider passing `allow_symbolic = true` to `structural_simplify(s
 
 If `sys` is non-autonomous (i.e. ``t`` appears explicitly in its equations), it is often desirable to also pass an algebraic equation relating the new and old independent variables (e.g. ``t = f(u(t))``).
 Otherwise the transformed system will be underdetermined and cannot be structurally simplified without additional changes.
+If an algebraic relation is not known, consider using `add_old_diff`.
 
 # Usage with hierarchical systems
 
@@ -91,16 +89,20 @@ By changing the independent variable, it can be reformulated for vertical positi
 ```julia
 julia> @variables x(t) y(t);
 
-julia> @named M = ODESystem([D(D(y)) ~ -9.81, D(x) ~ 10.0], t);
+julia> @named M = ODESystem([D(D(y)) ~ -9.81, D(D(x)) ~ 0.0], t);
+
+julia> M = change_independent_variable(M, x);
 
-julia> M′ = change_independent_variable(complete(M), x);
+julia> M = structural_simplify(M; allow_symbolic = true);
 
-julia> unknowns(M′)
-1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
+julia> unknowns(M)
+3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
+ xˍt(x)
  y(x)
+ yˍx(x)
 ```
 """
-function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, add_old_diff = false, simplify = true, fold = false, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_old_diff = false, simplify = true, fold = false, kwargs...)
     iv2_of_iv1 = unwrap(iv) # e.g. u(t)
     iv1 = get_iv(sys) # e.g. t
 
@@ -115,78 +117,45 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     iv1name = nameof(iv1) # e.g. :t
     iv2name = nameof(operation(iv2_of_iv1)) # e.g. :u
     iv2, = @independent_variables $iv2name # e.g. u
-    iv1_of_iv2, = @variables $iv1name(iv2) # inverse in case sys is autonomous; e.g. t(u)
-    iv1_of_iv2 = GlobalScope(iv1_of_iv2) # do not namespace old independent variable as new dependent variable
+    iv1_of_iv2, = GlobalScope.(@variables $iv1name(iv2)) # inverse, e.g. t(u), global because iv1 has no namespacing in sys
     D1 = Differential(iv1) # e.g. d/d(t)
-    D2 = Differential(iv2_of_iv1) # e.g. d/d(u(t))
+    div2_of_iv1 = GlobalScope(default_toterm(D1(iv2_of_iv1))) # e.g. uˍt(t)
+    div2_of_iv2 = substitute(div2_of_iv1, iv1 => iv2) # e.g. uˍt(u)
+    div2_of_iv2_of_iv1 = substitute(div2_of_iv2, iv2 => iv2_of_iv1) # e.g. uˍt(u(t))
+
+    # If requested, add a differential equation for the old independent variable as a function of the old one
+    if add_old_diff
+        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)
+    end
+    @set! sys.eqs = [get_eqs(sys); eqs] # add extra equations we derived before starting transformation process
+    @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) # add dummy variables and old independent variable as a function of the new one
 
-    # 1) Utility that performs the chain rule on an expression, e.g. (d/dt)(f(t)) -> (d/dt)(f(u(t))) -> df(u(t))/du(t) * du(t)/dt
-    function chain_rule(ex)
+    # Create a utility that performs the chain rule on an expression, followed by insertion of the new independent variable
+    #   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)
+    # Then use it to transform everything in the system!
+    function transform(ex)
+        # 1) Replace the argument of every function; e.g. f(t) -> f(u(t))
         for var in vars(ex; op = Nothing) # loop over all variables in expression (op = Nothing prevents interpreting "D(f(t))" as one big variable)
             is_function_of_iv1 = iscall(var) && isequal(only(arguments(var)), iv1) # is the expression of the form f(t)?
-            if is_function_of_iv1 && !isequal(var, iv2_of_iv1) # substitute f(t) -> f(u(t)), but not u(t) -> u(u(t))
+            if is_function_of_iv1 && !isequal(var, iv2_of_iv1) # prevent e.g. u(t) -> u(u(t))
                 var_of_iv1 = var # e.g. f(t)
-                var_of_iv2 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(u(t))
-                ex = substitute(ex, var_of_iv1 => var_of_iv2)
+                var_of_iv2_of_iv1 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(u(t))
+                ex = substitute(ex, var_of_iv1 => var_of_iv2_of_iv1; fold)
             end
         end
-        ex = expand_derivatives(ex, simplify) # expand chain rule, e.g. (d/dt)(f(u(t)))) -> df(u(t))/du(t) * du(t)/dt
-        return ex
-    end
-
-    # 2) Find e.g. du/dt in equations, then calculate e.g. d²u/dt², ...
-    eqs = [eqs; get_eqs(sys)] # all equations (system-defined + user-provided) we may use
-    idxs = findall(eq -> isequal(eq.lhs, D1(iv2_of_iv1)), eqs)
-    if length(idxs) != 1
-        error("Exactly one equation for $D1($iv2_of_iv1) was not specified! Got $(length(idxs)) equations:\n", join(eqs[idxs], '\n'))
-    end
-    div2_of_iv1_eq = popat!(eqs, only(idxs)) # get and remove e.g. du/dt = ... (may be added back later as a dummy)
-    div2_of_iv1 = chain_rule(div2_of_iv1_eq.rhs)
-    if isequal(div2_of_iv1, 0) # e.g. du/dt ~ 0
-        error("Independent variable transformation $(div2_of_iv1_eq) is singular!")
-    end
-    ddiv2_of_iv1 = chain_rule(D1(div2_of_iv1)) # TODO: implement higher orders (order >= 3) derivatives with a loop
-
-    # 3) If requested, insert extra dummy equations for e.g. du/dt, d²u/dt², ...
-    #    Otherwise, replace all these derivatives by their explicit expressions
-    unks = get_unknowns(sys)
-    if dummies
-        div2 = substitute(default_toterm(D1(iv2_of_iv1)), iv1 => iv2) # e.g. uˍt(u)
-        ddiv2 = substitute(default_toterm(D1(D1(iv2_of_iv1))), iv1 => iv2) # e.g. uˍtt(u)
-        div2, ddiv2 = GlobalScope.([div2, ddiv2]) # do not namespace dummies in new system
-        eqs = [[div2 ~ div2_of_iv1, ddiv2 ~ ddiv2_of_iv1]; eqs] # add dummy equations
-        unks = [unks; [div2, ddiv2]] # add dummy variables
-    else
-        div2 = div2_of_iv1
-        ddiv2 = ddiv2_of_iv1
-    end
-
-    # 4) If requested, add a differential equation for the old independent variable as a function of the old one
-    if add_old_diff
-        div1lhs = substitute(D2(iv1_of_iv2), iv2 => iv2_of_iv1) # e.g. (d/du(t))(t(u(t))); will be transformed to (d/du)(t(u)) by chain rule
-        div1rhs = 1 / div2 # du/dt = 1/(dt/du) (https://en.wikipedia.org/wiki/Inverse_function_rule)
-        eqs = [eqs; div1lhs ~ div1rhs]
-        unks = [unks; iv1_of_iv2]
-    end
-    @set! sys.eqs = eqs # add extra equations we derived before starting transformation process
-    @set! sys.unknowns = unks # add dummy variables and old independent variable as a function of the new one
-
-    # 5) Transform everything from old to new independent variable, e.g. t -> u.
-    #    Substitution order matters! Must begin with highest order to get D(D(u(t))) -> u_tt(u).
-    #    If we had started with the lowest order, we would get D(D(u(t))) -> D(u_t(u)) -> 0!
-    iv1_to_iv2_subs = [ # a vector ensures substitution order
-        D1(D1(iv2_of_iv1)) => ddiv2 # order 2, e.g. D(D(u(t))) -> u_tt(u) or explicit expression
-        D1(iv2_of_iv1) => div2 # order 1, e.g. D(u(t)) -> u_t(u) or explicit expression
-        iv2_of_iv1 => iv2 # order 0, e.g. u(t) -> u
-        iv1 => iv1_of_iv2 # in case sys was autonomous, e.g. t -> t(u)
-    ]
-    function transform(ex)
-        ex = chain_rule(ex)
-        for sub in iv1_to_iv2_subs
-            ex = substitute(ex, sub; fold)
+        # 2) Expand with chain rule until nothing changes anymore
+        orgex = nothing
+        while !isequal(ex, orgex)
+            orgex = ex # save original
+            ex = expand_derivatives(ex, simplify) # expand chain rule, e.g. (d/dt)(f(u(t)))) -> df(u(t))/du(t) * du(t)/dt
+            ex = substitute(ex, D1(iv2_of_iv1) => div2_of_iv2_of_iv1; fold) # e.g. du(t)/dt -> uˍt(u(t))
         end
+        # 3) Set new independent variable
+        ex = substitute(ex, iv2_of_iv1 => iv2; fold) # set e.g. u(t) -> u everywhere
+        ex = substitute(ex, iv1 => iv1_of_iv2; fold) # set e.g. t -> t(u) everywhere
         return ex
     end
+
     function transform(sys::AbstractODESystem)
         eqs = map(transform, get_eqs(sys))
         unknowns = map(transform, get_unknowns(sys))
@@ -203,7 +172,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
         wascomplete = iscomplete(sys) # save before reconstructing system
         sys = typeof(sys)( # recreate system with transformed fields
             eqs, iv2, unknowns, ps; observed, initialization_eqs, parameter_dependencies, defaults, guesses,
-            assertions, name = nameof(sys), description = description(sys), kwargs...
+            assertions, name = nameof(sys), description = description(sys)
         )
         systems = map(transform, systems) # recurse through subsystems
         sys = compose(sys, systems) # rebuild hierarchical system
@@ -213,5 +182,6 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
         end
         return sys
     end
+
     return transform(sys)
 end

test/basic_transformations.jl

@@ -30,8 +30,9 @@ end
     eqs1 = [D(D(x)) ~ D(x) + x, D(y) ~ 1]
     M1 = ODESystem(eqs1, t; name = :M)
     M2 = change_independent_variable(M1, y)
-    eqs2 = substitute(equations(M2), M2.y => M1.t) # system should be equivalent when parametrized with y (since D(y) ~ 1), so substitute back ...
-    @test eqs1[1] == only(eqs2) # ... and check that the equations are unmodified
+    @variables y x(y) yˍt(y)
+    Dy = Differential(y)
+    @test Set(equations(M2)) == Set([yˍt^2*(Dy^2)(x) + yˍt*Dy(yˍt)*Dy(x) ~ x + Dy(x)*yˍt, yˍt ~ 1])
 end
 
 @testset "Change independent variable" begin
@@ -44,7 +45,7 @@ end
     ]
     initialization_eqs = [x ~ 1.0, y ~ 1.0, D(y) ~ 0.0]
     M1 = ODESystem(eqs, t; initialization_eqs, name = :M)
-    M2 = change_independent_variable(M1, s; dummies = true)
+    M2 = change_independent_variable(M1, s)
 
     M1 = structural_simplify(M1; allow_symbolic = true)
     M2 = structural_simplify(M2; allow_symbolic = true)
@@ -77,16 +78,15 @@ end
     M1 = compose(M1, r, m, Λ)
 
     # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
-    M2 = change_independent_variable(M1, M1.a; dummies = true)
+    M2 = change_independent_variable(M1, M1.a)
     M2c = complete(M2) # just for the following equation comparison (without namespacing)
-    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, aˍt, aˍtt = M2c.a, M2c.ȧ, M2c.Ω, M2c.r.Ω, M2c.m.Ω, M2c.Λ.Ω, M2c.ϕ, M2c.aˍt, M2c.aˍtt
+    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, aˍt = M2c.a, M2c.ȧ, M2c.Ω, M2c.r.Ω, M2c.m.Ω, M2c.Λ.Ω, M2c.ϕ, M2c.aˍt
     Da = Differential(a)
     @test Set(equations(M2)) == Set([
-        aˍt ~ ȧ # 1st order dummy equation
-        aˍtt ~ Da(ȧ) * aˍt # 2nd order dummy equation
+        aˍt ~ ȧ # dummy equation
         Ω ~ Ωr + Ωm + ΩΛ
         ȧ ~ √(Ω) * a^2
-        aˍtt*Da(ϕ) + aˍt^2*(Da^2)(ϕ) ~ -3*aˍt^2/a*Da(ϕ)
+        Da(aˍt)*Da(ϕ)*aˍt + aˍt^2*(Da^2)(ϕ) ~ -3*aˍt^2/a*Da(ϕ)
         aˍt*Da(Ωr) ~ -4*Ωr*aˍt/a
         aˍt*Da(Ωm) ~ -3*Ωm*aˍt/a
         aˍt*Da(ΩΛ) ~ 0
@@ -104,13 +104,13 @@ end
 @testset "Change independent variable (simple)" begin
     @variables x(t) y1(t) # y(t)[1:1] # TODO: use array variables y(t)[1:2] when fixed: https://github.com/JuliaSymbolics/Symbolics.jl/issues/1383
     Mt = ODESystem([D(x) ~ 2*x, D(y1) ~ y1], t; name = :M)
-    Mx = change_independent_variable(Mt, x; dummies = true)
+    Mx = change_independent_variable(Mt, x)
     @variables x xˍt(x) xˍtt(x) y1(x) # y(x)[1:1] # TODO: array variables
     Dx = Differential(x)
-    @test (Set(equations(Mx)) == Set([xˍt ~ 2*x, xˍtt ~ 2*xˍt, xˍt*Dx(y1) ~ y1]))
+    @test (Set(equations(Mx)) == Set([xˍt ~ 2*x, xˍt*Dx(y1) ~ y1]))
 end
 
-@testset "Change independent variable (free fall)" begin
+@testset "Change independent variable (free fall with 1st order horizontal equation)" begin
     @variables x(t) y(t)
     @parameters g v # gravitational acceleration and constant horizontal velocity
     Mt = ODESystem([D(D(y)) ~ -g, D(x) ~ v], t; name = :M) # gives (x, y) as function of t, ...
@@ -122,6 +122,18 @@ end
     @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.t)^2/2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2)
 end
 
+@testset "Change independent variable (free fall with 2nd order horizontal equation)" begin
+    @variables x(t) y(t)
+    @parameters g # gravitational acceleration
+    Mt = ODESystem([D(D(y)) ~ -g, D(D(x)) ~ 0], t; name = :M) # gives (x, y) as function of t, ...
+    Mx = change_independent_variable(Mt, x; add_old_diff = true) # ... but we want y as a function of x
+    Mx = structural_simplify(Mx; allow_symbolic = true)
+    Dx = Differential(Mx.x)
+    prob = ODEProblem(Mx, [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0, Mx.xˍt => 10.0], (0.0, 20.0), [g => 9.81]) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities
+    sol = solve(prob, Tsit5(); reltol = 1e-5)
+    @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.t)^2/2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2)
+end
+
 @testset "Change independent variable (crazy analytical example)" begin
     @independent_variables t
     D = Differential(t)
@@ -130,16 +142,15 @@ end
         D(D(y)) ~ D(x)^2 + D(y^3) |> expand_derivatives # expand D(y^3) # TODO: make this test 3rd order
         D(x) ~ x^4 + y^5 + t^6
     ], t; name = :M)
-    M2 = change_independent_variable(M1, x; dummies = true)
+    M2 = change_independent_variable(M1, x)
 
     # Compare to pen-and-paper result
     @independent_variables x
     Dx = Differential(x)
     @variables xˍt(x) xˍtt(x) y(x) t(x)
     @test Set(equations(M2)) == Set([
-        xˍt^2*(Dx^2)(y) + xˍtt*Dx(y) ~ xˍt^2 + 3*y^2*Dx(y)*xˍt # from D(D(y))
-        xˍt ~ x^4 + y^5 + t^6 # 1st order dummy equation
-        xˍtt ~ 4*x^3*xˍt + 5*y^4*Dx(y)*xˍt + 6*t^5 # 2nd order dummy equation
+        xˍt^2*(Dx^2)(y) + xˍt*Dx(xˍt)*Dx(y) ~ xˍt^2 + 3*y^2*Dx(y)*xˍt # from D(D(y))
+        xˍt ~ x^4 + y^5 + t^6 # dummy equation
     ])
 end
 
@@ -157,11 +168,12 @@ end
 
     # Ensure that interpolations are called with the same variables
     M2 = change_independent_variable(M1, x, [t ~ √(x)])
-    @variables x y(x) t(x)
+    @variables x xˍt(x) y(x) t(x)
     Dx = Differential(x)
     @test Set(equations(M2)) == Set([
         t ~ √(x)
-        2*t*Dx(y) ~ 1*fc(t) + 2*fc(x) + 3*fc(y) + 1*callme(f, t) + 2*callme(f, x) + 3*callme(f, y)
+        xˍt ~ 2*t
+        xˍt*Dx(y) ~ 1*fc(t) + 2*fc(x) + 3*fc(y) + 1*callme(f, t) + 2*callme(f, x) + 3*callme(f, y)
     ])
 
     _f = LinearInterpolation([1.0, 1.0], [-100.0, +100.0]) # constant value 1
@@ -173,18 +185,10 @@ end
 
 @testset "Change independent variable (errors)" begin
     @variables x(t) y z(y) w(t) v(t)
-    M = ODESystem([D(x) ~ 0, v ~ x], t; name = :M)
-    @test_throws "singular" change_independent_variable(M, x)
+    M = ODESystem([D(x) ~ 1, v ~ x], t; name = :M)
     @test_throws "structurally simplified" change_independent_variable(structural_simplify(M), y)
-    @test_throws "Got 0 equations:" change_independent_variable(M, w)
-    @test_throws "Got 0 equations:" change_independent_variable(M, v)
-    M = ODESystem([2 * D(x) ~ 1, v ~ x], t; name = :M) # TODO: allow equations like this
-    @test_throws "Got 0 equations:" change_independent_variable(M, x)
-    M = ODESystem([D(x) ~ 1, v ~ 1], t; name = :M)
-    @test_throws "Got 2 equations:" change_independent_variable(M, x, [D(x) ~ 2])
     @test_throws "not a function of the independent variable" change_independent_variable(M, y)
     @test_throws "not a function of the independent variable" change_independent_variable(M, z)
-    M = ODESystem([D(x) ~ 0, v ~ x], t; name = :M)
     @variables x(..) # require explicit argument
     M = ODESystem([D(x(t)) ~ x(t-1)], t; name = :M)
     @test_throws "DDE" change_independent_variable(M, x(t))