From acec3cb7e2ae3037d239a4e09ef57d73e6a72559 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 23 Feb 2025 17:35:10 +0100
Subject: [PATCH 01/32] Change independent variable of simple 1st order system

---
 src/systems/diffeqs/basic_transformations.jl | 31 ++++++++++++++++++++
 test/basic_transformations.jl                | 25 ++++++++++++++++
 2 files changed, 56 insertions(+)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index e2be889c5e..7bf348911b 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -55,3 +55,34 @@ function liouville_transform(sys::AbstractODESystem)
     vars = [unknowns(sys); trJ]
     ODESystem(neweqs, t, vars, parameters(sys), checks = false)
 end
+
+# TODO: handle case when new iv is a variable already in the system
+function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2; kwargs...)
+    iv1 = ModelingToolkit.get_iv(sys) # old independent variable
+    iv2 = iv # new independent variable
+
+    name2 = nameof(iv2)
+    iv2func, = @variables $name2(iv1)
+
+    eqs = ModelingToolkit.get_eqs(sys) |> copy # don't modify original system
+    vars = []
+    for (i, eq) in enumerate(eqs)
+        vars = Symbolics.get_variables(eq)
+        for var1 in vars
+            if Symbolics.iscall(var1) # skip e.g. constants
+                name = nameof(operation(var1))
+                var2, = @variables $name(iv2func)
+                eq = substitute(eq, var1 => var2; fold = false)
+            end
+        end
+        eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
+        div1_div2 = Differential(iv2)(iv1_of_iv2) |> expand_derivatives
+        div2_div1 = 1 / div1_div2
+        eq = substitute(eq, Differential(iv1)(iv2func) => div2_div1) # substitute in d(iv1)/d(iv2)
+        eq = substitute(eq, iv2func => iv2) # make iv2 independent
+        eqs[i] = eq
+    end
+
+    sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
+    return sys2
+end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index d9b59408e6..512510011c 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -32,3 +32,28 @@ u0 = [x => 1.0,
 prob = ODEProblem(sys2, u0, tspan, p, jac = true)
 sol = solve(prob, Tsit5())
 @test sol[end, end] ≈ 1.0742818931017244
+
+@testset "Change independent variable" begin
+    # Autonomous 1st order (t doesn't appear in the equations)
+    @independent_variables t
+    @variables x(t) y(t) z(t)
+    eqs = [
+        D(x) ~ y
+        D(y) ~ 2*D(x)
+        z ~ x + D(y)
+    ]
+    @named sys1 = ODESystem(eqs, t)
+
+    @independent_variables s
+    @named sys2 = ModelingToolkit.change_independent_variable(sys1, s, s^2)
+
+    sys1 = structural_simplify(sys1)
+    sys2 = structural_simplify(sys2)
+    prob1 = ODEProblem(sys1, unknowns(sys1) .=> 1.0, (0.0, 1.0))
+    prob2 = ODEProblem(sys2, unknowns(sys2) .=> 1.0, (0.0, 1.0))
+    sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
+    sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
+    ts = range(0.0, 1.0, length = 50)
+    ss = .√(ts)
+    @test isapprox(sol1(ts), sol2(ss); atol = 1e-8)
+end

From 3600b0de65323f32ee6f47b422f810fc15ed1dc8 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 23 Feb 2025 18:00:12 +0100
Subject: [PATCH 02/32] Change independent variable of 2nd-order systems

---
 src/systems/diffeqs/basic_transformations.jl |  8 ++++----
 test/basic_transformations.jl                | 11 ++++++-----
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 7bf348911b..60eb46483e 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -57,7 +57,7 @@ function liouville_transform(sys::AbstractODESystem)
 end
 
 # TODO: handle case when new iv is a variable already in the system
-function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2; kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2, iv2_of_iv1; kwargs...)
     iv1 = ModelingToolkit.get_iv(sys) # old independent variable
     iv2 = iv # new independent variable
 
@@ -76,10 +76,10 @@ function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2; kwa
             end
         end
         eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
-        div1_div2 = Differential(iv2)(iv1_of_iv2) |> expand_derivatives
-        div2_div1 = 1 / div1_div2
-        eq = substitute(eq, Differential(iv1)(iv2func) => div2_div1) # substitute in d(iv1)/d(iv2)
+        eq = substitute(eq, Differential(iv1)(iv2func) => Differential(iv1)(iv2_of_iv1)) # substitute in d(iv2)/d(iv1)
+        eq = expand_derivatives(eq)
         eq = substitute(eq, iv2func => iv2) # make iv2 independent
+        eq = substitute(eq, iv1 => iv1_of_iv2) # substitute any remaining old ivars
         eqs[i] = eq
     end
 
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 512510011c..efa799a0ee 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -39,21 +39,22 @@ sol = solve(prob, Tsit5())
     @variables x(t) y(t) z(t)
     eqs = [
         D(x) ~ y
-        D(y) ~ 2*D(x)
+        D(D(y)) ~ 2 * x * D(y)
         z ~ x + D(y)
     ]
     @named sys1 = ODESystem(eqs, t)
 
     @independent_variables s
-    @named sys2 = ModelingToolkit.change_independent_variable(sys1, s, s^2)
+    @named sys2 = ModelingToolkit.change_independent_variable(sys1, s, s^2, √(t))
 
     sys1 = structural_simplify(sys1)
     sys2 = structural_simplify(sys2)
-    prob1 = ODEProblem(sys1, unknowns(sys1) .=> 1.0, (0.0, 1.0))
-    prob2 = ODEProblem(sys2, unknowns(sys2) .=> 1.0, (0.0, 1.0))
+    prob1 = ODEProblem(sys1, [sys1.x => 1.0, sys1.y => 1.0, Differential(t)(sys1.y) => 0.0], (1.0, 4.0))
+    prob2 = ODEProblem(sys2, [sys2.x => 1.0, sys2.y => 1.0, Differential(s)(sys2.y) => 0.0], (1.0, 2.0))
     sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     ts = range(0.0, 1.0, length = 50)
     ss = .√(ts)
-    @test isapprox(sol1(ts), sol2(ss); atol = 1e-8)
+    @test all(isapprox.(sol1(ts, idxs=sys1.x), sol2(ss, idxs=sys2.x); atol = 1e-7)) &&
+          all(isapprox.(sol1(ts, idxs=sys1.y), sol2(ss, idxs=sys2.y); atol = 1e-7))
 end

From acd0fabc7aea69e7702e98623f944a7171b85138 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sat, 1 Mar 2025 13:28:27 +0100
Subject: [PATCH 03/32] Change independent variable to a dependent one

---
 src/systems/diffeqs/basic_transformations.jl | 38 ++++++++++++++++++++
 test/basic_transformations.jl                | 21 +++++++++++
 2 files changed, 59 insertions(+)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 60eb46483e..01a103adab 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -86,3 +86,41 @@ function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2, iv2
     sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
     return sys2
 end
+
+function change_independent_variable(sys::AbstractODESystem, iv; kwargs...)
+    iv1 = get_iv(sys) # e.g. t
+    iv2func = iv # e.g. a(t)
+    iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
+    iv2, = @independent_variables $iv2name
+
+    # TODO: find iv2func in system and replace it with some dummy variable
+    # TODO: not just to 1st, but to all orders
+
+    div2name = Symbol(iv2name, :_t)
+    div2, = @variables $div2name(iv2)
+
+    eqs = ModelingToolkit.get_eqs(sys) |> copy # don't modify original system
+    vars = []
+    for (i, eq) in enumerate(eqs)
+        vars = Symbolics.get_variables(eq)
+        for var1 in vars
+            if Symbolics.iscall(var1) && !isequal(var1, iv2func) # && isequal(only(arguments(var1)), iv1) # skip e.g. constants
+                name = nameof(operation(var1))
+                var2, = @variables $name(iv2func)
+                eq = substitute(eq, var1 => var2; fold = false)
+            end
+        end
+        eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
+        #eq = substitute(eq, Differential(iv1)(iv2func) => Differential(iv1)(iv2_of_iv1)) # substitute in d(iv2)/d(iv1)
+        #eq = expand_derivatives(eq)
+        eq = substitute(eq, Differential(iv1)(Differential(iv1)(iv2func)) => div2) # e.g. D(a(t)) => a_t(t) # TODO: more orders
+        eq = substitute(eq, Differential(iv1)(iv2func) => div2) # e.g. D(a(t)) => a_t(t) # TODO: more orders
+        eq = substitute(eq, iv2func => iv2) # make iv2 independent
+        #eq = substitute(eq, iv1 => iv1_of_iv2) # substitute any remaining old ivars
+        eqs[i] = eq
+        println(eq)
+    end
+
+    sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
+    return sys2
+end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index efa799a0ee..3c7a4e6755 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -58,3 +58,24 @@ sol = solve(prob, Tsit5())
     @test all(isapprox.(sol1(ts, idxs=sys1.x), sol2(ss, idxs=sys2.x); atol = 1e-7)) &&
           all(isapprox.(sol1(ts, idxs=sys1.y), sol2(ss, idxs=sys2.y); atol = 1e-7))
 end
+
+#@testset "Change independent variable (Friedmann equation)" begin
+    @independent_variables t
+    D = Differential(t)
+    @variables a(t) ρr(t) ρm(t) ρΛ(t) ρ(t) ϕ(t)
+    @parameters Ωr0 Ωm0 ΩΛ0
+    eqs = [
+        ρr ~ 3/(8*Num(π)) * Ωr0 / a^4
+        ρm ~ 3/(8*Num(π)) * Ωm0 / a^3
+        ρΛ ~ 3/(8*Num(π)) * ΩΛ0
+        ρ ~ ρr + ρm + ρΛ
+        D(a) ~ √(8*Num(π)/3*ρ*a^4)
+        D(D(ϕ)) ~ -3*D(a)/a*D(ϕ)
+    ]
+    @named M1 = ODESystem(eqs, t)
+    M1 = complete(M1)
+
+    M2 = ModelingToolkit.change_independent_variable(M1, M1.a)
+
+    M2 = structural_simplify(M2; allow_symbolic = true)
+#end

From 865ef02623a9fcdb81b6227fa64a7aa4800a3744 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sat, 1 Mar 2025 23:45:22 +0100
Subject: [PATCH 04/32] Merge functions for changing independent variable into
 one

---
 src/systems/diffeqs/basic_transformations.jl | 83 +++++++++-----------
 test/basic_transformations.jl                | 24 +++---
 2 files changed, 50 insertions(+), 57 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 01a103adab..841c0bda6d 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -56,52 +56,28 @@ function liouville_transform(sys::AbstractODESystem)
     ODESystem(neweqs, t, vars, parameters(sys), checks = false)
 end
 
-# TODO: handle case when new iv is a variable already in the system
-function change_independent_variable(sys::AbstractODESystem, iv, iv1_of_iv2, iv2_of_iv1; kwargs...)
-    iv1 = ModelingToolkit.get_iv(sys) # old independent variable
-    iv2 = iv # new independent variable
-
-    name2 = nameof(iv2)
-    iv2func, = @variables $name2(iv1)
-
-    eqs = ModelingToolkit.get_eqs(sys) |> copy # don't modify original system
-    vars = []
-    for (i, eq) in enumerate(eqs)
-        vars = Symbolics.get_variables(eq)
-        for var1 in vars
-            if Symbolics.iscall(var1) # skip e.g. constants
-                name = nameof(operation(var1))
-                var2, = @variables $name(iv2func)
-                eq = substitute(eq, var1 => var2; fold = false)
-            end
-        end
-        eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
-        eq = substitute(eq, Differential(iv1)(iv2func) => Differential(iv1)(iv2_of_iv1)) # substitute in d(iv2)/d(iv1)
-        eq = expand_derivatives(eq)
-        eq = substitute(eq, iv2func => iv2) # make iv2 independent
-        eq = substitute(eq, iv1 => iv1_of_iv2) # substitute any remaining old ivars
-        eqs[i] = eq
-    end
-
-    sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
-    return sys2
-end
-
-function change_independent_variable(sys::AbstractODESystem, iv; kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, kwargs...)
     iv1 = get_iv(sys) # e.g. t
-    iv2func = iv # e.g. a(t)
-    iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
-    iv2, = @independent_variables $iv2name
-
-    # TODO: find iv2func in system and replace it with some dummy variable
-    # TODO: not just to 1st, but to all orders
+    iv2name = nameof(iv) # TODO: handle namespacing?
+    iv2, = @independent_variables $iv2name # e.g. a
+    iv2func, = @variables $iv2name(iv1) # e.g. a(t)
+    D1 = Differential(iv1)
+    D2 = Differential(iv2)
 
     div2name = Symbol(iv2name, :_t)
-    div2, = @variables $div2name(iv2)
+    div2, = @variables $div2name(iv2) # e.g. a_t(a)
+
+    ddiv2name = Symbol(iv2name, :_tt)
+    ddiv2, = @variables $ddiv2name(iv2) # e.g. a_tt(a)
 
     eqs = ModelingToolkit.get_eqs(sys) |> copy # don't modify original system
+    !isnothing(eq) && push!(eqs, eq)
     vars = []
+    div2_div1 = nothing
     for (i, eq) in enumerate(eqs)
+        verbose && println("1. ", eq)
+
+        # 1) Substitute f(t₁) => f(t₂(t₁)) in all variables
         vars = Symbolics.get_variables(eq)
         for var1 in vars
             if Symbolics.iscall(var1) && !isequal(var1, iv2func) # && isequal(only(arguments(var1)), iv1) # skip e.g. constants
@@ -110,17 +86,32 @@ function change_independent_variable(sys::AbstractODESystem, iv; kwargs...)
                 eq = substitute(eq, var1 => var2; fold = false)
             end
         end
+        verbose && println("2. ", eq)
+
+        # 2) Substitute in dummy variables for dⁿt₂/dt₁ⁿ
         eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
-        #eq = substitute(eq, Differential(iv1)(iv2func) => Differential(iv1)(iv2_of_iv1)) # substitute in d(iv2)/d(iv1)
-        #eq = expand_derivatives(eq)
-        eq = substitute(eq, Differential(iv1)(Differential(iv1)(iv2func)) => div2) # e.g. D(a(t)) => a_t(t) # TODO: more orders
-        eq = substitute(eq, Differential(iv1)(iv2func) => div2) # e.g. D(a(t)) => a_t(t) # TODO: more orders
-        eq = substitute(eq, iv2func => iv2) # make iv2 independent
-        #eq = substitute(eq, iv1 => iv1_of_iv2) # substitute any remaining old ivars
+        verbose && println("3. ", eq)
+        eq = substitute(eq, D1(D1(iv2func)) => ddiv2) # order 2 # TODO: more orders
+        eq = substitute(eq, D1(iv2func) => div2) # order 1; e.g. D(a(t)) => a_t(t)
+        eq = substitute(eq, iv2func => iv2) # order 0; make iv2 independent
+        verbose && println("4. ", eq)
+        verbose && println()
+
         eqs[i] = eq
-        println(eq)
+
+        if isequal(eq.lhs, div2)
+            div2_div1 = eq.rhs
+        end
     end
 
+    isnothing(div2_div1) && error("No equation for $D1($iv2func) was specified.")
+    verbose && println("Found $div2 = $div2_div1")
+
+    # 3) Add equations for dummy variables
+    div1_div2 = 1 / div2_div1
+    push!(eqs, ddiv2 ~ expand_derivatives(-D2(div1_div2) / div1_div2^3)) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
+
+    # 4) Recreate system with new equations
     sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
     return sys2
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 3c7a4e6755..34f350222c 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -34,22 +34,23 @@ sol = solve(prob, Tsit5())
 @test sol[end, end] ≈ 1.0742818931017244
 
 @testset "Change independent variable" begin
-    # Autonomous 1st order (t doesn't appear in the equations)
     @independent_variables t
-    @variables x(t) y(t) z(t)
+    @variables x(t) y(t) z(t) s(t)
     eqs = [
         D(x) ~ y
         D(D(y)) ~ 2 * x * D(y)
         z ~ x + D(y)
+        D(s) ~ 1 / (2*s)
     ]
     @named sys1 = ODESystem(eqs, t)
+    sys1 = complete(sys1)
 
     @independent_variables s
-    @named sys2 = ModelingToolkit.change_independent_variable(sys1, s, s^2, √(t))
+    sys2 = ModelingToolkit.change_independent_variable(sys1, s)
 
-    sys1 = structural_simplify(sys1)
-    sys2 = structural_simplify(sys2)
-    prob1 = ODEProblem(sys1, [sys1.x => 1.0, sys1.y => 1.0, Differential(t)(sys1.y) => 0.0], (1.0, 4.0))
+    sys1 = structural_simplify(sys1; allow_symbolic = true)
+    sys2 = structural_simplify(sys2; allow_symbolic = true)
+    prob1 = ODEProblem(sys1, [sys1.x => 1.0, sys1.y => 1.0, Differential(t)(sys1.y) => 0.0, sys1.s => 1.0], (1.0, 4.0))
     prob2 = ODEProblem(sys2, [sys2.x => 1.0, sys2.y => 1.0, Differential(s)(sys2.y) => 0.0], (1.0, 2.0))
     sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
@@ -59,10 +60,10 @@ sol = solve(prob, Tsit5())
           all(isapprox.(sol1(ts, idxs=sys1.y), sol2(ss, idxs=sys2.y); atol = 1e-7))
 end
 
-#@testset "Change independent variable (Friedmann equation)" begin
+@testset "Change independent variable (Friedmann equation)" begin
     @independent_variables t
     D = Differential(t)
-    @variables a(t) ρr(t) ρm(t) ρΛ(t) ρ(t) ϕ(t)
+    @variables a(t) ρr(t) ρm(t) ρΛ(t) ρ(t) P(t) ϕ(t)
     @parameters Ωr0 Ωm0 ΩΛ0
     eqs = [
         ρr ~ 3/(8*Num(π)) * Ωr0 / a^4
@@ -75,7 +76,8 @@ end
     @named M1 = ODESystem(eqs, t)
     M1 = complete(M1)
 
-    M2 = ModelingToolkit.change_independent_variable(M1, M1.a)
-
+    @independent_variables a
+    M2 = ModelingToolkit.change_independent_variable(M1, a)
     M2 = structural_simplify(M2; allow_symbolic = true)
-#end
+    @test length(unknowns(M2)) == 2
+end

From af3ae69d0a307822a2212462eccf63a843451ffb Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 12:03:16 +0100
Subject: [PATCH 05/32] Fix broken Liouville transform test/documentation

---
 src/systems/diffeqs/basic_transformations.jl | 24 ++++-------
 test/basic_transformations.jl                | 44 ++++++++------------
 2 files changed, 27 insertions(+), 41 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 841c0bda6d..4640aa9b5f 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -14,26 +14,20 @@ the final value of `trJ` is the probability of ``u(t)``.
 Example:
 
 ```julia
-using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit, OrdinaryDiffEq
 
 @independent_variables t
 @parameters α β γ δ
 @variables x(t) y(t)
 D = Differential(t)
+eqs = [D(x) ~ α*x - β*x*y, D(y) ~ -δ*y + γ*x*y]
+@named sys = ODESystem(eqs, t)
 
-eqs = [D(x) ~ α*x - β*x*y,
-       D(y) ~ -δ*y + γ*x*y]
-
-sys = ODESystem(eqs)
 sys2 = liouville_transform(sys)
-@variables trJ
-
-u0 = [x => 1.0,
-      y => 1.0,
-      trJ => 1.0]
-
-prob = ODEProblem(complete(sys2),u0,tspan,p)
-sol = solve(prob,Tsit5())
+sys2 = complete(sys2)
+u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0]
+prob = ODEProblem(sys2, u0, tspan, p)
+sol = solve(prob, Tsit5())
 ```
 
 Where `sol[3,:]` is the evolution of `trJ` over time.
@@ -46,14 +40,14 @@ Optimal Transport Approach
 Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya
 https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf
 """
-function liouville_transform(sys::AbstractODESystem)
+function liouville_transform(sys::AbstractODESystem; kwargs...)
     t = get_iv(sys)
     @variables trJ
     D = ModelingToolkit.Differential(t)
     neweq = D(trJ) ~ trJ * -tr(calculate_jacobian(sys))
     neweqs = [equations(sys); neweq]
     vars = [unknowns(sys); trJ]
-    ODESystem(neweqs, t, vars, parameters(sys), checks = false)
+    ODESystem(neweqs, t, vars, parameters(sys); checks = false, name = nameof(sys), kwargs...)
 end
 
 function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, kwargs...)
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 34f350222c..2436dc7fc4 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -1,37 +1,29 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
 @independent_variables t
-@parameters α β γ δ
-@variables x(t) y(t)
 D = Differential(t)
 
-eqs = [D(x) ~ α * x - β * x * y,
-    D(y) ~ -δ * y + γ * x * y]
+@testset "Liouville transform" begin
+    @parameters α β γ δ
+    @variables x(t) y(t)
+    eqs = [D(x) ~ α * x - β * x * y, D(y) ~ -δ * y + γ * x * y]
+    @named sys = ODESystem(eqs, t)
+    sys = complete(sys)
 
-sys = ODESystem(eqs)
+    u0 = [x => 1.0, y => 1.0]
+    p = [α => 1.5, β => 1.0, δ => 3.0, γ => 1.0]
+    tspan = (0.0, 10.0)
+    prob = ODEProblem(sys, u0, tspan, p)
+    sol = solve(prob, Tsit5())
 
-u0 = [x => 1.0,
-    y => 1.0]
+    sys2 = liouville_transform(sys)
+    sys2 = complete(sys2)
 
-p = [α => 1.5,
-    β => 1.0,
-    δ => 3.0,
-    γ => 1.0]
-
-tspan = (0.0, 10.0)
-prob = ODEProblem(sys, u0, tspan, p)
-sol = solve(prob, Tsit5())
-
-sys2 = liouville_transform(sys)
-@variables trJ
-
-u0 = [x => 1.0,
-    y => 1.0,
-    trJ => 1.0]
-
-prob = ODEProblem(sys2, u0, tspan, p, jac = true)
-sol = solve(prob, Tsit5())
-@test sol[end, end] ≈ 1.0742818931017244
+    u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0]
+    prob = ODEProblem(sys2, u0, tspan, p, jac = true)
+    sol = solve(prob, Tsit5())
+    @test sol[end, end] ≈ 1.0742818931017244
+end
 
 @testset "Change independent variable" begin
     @independent_variables t

From cbceec481a050189015189edc56d4f879a112879 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 12:12:47 +0100
Subject: [PATCH 06/32] Specify new independent variable as a dependent
 variable in the old system

---
 src/systems/diffeqs/basic_transformations.jl |  2 +-
 test/basic_transformations.jl                | 24 ++++++++------------
 2 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 4640aa9b5f..19c997e787 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -52,7 +52,7 @@ end
 
 function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, kwargs...)
     iv1 = get_iv(sys) # e.g. t
-    iv2name = nameof(iv) # TODO: handle namespacing?
+    iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
     iv2, = @independent_variables $iv2name # e.g. a
     iv2func, = @variables $iv2name(iv1) # e.g. a(t)
     D1 = Differential(iv1)
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 2436dc7fc4..55f1612fd3 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -26,7 +26,6 @@ D = Differential(t)
 end
 
 @testset "Change independent variable" begin
-    @independent_variables t
     @variables x(t) y(t) z(t) s(t)
     eqs = [
         D(x) ~ y
@@ -34,26 +33,22 @@ end
         z ~ x + D(y)
         D(s) ~ 1 / (2*s)
     ]
-    @named sys1 = ODESystem(eqs, t)
-    sys1 = complete(sys1)
+    M1 = ODESystem(eqs, t; name = :M) |> complete
+    M2 = ModelingToolkit.change_independent_variable(M1, M1.s)
 
-    @independent_variables s
-    sys2 = ModelingToolkit.change_independent_variable(sys1, s)
-
-    sys1 = structural_simplify(sys1; allow_symbolic = true)
-    sys2 = structural_simplify(sys2; allow_symbolic = true)
-    prob1 = ODEProblem(sys1, [sys1.x => 1.0, sys1.y => 1.0, Differential(t)(sys1.y) => 0.0, sys1.s => 1.0], (1.0, 4.0))
-    prob2 = ODEProblem(sys2, [sys2.x => 1.0, sys2.y => 1.0, Differential(s)(sys2.y) => 0.0], (1.0, 2.0))
+    M1 = structural_simplify(M1; allow_symbolic = true)
+    M2 = structural_simplify(M2; allow_symbolic = true)
+    prob1 = ODEProblem(M1, [M1.x => 1.0, M1.y => 1.0, Differential(M1.t)(M1.y) => 0.0, M1.s => 1.0], (1.0, 4.0))
+    prob2 = ODEProblem(M2, [M2.x => 1.0, M2.y => 1.0, Differential(M2.s)(M2.y) => 0.0], (1.0, 2.0))
     sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     ts = range(0.0, 1.0, length = 50)
     ss = .√(ts)
-    @test all(isapprox.(sol1(ts, idxs=sys1.x), sol2(ss, idxs=sys2.x); atol = 1e-7)) &&
-          all(isapprox.(sol1(ts, idxs=sys1.y), sol2(ss, idxs=sys2.y); atol = 1e-7))
+    @test all(isapprox.(sol1(ts, idxs=M1.x), sol2(ss, idxs=M2.x); atol = 1e-7)) &&
+          all(isapprox.(sol1(ts, idxs=M1.y), sol2(ss, idxs=M2.y); atol = 1e-7))
 end
 
 @testset "Change independent variable (Friedmann equation)" begin
-    @independent_variables t
     D = Differential(t)
     @variables a(t) ρr(t) ρm(t) ρΛ(t) ρ(t) P(t) ϕ(t)
     @parameters Ωr0 Ωm0 ΩΛ0
@@ -68,8 +63,7 @@ end
     @named M1 = ODESystem(eqs, t)
     M1 = complete(M1)
 
-    @independent_variables a
-    M2 = ModelingToolkit.change_independent_variable(M1, a)
+    M2 = ModelingToolkit.change_independent_variable(M1, M1.a)
     M2 = structural_simplify(M2; allow_symbolic = true)
     @test length(unknowns(M2)) == 2
 end

From 75fe0f062dce5a3d89c71b8feffd3fd4ed1ff751 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 14:51:59 +0100
Subject: [PATCH 07/32] Optionally simplify dummy derivative expressions when
 changing independent variable

---
 src/systems/diffeqs/basic_transformations.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 19c997e787..a771f0e5f0 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -50,7 +50,7 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
     ODESystem(neweqs, t, vars, parameters(sys); checks = false, name = nameof(sys), kwargs...)
 end
 
-function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, simplify = true, kwargs...)
     iv1 = get_iv(sys) # e.g. t
     iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
     iv2, = @independent_variables $iv2name # e.g. a
@@ -103,7 +103,11 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; v
 
     # 3) Add equations for dummy variables
     div1_div2 = 1 / div2_div1
-    push!(eqs, ddiv2 ~ expand_derivatives(-D2(div1_div2) / div1_div2^3)) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
+    ddiv1_ddiv2 = expand_derivatives(-D2(div1_div2) / div1_div2^3)
+    if simplify
+        ddiv1_ddiv2 = Symbolics.simplify(ddiv1_ddiv2)
+    end
+    push!(eqs, ddiv2 ~ ddiv1_ddiv2) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
 
     # 4) Recreate system with new equations
     sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)

From b8c9839bbfc071b7b32c7a1a0386b0a69a86f135 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 15:34:25 +0100
Subject: [PATCH 08/32] Improve change of independent variable tests

---
 test/basic_transformations.jl | 33 ++++++++++++++++++++++++++++-----
 1 file changed, 28 insertions(+), 5 deletions(-)

diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 55f1612fd3..5f17abdd52 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -50,20 +50,43 @@ end
 
 @testset "Change independent variable (Friedmann equation)" begin
     D = Differential(t)
-    @variables a(t) ρr(t) ρm(t) ρΛ(t) ρ(t) P(t) ϕ(t)
+    @variables a(t) ȧ(t) ρr(t) ρm(t) ρΛ(t) ρ(t) P(t) ϕ(t)
     @parameters Ωr0 Ωm0 ΩΛ0
     eqs = [
         ρr ~ 3/(8*Num(π)) * Ωr0 / a^4
         ρm ~ 3/(8*Num(π)) * Ωm0 / a^3
         ρΛ ~ 3/(8*Num(π)) * ΩΛ0
         ρ ~ ρr + ρm + ρΛ
-        D(a) ~ √(8*Num(π)/3*ρ*a^4)
+        ȧ ~ √(8*Num(π)/3*ρ*a^4)
+        D(a) ~ ȧ
         D(D(ϕ)) ~ -3*D(a)/a*D(ϕ)
     ]
-    @named M1 = ODESystem(eqs, t)
-    M1 = complete(M1)
+    M1 = ODESystem(eqs, t; name = :M) |> complete
 
-    M2 = ModelingToolkit.change_independent_variable(M1, M1.a)
+    # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
+    M2 = ModelingToolkit.change_independent_variable(M1, M1.a) |> complete #, D(b) ~ D(a)/a; verbose = true)
+    @variables b(M2.a)
+    M3 = ModelingToolkit.change_independent_variable(M2, b, Differential(M2.a)(b) ~ exp(-b))
     M2 = structural_simplify(M2; allow_symbolic = true)
+    M3 = structural_simplify(M3; allow_symbolic = true)
     @test length(unknowns(M2)) == 2
 end
+
+@testset "Change independent variable (simple)" begin
+    @variables x(t)
+    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
+    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x)
+    @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2x, x_tt ~ 4x]))
+end
+
+@testset "Change independent variable (free fall)" 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) |> complete # gives (x, y) as function of t, ...
+    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x) # ... 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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 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.x/v)^2/2]; atol = 1e-10)) # eliminate t from anal solution x(t) = v*t, y(t) = v*t - g*t^2/2
+end

From dfc52694717c65710585e5aa96b24cc0ce0158f7 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 15:58:03 +0100
Subject: [PATCH 09/32] Explicitly insert dummy equations into system, if
 requested

---
 src/systems/diffeqs/basic_transformations.jl | 12 ++++++++++--
 test/basic_transformations.jl                | 15 ++++++++++++---
 2 files changed, 22 insertions(+), 5 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index a771f0e5f0..f0c8702e92 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -50,7 +50,7 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
     ODESystem(neweqs, t, vars, parameters(sys); checks = false, name = nameof(sys), kwargs...)
 end
 
-function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, simplify = true, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, simplify = true, dummies = false, kwargs...)
     iv1 = get_iv(sys) # e.g. t
     iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
     iv2, = @independent_variables $iv2name # e.g. a
@@ -109,7 +109,15 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; v
     end
     push!(eqs, ddiv2 ~ ddiv1_ddiv2) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
 
-    # 4) Recreate system with new equations
+    # 4) If requested, instead remove and insert dummy equations
+    if !dummies
+        dummyidxs = findall(eq -> isequal(eq.lhs, div2) || isequal(eq.lhs, ddiv2), eqs)
+        dummyeqs = splice!(eqs, dummyidxs) # return and remove dummy equations
+        dummysubs = Dict(eq.lhs => eq.rhs for eq in dummyeqs)
+        eqs = substitute.(eqs, Ref(dummysubs)) # don't iterate over dummysubs
+    end
+
+    # 5) Recreate system with new equations
     sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
     return sys2
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 5f17abdd52..add54a62c8 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -25,6 +25,15 @@ D = Differential(t)
     @test sol[end, end] ≈ 1.0742818931017244
 end
 
+@testset "Change independent variable (trivial)" begin
+    @variables x(t) y(t)
+    eqs1 = [D(D(x)) ~ D(x) + x, D(y) ~ 1]
+    M1 = ODESystem(eqs1, t; name = :M) |> complete
+    M2 = ModelingToolkit.change_independent_variable(M1, 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
+end
+
 @testset "Change independent variable" begin
     @variables x(t) y(t) z(t) s(t)
     eqs = [
@@ -75,7 +84,7 @@ end
 @testset "Change independent variable (simple)" begin
     @variables x(t)
     Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
-    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x)
+    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x; dummies = true)
     @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2x, x_tt ~ 4x]))
 end
 
@@ -83,10 +92,10 @@ end
     @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) |> complete # gives (x, y) as function of t, ...
-    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x) # ... but we want y as a function of x
+    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x; dummies = false) # ... 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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 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.x/v)^2/2]; atol = 1e-10)) # eliminate t from anal solution x(t) = v*t, y(t) = v*t - g*t^2/2
+    @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.x/v)^2/2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2)
 end

From ae9c174b06e4a5eb959767861576f81af08b02cb Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 18:09:03 +0100
Subject: [PATCH 10/32] Add and test errors when changing independent variable

---
 src/systems/diffeqs/basic_transformations.jl | 24 +++++++++++++++++---
 test/basic_transformations.jl                | 15 ++++++++++++
 2 files changed, 36 insertions(+), 3 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index f0c8702e92..4991ef73f1 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -51,10 +51,24 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
 end
 
 function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, simplify = true, dummies = false, kwargs...)
+    if !iscomplete(sys)
+        error("Cannot change independent variable of incomplete system $(nameof(sys))")
+    elseif isscheduled(sys)
+        error("Cannot change independent variable of structurally simplified system $(nameof(sys))")
+    elseif !isempty(get_systems(sys))
+        error("Cannot change independent variable of hierarchical system $(nameof(sys)). Flatten it first.") # TODO: implement
+    end
+
+    iv = unwrap(iv)
     iv1 = get_iv(sys) # e.g. t
-    iv2name = nameof(operation(unwrap(iv))) # TODO: handle namespacing?
+
+    if !iscall(iv) || !isequal(only(arguments(iv)), iv1)
+        error("New independent variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))")
+    end
+
+    iv2func = iv # e.g. a(t)
+    iv2name = nameof(operation(iv))
     iv2, = @independent_variables $iv2name # e.g. a
-    iv2func, = @variables $iv2name(iv1) # e.g. a(t)
     D1 = Differential(iv1)
     D2 = Differential(iv2)
 
@@ -98,8 +112,12 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; v
         end
     end
 
-    isnothing(div2_div1) && error("No equation for $D1($iv2func) was specified.")
     verbose && println("Found $div2 = $div2_div1")
+    if isnothing(div2_div1)
+        error("No equation for $D1($iv2func) was specified.")
+    elseif isequal(div2_div1, 0)
+        error("Cannot change independent variable from $iv1 to $iv2 with singular transformation $(div2 ~ div2_div1).")
+    end
 
     # 3) Add equations for dummy variables
     div1_div2 = 1 / div2_div1
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index add54a62c8..992bbfcc48 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -99,3 +99,18 @@ end
     sol = solve(prob, Tsit5(); reltol = 1e-5)
     @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.x/v)^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 (errors)" begin
+    @variables x(t) y z(y) w(t) v(t)
+    M = ODESystem([D(x) ~ 0, v ~ x], t; name = :M)
+    @test_throws "incomplete" ModelingToolkit.change_independent_variable(M, M.x)
+    M = complete(M)
+    @test_throws "singular" ModelingToolkit.change_independent_variable(M, M.x)
+    @test_throws "structurally simplified" ModelingToolkit.change_independent_variable(structural_simplify(M), y)
+    @test_throws "No equation" ModelingToolkit.change_independent_variable(M, w)
+    @test_throws "No equation" ModelingToolkit.change_independent_variable(M, v)
+    @test_throws "not a function of the independent variable" ModelingToolkit.change_independent_variable(M, y)
+    @test_throws "not a function of the independent variable" ModelingToolkit.change_independent_variable(M, z)
+    M = compose(M, M)
+    @test_throws "hierarchical" ModelingToolkit.change_independent_variable(M, M.x)
+end

From 00b27111ddd3ec9b00ebd620cab0c3e5b6fd8116 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 2 Mar 2025 19:45:22 +0100
Subject: [PATCH 11/32] Export and document change_independent_variable

---
 docs/src/systems/ODESystem.md                |  1 +
 src/ModelingToolkit.jl                       |  2 +-
 src/systems/diffeqs/basic_transformations.jl | 41 +++++++++++++++++++-
 test/basic_transformations.jl                | 36 ++++++++---------
 4 files changed, 60 insertions(+), 20 deletions(-)

diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md
index 6cc34725c4..4be0f6e263 100644
--- a/docs/src/systems/ODESystem.md
+++ b/docs/src/systems/ODESystem.md
@@ -30,6 +30,7 @@ ODESystem
 structural_simplify
 ode_order_lowering
 dae_index_lowering
+change_independent_variable
 liouville_transform
 alias_elimination
 tearing
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index e6b8130af3..51156ba993 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -255,7 +255,7 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        tunable_parameters, isirreducible, getdescription, hasdescription,
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
-export ode_order_lowering, dae_order_lowering, liouville_transform
+export ode_order_lowering, dae_order_lowering, liouville_transform, change_independent_variable
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 4991ef73f1..665f150407 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -50,7 +50,46 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
     ODESystem(neweqs, t, vars, parameters(sys); checks = false, name = nameof(sys), kwargs...)
 end
 
-function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; verbose = false, simplify = true, dummies = false, kwargs...)
+"""
+    function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; dummies = false, simplify = true, verbose = false, kwargs...)
+
+Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``f(t)``).
+An equation in `sys` must define the rate of change of the new independent variable (e.g. ``df(t)/dt``).
+Alternatively, `eq` can specify such an equation.
+
+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. ``df(t)/dt > 0`` or ``df(t)/dt < 0``).
+
+Keyword arguments
+=================
+If `dummies`, derivatives of the new independent variable are expressed through dummy equations; otherwise they are explicitly inserted into the equations.
+If `simplify`, these dummy expressions are simplified and often give a tidier transformation.
+If `verbose`, the function prints intermediate transformations of equations to aid debugging.
+Any additional keyword arguments `kwargs...` are forwarded to the constructor that rebuilds the system.
+
+Usage before structural simplification
+======================================
+The variable change must take place before structural simplification.
+Subsequently, consider passing `allow_symbolic = true` to `structural_simplify(sys)` to reduce the number of unknowns, with the understanding that the transformation is well-defined.
+
+Example
+=======
+Consider a free fall with constant horizontal velocity.
+The laws of physics naturally describes position as a function of time.
+By changing the independent variable, it can be reformulated for vertical position as a function of horizontal distance:
+```julia
+julia> @variables x(t) y(t);
+
+julia> @named M = ODESystem([D(D(y)) ~ -9.81, D(x) ~ 10.0], t);
+
+julia> M′ = change_independent_variable(complete(M), x);
+
+julia> unknowns(M′)
+1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
+ y(x)
+```
+"""
+function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; dummies = false, simplify = true, verbose = false, kwargs...)
     if !iscomplete(sys)
         error("Cannot change independent variable of incomplete system $(nameof(sys))")
     elseif isscheduled(sys)
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 992bbfcc48..2e7f0ce54b 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -29,7 +29,7 @@ end
     @variables x(t) y(t)
     eqs1 = [D(D(x)) ~ D(x) + x, D(y) ~ 1]
     M1 = ODESystem(eqs1, t; name = :M) |> complete
-    M2 = ModelingToolkit.change_independent_variable(M1, M1.y)
+    M2 = change_independent_variable(M1, 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
 end
@@ -43,10 +43,10 @@ end
         D(s) ~ 1 / (2*s)
     ]
     M1 = ODESystem(eqs, t; name = :M) |> complete
-    M2 = ModelingToolkit.change_independent_variable(M1, M1.s)
+    M2 = change_independent_variable(M1, M1.s)
 
-    M1 = structural_simplify(M1; allow_symbolic = true)
-    M2 = structural_simplify(M2; allow_symbolic = true)
+    M1 = structural_simplify(M1)
+    M2 = structural_simplify(M2)
     prob1 = ODEProblem(M1, [M1.x => 1.0, M1.y => 1.0, Differential(M1.t)(M1.y) => 0.0, M1.s => 1.0], (1.0, 4.0))
     prob2 = ODEProblem(M2, [M2.x => 1.0, M2.y => 1.0, Differential(M2.s)(M2.y) => 0.0], (1.0, 2.0))
     sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
@@ -73,18 +73,18 @@ end
     M1 = ODESystem(eqs, t; name = :M) |> complete
 
     # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
-    M2 = ModelingToolkit.change_independent_variable(M1, M1.a) |> complete #, D(b) ~ D(a)/a; verbose = true)
+    M2 = change_independent_variable(M1, M1.a) |> complete #, D(b) ~ D(a)/a; verbose = true)
     @variables b(M2.a)
-    M3 = ModelingToolkit.change_independent_variable(M2, b, Differential(M2.a)(b) ~ exp(-b))
+    M3 = change_independent_variable(M2, b, Differential(M2.a)(b) ~ exp(-b))
     M2 = structural_simplify(M2; allow_symbolic = true)
     M3 = structural_simplify(M3; allow_symbolic = true)
-    @test length(unknowns(M2)) == 2
+    @test length(unknowns(M2)) == 2 && length(unknowns(M3)) == 2
 end
 
 @testset "Change independent variable (simple)" begin
     @variables x(t)
-    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
-    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x; dummies = true)
+    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete # TODO: avoid complete. can avoid it if passing defined $variable directly to change_independent_variable
+    Mx = change_independent_variable(Mt, Mt.x; dummies = true)
     @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2x, x_tt ~ 4x]))
 end
 
@@ -92,7 +92,7 @@ end
     @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) |> complete # gives (x, y) as function of t, ...
-    Mx = ModelingToolkit.change_independent_variable(Mt, Mt.x; dummies = false) # ... but we want y as a function of x
+    Mx = change_independent_variable(Mt, Mt.x; dummies = false) # ... 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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities
@@ -103,14 +103,14 @@ 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 "incomplete" ModelingToolkit.change_independent_variable(M, M.x)
+    @test_throws "incomplete" change_independent_variable(M, M.x)
     M = complete(M)
-    @test_throws "singular" ModelingToolkit.change_independent_variable(M, M.x)
-    @test_throws "structurally simplified" ModelingToolkit.change_independent_variable(structural_simplify(M), y)
-    @test_throws "No equation" ModelingToolkit.change_independent_variable(M, w)
-    @test_throws "No equation" ModelingToolkit.change_independent_variable(M, v)
-    @test_throws "not a function of the independent variable" ModelingToolkit.change_independent_variable(M, y)
-    @test_throws "not a function of the independent variable" ModelingToolkit.change_independent_variable(M, z)
+    @test_throws "singular" change_independent_variable(M, M.x)
+    @test_throws "structurally simplified" change_independent_variable(structural_simplify(M), y)
+    @test_throws "No equation" change_independent_variable(M, w)
+    @test_throws "No equation" change_independent_variable(M, v)
+    @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 = compose(M, M)
-    @test_throws "hierarchical" ModelingToolkit.change_independent_variable(M, M.x)
+    @test_throws "hierarchical" change_independent_variable(M, M.x)
 end

From ee212232247ef248c969ae6dcc84be27dbebfb02 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Mon, 3 Mar 2025 16:45:27 +0100
Subject: [PATCH 12/32] Handle autonomous systems and more fields

---
 src/systems/diffeqs/basic_transformations.jl | 66 +++++++++++---------
 test/basic_transformations.jl                | 21 ++++---
 2 files changed, 49 insertions(+), 38 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 665f150407..56d3ada089 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -51,7 +51,7 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
 end
 
 """
-    function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; dummies = false, simplify = true, verbose = false, kwargs...)
+    function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, verbose = false, kwargs...)
 
 Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``f(t)``).
 An equation in `sys` must define the rate of change of the new independent variable (e.g. ``df(t)/dt``).
@@ -89,7 +89,7 @@ julia> unknowns(M′)
  y(x)
 ```
 """
-function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; dummies = false, simplify = true, verbose = false, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, verbose = false, kwargs...)
     if !iscomplete(sys)
         error("Cannot change independent variable of incomplete system $(nameof(sys))")
     elseif isscheduled(sys)
@@ -100,9 +100,10 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; d
 
     iv = unwrap(iv)
     iv1 = get_iv(sys) # e.g. t
-
     if !iscall(iv) || !isequal(only(arguments(iv)), iv1)
         error("New independent variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))")
+    elseif !isautonomous(sys) && isempty(findall(eq -> isequal(eq.lhs, iv1), eqs))
+        error("System $(nameof(sys)) is autonomous in $iv1. An equation of the form $iv1 ~ F($iv) must be provided.")
     end
 
     iv2func = iv # e.g. a(t)
@@ -111,50 +112,53 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; d
     D1 = Differential(iv1)
     D2 = Differential(iv2)
 
+    iv1name = nameof(iv1) # e.g. t
+    iv1func, = @variables $iv1name(iv2) # e.g. t(a)
+
     div2name = Symbol(iv2name, :_t)
     div2, = @variables $div2name(iv2) # e.g. a_t(a)
 
     ddiv2name = Symbol(iv2name, :_tt)
     ddiv2, = @variables $ddiv2name(iv2) # e.g. a_tt(a)
 
-    eqs = ModelingToolkit.get_eqs(sys) |> copy # don't modify original system
-    !isnothing(eq) && push!(eqs, eq)
-    vars = []
-    div2_div1 = nothing
-    for (i, eq) in enumerate(eqs)
-        verbose && println("1. ", eq)
-
+    function transform(ex)
         # 1) Substitute f(t₁) => f(t₂(t₁)) in all variables
-        vars = Symbolics.get_variables(eq)
+        verbose && println("1. ", ex)
+        vars = Symbolics.get_variables(ex)
         for var1 in vars
             if Symbolics.iscall(var1) && !isequal(var1, iv2func) # && isequal(only(arguments(var1)), iv1) # skip e.g. constants
                 name = nameof(operation(var1))
                 var2, = @variables $name(iv2func)
-                eq = substitute(eq, var1 => var2; fold = false)
+                ex = substitute(ex, var1 => var2; fold = false)
             end
         end
-        verbose && println("2. ", eq)
 
         # 2) Substitute in dummy variables for dⁿt₂/dt₁ⁿ
-        eq = expand_derivatives(eq) # expand out with chain rule to get d(iv2)/d(iv1)
-        verbose && println("3. ", eq)
-        eq = substitute(eq, D1(D1(iv2func)) => ddiv2) # order 2 # TODO: more orders
-        eq = substitute(eq, D1(iv2func) => div2) # order 1; e.g. D(a(t)) => a_t(t)
-        eq = substitute(eq, iv2func => iv2) # order 0; make iv2 independent
-        verbose && println("4. ", eq)
-        verbose && println()
-
-        eqs[i] = eq
-
-        if isequal(eq.lhs, div2)
-            div2_div1 = eq.rhs
-        end
+        verbose && println("2. ", ex)
+        ex = expand_derivatives(ex) # expand out with chain rule to get d(iv2)/d(iv1)
+        verbose && println("3. ", ex)
+        ex = substitute(ex, D1(D1(iv2func)) => ddiv2) # order 2 # TODO: more orders
+        ex = substitute(ex, D1(iv2func) => div2) # order 1; e.g. D(a(t)) => a_t(t)
+        ex = substitute(ex, iv2func => iv2) # order 0; make iv2 independent
+        verbose && println("4. ", ex)
+        ex = substitute(ex, iv1 => iv1func) # if autonomous
+        verbose && println("5. ", ex)
+        return ex
     end
 
-    verbose && println("Found $div2 = $div2_div1")
-    if isnothing(div2_div1)
+    eqs = [get_eqs(sys); eqs]
+    eqs = map(transform, eqs)
+
+    initialization_eqs = map(transform, get_initialization_eqs(sys))
+
+    div2_div1_idxs = findall(eq -> isequal(eq.lhs, div2), eqs)
+    if length(div2_div1_idxs) == 0
         error("No equation for $D1($iv2func) was specified.")
-    elseif isequal(div2_div1, 0)
+    elseif length(div2_div1_idxs) >= 2
+        error("Multiple equations for $D1($iv2func) were specified.")
+    end
+    div2_div1 = eqs[only(div2_div1_idxs)].rhs
+    if isequal(div2_div1, 0)
         error("Cannot change independent variable from $iv1 to $iv2 with singular transformation $(div2 ~ div2_div1).")
     end
 
@@ -167,7 +171,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; d
     push!(eqs, ddiv2 ~ ddiv1_ddiv2) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
 
     # 4) If requested, instead remove and insert dummy equations
-    if !dummies
+    if !dummies # TODO: make dummies = false work with more fields!
         dummyidxs = findall(eq -> isequal(eq.lhs, div2) || isequal(eq.lhs, ddiv2), eqs)
         dummyeqs = splice!(eqs, dummyidxs) # return and remove dummy equations
         dummysubs = Dict(eq.lhs => eq.rhs for eq in dummyeqs)
@@ -175,6 +179,6 @@ function change_independent_variable(sys::AbstractODESystem, iv, eq = nothing; d
     end
 
     # 5) Recreate system with new equations
-    sys2 = typeof(sys)(eqs, iv2; name = nameof(sys), description = description(sys), kwargs...)
+    sys2 = typeof(sys)(eqs, iv2; initialization_eqs, name = nameof(sys), description = description(sys), kwargs...)
     return sys2
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 2e7f0ce54b..dc9fa21981 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -42,13 +42,14 @@ end
         z ~ x + D(y)
         D(s) ~ 1 / (2*s)
     ]
-    M1 = ODESystem(eqs, t; name = :M) |> complete
-    M2 = change_independent_variable(M1, M1.s)
+    initialization_eqs = [x ~ 1.0, y ~ 1.0, D(y) ~ 0.0]
+    M1 = ODESystem(eqs, t; initialization_eqs, name = :M) |> complete
+    M2 = change_independent_variable(M1, M1.s; dummies = true)
 
-    M1 = structural_simplify(M1)
-    M2 = structural_simplify(M2)
-    prob1 = ODEProblem(M1, [M1.x => 1.0, M1.y => 1.0, Differential(M1.t)(M1.y) => 0.0, M1.s => 1.0], (1.0, 4.0))
-    prob2 = ODEProblem(M2, [M2.x => 1.0, M2.y => 1.0, Differential(M2.s)(M2.y) => 0.0], (1.0, 2.0))
+    M1 = structural_simplify(M1; allow_symbolic = true)
+    M2 = structural_simplify(M2; allow_symbolic = true)
+    prob1 = ODEProblem(M1, [M1.s => 1.0], (1.0, 4.0), [])
+    prob2 = ODEProblem(M2, [], (1.0, 2.0), [])
     sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     ts = range(0.0, 1.0, length = 50)
@@ -75,7 +76,7 @@ end
     # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
     M2 = change_independent_variable(M1, M1.a) |> complete #, D(b) ~ D(a)/a; verbose = true)
     @variables b(M2.a)
-    M3 = change_independent_variable(M2, b, Differential(M2.a)(b) ~ exp(-b))
+    M3 = change_independent_variable(M2, b, [Differential(M2.a)(b) ~ exp(-b), M2.a ~ exp(b)])
     M2 = structural_simplify(M2; allow_symbolic = true)
     M3 = structural_simplify(M3; allow_symbolic = true)
     @test length(unknowns(M2)) == 2 && length(unknowns(M3)) == 2
@@ -100,6 +101,12 @@ end
     @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.x/v)^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 (autonomous system)" begin
+    M = ODESystem([D(x) ~ t], t; name = :M) |> complete # non-autonomous
+    @test_throws "t ~ F(x(t)) must be provided" change_independent_variable(M, M.x)
+    @test_nowarn change_independent_variable(M, M.x, [t ~ 2*x])
+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)

From 9444d9395ac795042023687822f7a7a878db8be0 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Mon, 3 Mar 2025 16:48:18 +0100
Subject: [PATCH 13/32] Add tutorial for changing independent variable

---
 docs/pages.jl                                 |   1 +
 .../tutorials/change_independent_variable.md  | 108 ++++++++++++++++++
 2 files changed, 109 insertions(+)
 create mode 100644 docs/src/tutorials/change_independent_variable.md

diff --git a/docs/pages.jl b/docs/pages.jl
index f6c49a0de3..829c0d1c08 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -10,6 +10,7 @@ pages = [
         "tutorials/stochastic_diffeq.md",
         "tutorials/discrete_system.md",
         "tutorials/parameter_identifiability.md",
+        "tutorials/change_independent_variable.md",
         "tutorials/bifurcation_diagram_computation.md",
         "tutorials/SampledData.md",
         "tutorials/domain_connections.md",
diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
new file mode 100644
index 0000000000..97ca77fef8
--- /dev/null
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -0,0 +1,108 @@
+# Changing the independent variable of ODEs
+
+Ordinary differential equations describe the rate of change of some dependent variables with respect to one independent variable.
+For the modeler it is often most natural to write down the equations with a particular independent variable, say time $t$.
+In many cases there are good reasons for reparametrizing ODEs in terms of a different independent variable:
+
+1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
+2. Some differential equations vary more nicely (e.g. less stiff or better behaved) with respect to one independent variable than another.
+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).
+
+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.
+This is mechanical and error-prone.
+ModelingToolkit provides the utility function [`change_independent_variable`](@ref) that automates this process.
+
+## 1. Get one dependent variable as a function of another
+
+Consider a projectile shot with some initial velocity in a gravitational field.
+```@example changeivar
+using ModelingToolkit
+@independent_variables t
+D = Differential(t)
+@variables x(t) y(t)
+@parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
+M1 = ODESystem([
+    D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
+], t; initialization_eqs = [
+    #x ~ 0.0, # TODO: handle?
+    y ~ 0.0, D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
+], name = :M) |> complete
+M1s = structural_simplify(M1)
+```
+This is the standard parametrization that arises naturally from kinematics and Newton's laws.
+It expresses the position $(x(t), y(t))$ as a function of time $t$.
+But suppose we want to determine whether the projectile hits a target 10 meters away.
+There are at least three ways of answering this:
+* 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.
+* 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.
+* Solve the ODE for $y(x)$ and evaluate it at 10 meters.
+
+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, M1.x; dummies = true) |> complete
+@assert M2.x isa Num # hide
+M2s = structural_simplify(M2; allow_symbolic = true)
+```
+The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`.
+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.
+It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
+```@example changeivar
+using OrdinaryDiffEq, Plots
+prob = ODEProblem(M2s, [], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
+sol = solve(prob, Tsit5())
+plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
+```
+
+## 2. Reduce stiffness by changing to a logarithmic time axis
+
+In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe.
+In terms of conformal time $t$, they can be written
+```@example changeivar
+@variables a(t) Ω(t) Ωr(t) Ωm(t) ΩΛ(t)
+M1 = ODESystem([
+    D(a) ~ √(Ω) * a^2
+    Ω ~ Ωr + Ωm + ΩΛ
+    D(Ωm) ~ -3 * D(a)/a * Ωm
+    D(Ωr) ~ -4 * D(a)/a * Ωr
+    D(ΩΛ) ~ 0
+], t; initialization_eqs = [
+    ΩΛ + Ωr + Ωm ~ 1
+], name = :M) |> complete
+M1s = structural_simplify(M1)
+```
+Of course, we can solve this ODE as it is:
+```@example changeivar
+prob = ODEProblem(M1s, [M1.a => 1.0, M1.Ωr => 5e-5, M1.Ωm => 0.3], (0.0, -3.3), [])
+sol = solve(prob, Tsit5(); reltol = 1e-5)
+@assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide
+plot(sol, idxs = [M1.a, M1.Ωr/M1.Ω, M1.Ωm/M1.Ω, M1.ΩΛ/M1.Ω])
+```
+The solver becomes unstable due to stiffness.
+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.
+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.
+
+It is therefore common to write these ODEs in terms of $b = \ln a$.
+To do this, we will change the independent variable in two stages; from $t$ to $a$ to $b$.
+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) |> complete
+```
+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$:
+```@example changeivar
+a = M2.a
+Da = Differential(a)
+@variables b(a)
+M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)]) |> complete
+```
+We can now solve and plot the ODE in terms of $b$:
+```@example changeivar
+M3s = structural_simplify(M3; allow_symbolic = true)
+prob = ODEProblem(M3s, [M3.Ωr => 5e-5, M3.Ωm => 0.3], (0, -15), [])
+sol = solve(prob, Tsit5())
+@assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide
+plot(sol, idxs = [M3.Ωr/M3.Ω, M3.Ωm/M3.Ω, M3.ΩΛ/M3.Ω])
+```
+Notice that the variables vary "more nicely" with respect to $b$ than $t$, making the solver happier and avoiding numerical problems.

From 98318864dbfa9fc0ac0a4721ebb25433bd3524b8 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Mon, 3 Mar 2025 21:31:23 +0100
Subject: [PATCH 14/32] Reorder things and make universal function that
 transforms all ODESystem fields

---
 .../tutorials/change_independent_variable.md  |   6 +-
 src/systems/diffeqs/basic_transformations.jl  | 120 +++++++++++-------
 test/basic_transformations.jl                 |   4 +-
 3 files changed, 77 insertions(+), 53 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index 97ca77fef8..6274465264 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -23,9 +23,11 @@ D = Differential(t)
 @parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
 M1 = ODESystem([
     D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
-], t; initialization_eqs = [
+], t; defaults = [
+    y => 0.0
+], initialization_eqs = [
     #x ~ 0.0, # TODO: handle?
-    y ~ 0.0, D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
+    D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
 ], name = :M) |> complete
 M1s = structural_simplify(M1)
 ```
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 56d3ada089..396dbd51f5 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -110,20 +110,59 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     iv2name = nameof(operation(iv))
     iv2, = @independent_variables $iv2name # e.g. a
     D1 = Differential(iv1)
-    D2 = Differential(iv2)
 
     iv1name = nameof(iv1) # e.g. t
     iv1func, = @variables $iv1name(iv2) # e.g. t(a)
 
-    div2name = Symbol(iv2name, :_t)
-    div2, = @variables $div2name(iv2) # e.g. a_t(a)
+    eqs = [get_eqs(sys); eqs] # copies system equations to avoid modifying original system
 
-    ddiv2name = Symbol(iv2name, :_tt)
-    ddiv2, = @variables $ddiv2name(iv2) # e.g. a_tt(a)
+    # 1) Find and compute all necessary expressions for e.g. df/dt, d²f/dt², ...
+    # 1.1) Find the 1st order derivative of the new independent variable (e.g. da(t)/dt = ...), ...
+    div2_div1_idxs = findall(eq -> isequal(eq.lhs, D1(iv2func)), eqs) # index of e.g. da/dt = ...
+    if length(div2_div1_idxs) != 1
+        error("Exactly one equation for $D1($iv2func) was not specified.")
+    end
+    div2_div1_eq = popat!(eqs, only(div2_div1_idxs)) # get and remove e.g. df/dt = ... (may be added back later)
+    div2_div1 = div2_div1_eq.rhs
+    if isequal(div2_div1, 0)
+        error("Cannot change independent variable from $iv1 to $iv2 with singular transformation $div2_div1_eq.")
+    end
+    # 1.2) ... then compute the 2nd order derivative of the new independent variable
+    div1_div2 = 1 / div2_div1 # TODO: URL reference for clarity
+    ddiv2_ddiv1 = expand_derivatives(-Differential(iv2func)(div1_div2) / div1_div2^3, simplify) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders # TODO: pass simplify here
+    # 1.3) # TODO: handle higher orders (3+) derivatives ...
+
+    # 2) If requested, insert extra dummy equations for e.g. df/dt, d²f/dt², ...
+    #    Otherwise, replace all these derivatives by their explicit expressions
+    if dummies
+        div2name = Symbol(iv2name, :_t) # TODO: not always t
+        div2, = @variables $div2name(iv2) # e.g. a_t(a)
+        ddiv2name = Symbol(iv2name, :_tt) # TODO: not always t
+        ddiv2, = @variables $ddiv2name(iv2) # e.g. a_tt(a)
+        eqs = [eqs; [div2 ~ div2_div1, ddiv2 ~ ddiv2_ddiv1]] # add dummy equations
+        derivsubs = [D1(D1(iv2func)) => ddiv2, D1(iv2func) => div2] # order is crucial!
+    else
+        derivsubs = [D1(D1(iv2func)) => ddiv2_ddiv1, D1(iv2func) => div2_div1] # order is crucial!
+    end
+    derivsubs = [derivsubs; [iv2func => iv2, iv1 => iv1func]]
+
+    if verbose
+        # Explain what we just did
+        println("Order 1 (found):    $div2_div1_eq")
+        println("Order 2 (computed): $(D1(div2_div1_eq.lhs) ~ ddiv2_ddiv1)")
+        println()
+        println("Substitutions will be made in this order:")
+        for (n, sub) in enumerate(derivsubs)
+            println("$n: $(sub[1]) => $(sub[2])")
+        end
+        println()
+    end
 
+    # 3) Define a transformation function that performs the change of variable on any expression/equation
     function transform(ex)
-        # 1) Substitute f(t₁) => f(t₂(t₁)) in all variables
-        verbose && println("1. ", ex)
+        verbose && println("Step 0: ", ex)
+
+        # Step 1: substitute f(t₁) => f(t₂(t₁)) in all variables in the expression
         vars = Symbolics.get_variables(ex)
         for var1 in vars
             if Symbolics.iscall(var1) && !isequal(var1, iv2func) # && isequal(only(arguments(var1)), iv1) # skip e.g. constants
@@ -132,53 +171,36 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
                 ex = substitute(ex, var1 => var2; fold = false)
             end
         end
+        verbose && println("Step 1: ", ex)
 
-        # 2) Substitute in dummy variables for dⁿt₂/dt₁ⁿ
-        verbose && println("2. ", ex)
+        # Step 2: expand out all chain rule derivatives
         ex = expand_derivatives(ex) # expand out with chain rule to get d(iv2)/d(iv1)
-        verbose && println("3. ", ex)
-        ex = substitute(ex, D1(D1(iv2func)) => ddiv2) # order 2 # TODO: more orders
-        ex = substitute(ex, D1(iv2func) => div2) # order 1; e.g. D(a(t)) => a_t(t)
-        ex = substitute(ex, iv2func => iv2) # order 0; make iv2 independent
-        verbose && println("4. ", ex)
-        ex = substitute(ex, iv1 => iv1func) # if autonomous
-        verbose && println("5. ", ex)
+        verbose && println("Step 2: ", ex)
+
+        # Step 3: substitute d²f/dt², df/dt, ... (to dummy variables or explicit expressions, depending on dummies)
+        for sub in derivsubs
+            ex = substitute(ex, sub)
+        end
+        verbose && println("Step 3: ", ex)
+        verbose && println()
+
         return ex
     end
 
-    eqs = [get_eqs(sys); eqs]
+    # 4) Transform all fields
     eqs = map(transform, eqs)
-
+    observed = map(transform, get_observed(sys))
     initialization_eqs = map(transform, get_initialization_eqs(sys))
-
-    div2_div1_idxs = findall(eq -> isequal(eq.lhs, div2), eqs)
-    if length(div2_div1_idxs) == 0
-        error("No equation for $D1($iv2func) was specified.")
-    elseif length(div2_div1_idxs) >= 2
-        error("Multiple equations for $D1($iv2func) were specified.")
-    end
-    div2_div1 = eqs[only(div2_div1_idxs)].rhs
-    if isequal(div2_div1, 0)
-        error("Cannot change independent variable from $iv1 to $iv2 with singular transformation $(div2 ~ div2_div1).")
-    end
-
-    # 3) Add equations for dummy variables
-    div1_div2 = 1 / div2_div1
-    ddiv1_ddiv2 = expand_derivatives(-D2(div1_div2) / div1_div2^3)
-    if simplify
-        ddiv1_ddiv2 = Symbolics.simplify(ddiv1_ddiv2)
-    end
-    push!(eqs, ddiv2 ~ ddiv1_ddiv2) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders
-
-    # 4) If requested, instead remove and insert dummy equations
-    if !dummies # TODO: make dummies = false work with more fields!
-        dummyidxs = findall(eq -> isequal(eq.lhs, div2) || isequal(eq.lhs, ddiv2), eqs)
-        dummyeqs = splice!(eqs, dummyidxs) # return and remove dummy equations
-        dummysubs = Dict(eq.lhs => eq.rhs for eq in dummyeqs)
-        eqs = substitute.(eqs, Ref(dummysubs)) # don't iterate over dummysubs
-    end
-
-    # 5) Recreate system with new equations
-    sys2 = typeof(sys)(eqs, iv2; initialization_eqs, name = nameof(sys), description = description(sys), kwargs...)
-    return sys2
+    parameter_dependencies = map(transform, get_parameter_dependencies(sys))
+    defaults = Dict(transform(var) => transform(val) for (var, val) in get_defaults(sys))
+    guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
+    assertions = Dict(transform(condition) => msg for (condition, msg) in get_assertions(sys))
+    # TODO: handle subsystems
+
+    # 5) Recreate system with transformed fields
+    return typeof(sys)(
+        eqs, iv2;
+        observed, initialization_eqs, parameter_dependencies, defaults, guesses, assertions,
+        name = nameof(sys), description = description(sys), kwargs...
+    ) |> complete # original system had to be complete
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index dc9fa21981..99a252b5b8 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -114,8 +114,8 @@ end
     M = complete(M)
     @test_throws "singular" change_independent_variable(M, M.x)
     @test_throws "structurally simplified" change_independent_variable(structural_simplify(M), y)
-    @test_throws "No equation" change_independent_variable(M, w)
-    @test_throws "No equation" change_independent_variable(M, v)
+    @test_throws "not specified" change_independent_variable(M, w)
+    @test_throws "not specified" change_independent_variable(M, v)
     @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 = compose(M, M)

From 1ac7c7d3c221b115e58ce4ed52d010926da22b90 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Tue, 4 Mar 2025 22:37:37 +0100
Subject: [PATCH 15/32] Clean up independent variable change implementation

---
 .../tutorials/change_independent_variable.md  |   2 +-
 src/systems/diffeqs/basic_transformations.jl  | 167 ++++++++----------
 test/basic_transformations.jl                 |  55 ++++--
 3 files changed, 115 insertions(+), 109 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index 6274465264..0d4d8bcac8 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -26,7 +26,7 @@ M1 = ODESystem([
 ], t; defaults = [
     y => 0.0
 ], initialization_eqs = [
-    #x ~ 0.0, # TODO: handle?
+    #x ~ 0.0, # TODO: handle? # hide
     D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
 ], name = :M) |> complete
 M1s = structural_simplify(M1)
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 396dbd51f5..11a5514f42 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -51,32 +51,37 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
 end
 
 """
-    function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, verbose = false, kwargs...)
+    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = 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. ``f(t)``).
-An equation in `sys` must define the rate of change of the new independent variable (e.g. ``df(t)/dt``).
-Alternatively, `eq` can specify such an equation.
+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. ``df(t)/dt > 0`` or ``df(t)/dt < 0``).
+This is satisfied if one is a strictly increasing function of the other (e.g. ``du(t)/dt > 0`` or ``du(t)/dt < 0``).
 
 Keyword arguments
 =================
-If `dummies`, derivatives of the new independent variable are expressed through dummy equations; otherwise they are explicitly inserted into the equations.
+If `dummies`, derivatives of the new independent variable with respect to the old one are expressed through dummy equations; otherwise they are explicitly inserted into the equations.
 If `simplify`, these dummy expressions are simplified and often give a tidier transformation.
-If `verbose`, the function prints intermediate transformations of equations to aid debugging.
-Any additional keyword arguments `kwargs...` are forwarded to the constructor that rebuilds the system.
+If `fold`, internal substitutions will evaluate numerical expressions.
+Additional keyword arguments `kwargs...` are forwarded to the constructor that rebuilds `sys`.
 
 Usage before structural simplification
 ======================================
 The variable change must take place before structural simplification.
 Subsequently, consider passing `allow_symbolic = true` to `structural_simplify(sys)` to reduce the number of unknowns, with the understanding that the transformation is well-defined.
 
+Usage with non-autonomous systems
+=================================
+If `sys` is 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.
+
 Example
 =======
 Consider a free fall with constant horizontal velocity.
-The laws of physics naturally describes position as a function of time.
-By changing the independent variable, it can be reformulated for vertical position as a function of horizontal distance:
+Physics naturally describes position as a function of time.
+By changing the independent variable, it can be reformulated for vertical position as a function of horizontal position:
 ```julia
 julia> @variables x(t) y(t);
 
@@ -89,118 +94,94 @@ julia> unknowns(M′)
  y(x)
 ```
 """
-function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, verbose = false, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, fold = false, kwargs...)
+    iv2_of_iv1 = unwrap(iv) # e.g. u(t)
+    iv1 = get_iv(sys) # e.g. t
+
     if !iscomplete(sys)
-        error("Cannot change independent variable of incomplete system $(nameof(sys))")
+        error("System $(nameof(sys)) is incomplete. Complete it first!")
     elseif isscheduled(sys)
-        error("Cannot change independent variable of structurally simplified system $(nameof(sys))")
+        error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!")
     elseif !isempty(get_systems(sys))
-        error("Cannot change independent variable of hierarchical system $(nameof(sys)). Flatten it first.") # TODO: implement
+        error("System $(nameof(sys)) is hierarchical. Flatten it first!") # TODO: implement and allow?
+    elseif !iscall(iv2_of_iv1) || !isequal(only(arguments(iv2_of_iv1)), iv1)
+        error("Variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))!")
     end
 
-    iv = unwrap(iv)
-    iv1 = get_iv(sys) # e.g. t
-    if !iscall(iv) || !isequal(only(arguments(iv)), iv1)
-        error("New independent variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))")
-    elseif !isautonomous(sys) && isempty(findall(eq -> isequal(eq.lhs, iv1), eqs))
-        error("System $(nameof(sys)) is autonomous in $iv1. An equation of the form $iv1 ~ F($iv) must be provided.")
+    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)
+    D1 = Differential(iv1) # e.g. d/d(t)
+    D2 = Differential(iv2_of_iv1) # e.g. d/d(u(t))
+
+    # 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)
+        vars = get_variables(ex)
+        for var_of_iv1 in vars # loop through e.g. f(t)
+            if iscall(var_of_iv1) && !isequal(var_of_iv1, iv2_of_iv1) # handle e.g. f(t) -> f(u(t)), but not u(t) -> u(u(t))
+                varname = nameof(operation(var_of_iv1)) # e.g. :f
+                var_of_iv2, = @variables $varname(iv2_of_iv1) # e.g. f(u(t))
+                ex = substitute(ex, var_of_iv1 => var_of_iv2; fold) # e.g. f(t) -> f(u(t))
+            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
 
-    iv2func = iv # e.g. a(t)
-    iv2name = nameof(operation(iv))
-    iv2, = @independent_variables $iv2name # e.g. a
-    D1 = Differential(iv1)
-
-    iv1name = nameof(iv1) # e.g. t
-    iv1func, = @variables $iv1name(iv2) # e.g. t(a)
-
-    eqs = [get_eqs(sys); eqs] # copies system equations to avoid modifying original system
-
-    # 1) Find and compute all necessary expressions for e.g. df/dt, d²f/dt², ...
-    # 1.1) Find the 1st order derivative of the new independent variable (e.g. da(t)/dt = ...), ...
-    div2_div1_idxs = findall(eq -> isequal(eq.lhs, D1(iv2func)), eqs) # index of e.g. da/dt = ...
-    if length(div2_div1_idxs) != 1
-        error("Exactly one equation for $D1($iv2func) was not specified.")
+    # 2) Find e.g. du/dt in equations, then calculate e.g. d²u/dt², ...
+    eqs = [get_eqs(sys); eqs] # 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!")
     end
-    div2_div1_eq = popat!(eqs, only(div2_div1_idxs)) # get and remove e.g. df/dt = ... (may be added back later)
-    div2_div1 = div2_div1_eq.rhs
-    if isequal(div2_div1, 0)
-        error("Cannot change independent variable from $iv1 to $iv2 with singular transformation $div2_div1_eq.")
+    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
-    # 1.2) ... then compute the 2nd order derivative of the new independent variable
-    div1_div2 = 1 / div2_div1 # TODO: URL reference for clarity
-    ddiv2_ddiv1 = expand_derivatives(-Differential(iv2func)(div1_div2) / div1_div2^3, simplify) # e.g. https://math.stackexchange.com/questions/249253/second-derivative-of-the-inverse-function # TODO: higher orders # TODO: pass simplify here
-    # 1.3) # TODO: handle higher orders (3+) derivatives ...
+    ddiv2_of_iv1 = chain_rule(D1(div2_of_iv1)) # TODO: implement higher orders (order >= 3) derivatives with a loop
 
-    # 2) If requested, insert extra dummy equations for e.g. df/dt, d²f/dt², ...
+    # 3) If requested, insert extra dummy equations for e.g. du/dt, d²u/dt², ...
     #    Otherwise, replace all these derivatives by their explicit expressions
     if dummies
-        div2name = Symbol(iv2name, :_t) # TODO: not always t
-        div2, = @variables $div2name(iv2) # e.g. a_t(a)
-        ddiv2name = Symbol(iv2name, :_tt) # TODO: not always t
-        ddiv2, = @variables $ddiv2name(iv2) # e.g. a_tt(a)
-        eqs = [eqs; [div2 ~ div2_div1, ddiv2 ~ ddiv2_ddiv1]] # add dummy equations
-        derivsubs = [D1(D1(iv2func)) => ddiv2, D1(iv2func) => div2] # order is crucial!
+        div2name = Symbol(iv2name, :_, iv1name) # e.g. :u_t # TODO: customize
+        ddiv2name = Symbol(div2name, iv1name) # e.g. :u_tt # TODO: customize
+        div2, ddiv2 = @variables $div2name(iv2) $ddiv2name(iv2) # e.g. u_t(u), u_tt(u)
+        eqs = [eqs; [div2 ~ div2_of_iv1, ddiv2 ~ ddiv2_of_iv1]] # add dummy equations
     else
-        derivsubs = [D1(D1(iv2func)) => ddiv2_ddiv1, D1(iv2func) => div2_div1] # order is crucial!
-    end
-    derivsubs = [derivsubs; [iv2func => iv2, iv1 => iv1func]]
-
-    if verbose
-        # Explain what we just did
-        println("Order 1 (found):    $div2_div1_eq")
-        println("Order 2 (computed): $(D1(div2_div1_eq.lhs) ~ ddiv2_ddiv1)")
-        println()
-        println("Substitutions will be made in this order:")
-        for (n, sub) in enumerate(derivsubs)
-            println("$n: $(sub[1]) => $(sub[2])")
-        end
-        println()
+        div2 = div2_of_iv1
+        ddiv2 = ddiv2_of_iv1
     end
 
-    # 3) Define a transformation function that performs the change of variable on any expression/equation
+    # 4) 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)
+        D1(iv2_of_iv1) => div2 # order 1, e.g. D(u(t)) -> u_t(u)
+        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)
-        verbose && println("Step 0: ", ex)
-
-        # Step 1: substitute f(t₁) => f(t₂(t₁)) in all variables in the expression
-        vars = Symbolics.get_variables(ex)
-        for var1 in vars
-            if Symbolics.iscall(var1) && !isequal(var1, iv2func) # && isequal(only(arguments(var1)), iv1) # skip e.g. constants
-                name = nameof(operation(var1))
-                var2, = @variables $name(iv2func)
-                ex = substitute(ex, var1 => var2; fold = false)
-            end
+        ex = chain_rule(ex)
+        for sub in iv1_to_iv2_subs
+            ex = substitute(ex, sub; fold)
         end
-        verbose && println("Step 1: ", ex)
-
-        # Step 2: expand out all chain rule derivatives
-        ex = expand_derivatives(ex) # expand out with chain rule to get d(iv2)/d(iv1)
-        verbose && println("Step 2: ", ex)
-
-        # Step 3: substitute d²f/dt², df/dt, ... (to dummy variables or explicit expressions, depending on dummies)
-        for sub in derivsubs
-            ex = substitute(ex, sub)
-        end
-        verbose && println("Step 3: ", ex)
-        verbose && println()
-
         return ex
     end
-
-    # 4) Transform all fields
-    eqs = map(transform, eqs)
+    eqs = map(transform, eqs) # we derived and added equations to eqs; they are not in get_eqs(sys)!
     observed = map(transform, get_observed(sys))
     initialization_eqs = map(transform, get_initialization_eqs(sys))
     parameter_dependencies = map(transform, get_parameter_dependencies(sys))
     defaults = Dict(transform(var) => transform(val) for (var, val) in get_defaults(sys))
     guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
     assertions = Dict(transform(condition) => msg for (condition, msg) in get_assertions(sys))
-    # TODO: handle subsystems
 
     # 5) Recreate system with transformed fields
     return typeof(sys)(
         eqs, iv2;
         observed, initialization_eqs, parameter_dependencies, defaults, guesses, assertions,
         name = nameof(sys), description = description(sys), kwargs...
-    ) |> complete # original system had to be complete
+    ) |> complete # input system must be complete, so complete the output system
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 99a252b5b8..cfb3bf5db6 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -59,34 +59,44 @@ end
 end
 
 @testset "Change independent variable (Friedmann equation)" begin
+    @independent_variables t
     D = Differential(t)
-    @variables a(t) ȧ(t) ρr(t) ρm(t) ρΛ(t) ρ(t) P(t) ϕ(t)
-    @parameters Ωr0 Ωm0 ΩΛ0
+    @variables a(t) ȧ(t) Ω(t) ϕ(t)
     eqs = [
-        ρr ~ 3/(8*Num(π)) * Ωr0 / a^4
-        ρm ~ 3/(8*Num(π)) * Ωm0 / a^3
-        ρΛ ~ 3/(8*Num(π)) * ΩΛ0
-        ρ ~ ρr + ρm + ρΛ
-        ȧ ~ √(8*Num(π)/3*ρ*a^4)
+        Ω ~ 123
         D(a) ~ ȧ
+        ȧ ~ √(Ω) * a^2
         D(D(ϕ)) ~ -3*D(a)/a*D(ϕ)
     ]
     M1 = ODESystem(eqs, t; name = :M) |> complete
 
     # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
-    M2 = change_independent_variable(M1, M1.a) |> complete #, D(b) ~ D(a)/a; verbose = true)
+    M2 = change_independent_variable(M1, M1.a; dummies = true)
+    @independent_variables a
+    @variables ȧ(a) Ω(a) ϕ(a) a_t(a) a_tt(a)
+    Da = Differential(a)
+    @test Set(equations(M2)) == Set([
+        a_tt*Da(ϕ) + a_t^2*(Da^2)(ϕ) ~ -3*a_t^2/a*Da(ϕ)
+        ȧ ~ √(Ω) * a^2
+        Ω ~ 123
+        a_t ~ ȧ # 1st order dummy equation
+        a_tt ~ Da(ȧ) * a_t # 2nd order dummy equation
+    ])
+
     @variables b(M2.a)
     M3 = change_independent_variable(M2, b, [Differential(M2.a)(b) ~ exp(-b), M2.a ~ exp(b)])
+
+    M1 = structural_simplify(M1)
     M2 = structural_simplify(M2; allow_symbolic = true)
     M3 = structural_simplify(M3; allow_symbolic = true)
-    @test length(unknowns(M2)) == 2 && length(unknowns(M3)) == 2
+    @test length(unknowns(M3)) == length(unknowns(M2)) == length(unknowns(M1)) - 1
 end
 
 @testset "Change independent variable (simple)" begin
     @variables x(t)
-    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete # TODO: avoid complete. can avoid it if passing defined $variable directly to change_independent_variable
+    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
     Mx = change_independent_variable(Mt, Mt.x; dummies = true)
-    @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2x, x_tt ~ 4x]))
+    @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2*x, x_tt ~ 2*x_t]))
 end
 
 @testset "Change independent variable (free fall)" begin
@@ -101,10 +111,25 @@ end
     @test all(isapprox.(sol[Mx.y], sol[Mx.x - g*(Mx.x/v)^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 (autonomous system)" begin
-    M = ODESystem([D(x) ~ t], t; name = :M) |> complete # non-autonomous
-    @test_throws "t ~ F(x(t)) must be provided" change_independent_variable(M, M.x)
-    @test_nowarn change_independent_variable(M, M.x, [t ~ 2*x])
+@testset "Change independent variable (crazy analytical example)" begin
+    @independent_variables t
+    D = Differential(t)
+    @variables x(t) y(t)
+    M1 = ODESystem([ # crazy non-autonomous non-linear 2nd order ODE
+        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) |> complete
+    M2 = change_independent_variable(M1, M1.x; dummies = true)
+
+    # 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
+    ])
 end
 
 @testset "Change independent variable (errors)" begin

From 34a6a4ee655c68b93904b7cf04ca3e9a3bab28b7 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 5 Mar 2025 12:25:06 +0100
Subject: [PATCH 16/32] Actually run basic transformations test

---
 test/runtests.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 0ae66d3755..8f1f2a0be4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -46,6 +46,7 @@ end
             @safetestset "Model Parsing Test" include("model_parsing.jl")
             @safetestset "Error Handling" include("error_handling.jl")
             @safetestset "StructuralTransformations" include("structural_transformation/runtests.jl")
+            @safetestset "Basic transformations" include("basic_transformations.jl")
             @safetestset "State Selection Test" include("state_selection.jl")
             @safetestset "Symbolic Event Test" include("symbolic_events.jl")
             @safetestset "Stream Connect Test" include("stream_connectors.jl")

From 636ad045cebbf88f6755c9bd4b481ed9ee5c748b Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Wed, 5 Mar 2025 21:41:05 +0100
Subject: [PATCH 17/32] Change independent variable of hierarchical systems

---
 .../tutorials/change_independent_variable.md  | 43 ++++++++-----
 src/systems/diffeqs/basic_transformations.jl  | 63 ++++++++++++-------
 test/basic_transformations.jl                 | 25 +++++---
 3 files changed, 83 insertions(+), 48 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index 0d4d8bcac8..e49ededacd 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -42,7 +42,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, M1.x; dummies = true) |> complete
+M2 = change_independent_variable(M1, M1.x; dummies = true)
 @assert M2.x isa Num # hide
 M2s = structural_simplify(M2; allow_symbolic = true)
 ```
@@ -56,29 +56,38 @@ sol = solve(prob, Tsit5())
 plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
 ```
 
-## 2. Reduce stiffness by changing to a logarithmic time axis
+## 2. Alleviating stiffness by changing to logarithmic time
 
 In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe.
 In terms of conformal time $t$, they can be written
 ```@example changeivar
-@variables a(t) Ω(t) Ωr(t) Ωm(t) ΩΛ(t)
-M1 = ODESystem([
+@variables a(t) Ω(t)
+a = GlobalScope(a) # global var needed by all species
+function species(w; kw...)
+    eqs = [D(Ω) ~ -3(1 + w) * D(a)/a * Ω]
+    return ODESystem(eqs, t, [Ω], []; kw...)
+end
+@named r = species(1//3) # radiation
+@named m = species(0) # matter
+@named Λ = species(-1) # dark energy / cosmological constant
+eqs = [
+    Ω ~ r.Ω + m.Ω + Λ.Ω # total energy density
     D(a) ~ √(Ω) * a^2
-    Ω ~ Ωr + Ωm + ΩΛ
-    D(Ωm) ~ -3 * D(a)/a * Ωm
-    D(Ωr) ~ -4 * D(a)/a * Ωr
-    D(ΩΛ) ~ 0
-], t; initialization_eqs = [
-    ΩΛ + Ωr + Ωm ~ 1
-], name = :M) |> complete
+]
+initialization_eqs = [
+    Λ.Ω + r.Ω + m.Ω ~ 1
+]
+M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M)
+M1 = compose(M1, r, m, Λ)
+M1 = complete(M1; flatten = false)
 M1s = structural_simplify(M1)
 ```
 Of course, we can solve this ODE as it is:
 ```@example changeivar
-prob = ODEProblem(M1s, [M1.a => 1.0, M1.Ωr => 5e-5, M1.Ωm => 0.3], (0.0, -3.3), [])
+prob = ODEProblem(M1s, [M1.a => 1.0, M1.r.Ω => 5e-5, M1.m.Ω => 0.3], (0.0, -3.3), [])
 sol = solve(prob, Tsit5(); reltol = 1e-5)
 @assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide
-plot(sol, idxs = [M1.a, M1.Ωr/M1.Ω, M1.Ωm/M1.Ω, M1.ΩΛ/M1.Ω])
+plot(sol, idxs = [M1.a, M1.r.Ω/M1.Ω, M1.m.Ω/M1.Ω, M1.Λ.Ω/M1.Ω])
 ```
 The solver becomes unstable due to stiffness.
 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.
@@ -89,7 +98,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) |> complete
+M2 = change_independent_variable(M1, M1.a; dummies = true)
 ```
 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$:
@@ -97,14 +106,14 @@ This means that to change the independent variable from $a$ to $b$, we must prov
 a = M2.a
 Da = Differential(a)
 @variables b(a)
-M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)]) |> complete
+M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)])
 ```
 We can now solve and plot the ODE in terms of $b$:
 ```@example changeivar
 M3s = structural_simplify(M3; allow_symbolic = true)
-prob = ODEProblem(M3s, [M3.Ωr => 5e-5, M3.Ωm => 0.3], (0, -15), [])
+prob = ODEProblem(M3s, [M3.r.Ω => 5e-5, M3.m.Ω => 0.3], (0, -15), [])
 sol = solve(prob, Tsit5())
 @assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide
-plot(sol, idxs = [M3.Ωr/M3.Ω, M3.Ωm/M3.Ω, M3.ΩΛ/M3.Ω])
+plot(sol, idxs = [M3.r.Ω/M3.Ω, M3.m.Ω/M3.Ω, M3.Λ.Ω/M3.Ω])
 ```
 Notice that the variables vary "more nicely" with respect to $b$ than $t$, making the solver happier and avoiding numerical problems.
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 11a5514f42..a8f3db11e4 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -77,6 +77,11 @@ Usage with non-autonomous systems
 If `sys` is 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.
 
+Usage with hierarchical systems
+===============================
+It is recommended that `iv` is a non-namespaced variable in `sys`.
+This means it can belong to the top-level system or be a variable in a subsystem declared with `GlobalScope`.
+
 Example
 =======
 Consider a free fall with constant horizontal velocity.
@@ -102,8 +107,6 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
         error("System $(nameof(sys)) is incomplete. Complete it first!")
     elseif isscheduled(sys)
         error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!")
-    elseif !isempty(get_systems(sys))
-        error("System $(nameof(sys)) is hierarchical. Flatten it first!") # TODO: implement and allow?
     elseif !iscall(iv2_of_iv1) || !isequal(only(arguments(iv2_of_iv1)), iv1)
         error("Variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))!")
     end
@@ -112,6 +115,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     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
     D1 = Differential(iv1) # e.g. d/d(t)
     D2 = Differential(iv2_of_iv1) # e.g. d/d(u(t))
 
@@ -119,10 +123,9 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     function chain_rule(ex)
         vars = get_variables(ex)
         for var_of_iv1 in vars # loop through e.g. f(t)
-            if iscall(var_of_iv1) && !isequal(var_of_iv1, iv2_of_iv1) # handle e.g. f(t) -> f(u(t)), but not u(t) -> u(u(t))
-                varname = nameof(operation(var_of_iv1)) # e.g. :f
-                var_of_iv2, = @variables $varname(iv2_of_iv1) # e.g. f(u(t))
-                ex = substitute(ex, var_of_iv1 => var_of_iv2; fold) # e.g. f(t) -> f(u(t))
+            if iscall(var_of_iv1) && isequal(only(arguments(var_of_iv1)), iv1) && !isequal(var_of_iv1, iv2_of_iv1) # handle e.g. f(t) -> f(u(t)), but not u(t) -> u(u(t))
+                var_of_iv2 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(t) -> f(u(t))
+                ex = substitute(ex, var_of_iv1 => var_of_iv2)
             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
@@ -130,7 +133,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     end
 
     # 2) Find e.g. du/dt in equations, then calculate e.g. d²u/dt², ...
-    eqs = [get_eqs(sys); eqs] # all equations (system-defined + user-provided) we may use
+    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!")
@@ -148,11 +151,14 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
         div2name = Symbol(iv2name, :_, iv1name) # e.g. :u_t # TODO: customize
         ddiv2name = Symbol(div2name, iv1name) # e.g. :u_tt # TODO: customize
         div2, ddiv2 = @variables $div2name(iv2) $ddiv2name(iv2) # e.g. u_t(u), u_tt(u)
-        eqs = [eqs; [div2 ~ div2_of_iv1, ddiv2 ~ ddiv2_of_iv1]] # add dummy equations
+        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
+        @set! sys.unknowns = [get_unknowns(sys); [div2, ddiv2]] # add dummy variables
     else
         div2 = div2_of_iv1
         ddiv2 = ddiv2_of_iv1
     end
+    @set! sys.eqs = eqs # add extra equations we derived before starting transformation process
 
     # 4) 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).
@@ -170,18 +176,31 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
         end
         return ex
     end
-    eqs = map(transform, eqs) # we derived and added equations to eqs; they are not in get_eqs(sys)!
-    observed = map(transform, get_observed(sys))
-    initialization_eqs = map(transform, get_initialization_eqs(sys))
-    parameter_dependencies = map(transform, get_parameter_dependencies(sys))
-    defaults = Dict(transform(var) => transform(val) for (var, val) in get_defaults(sys))
-    guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
-    assertions = Dict(transform(condition) => msg for (condition, msg) in get_assertions(sys))
-
-    # 5) Recreate system with transformed fields
-    return typeof(sys)(
-        eqs, iv2;
-        observed, initialization_eqs, parameter_dependencies, defaults, guesses, assertions,
-        name = nameof(sys), description = description(sys), kwargs...
-    ) |> complete # input system must be complete, so complete the output system
+    function transform(sys::AbstractODESystem)
+        eqs = map(transform, get_eqs(sys))
+        unknowns = map(transform, get_unknowns(sys))
+        unknowns = filter(var -> !isequal(var, iv2), unknowns) # remove e.g. u
+        ps = map(transform, get_ps(sys))
+        ps = filter(!isinitial, ps) # remove Initial(...) # # TODO: shouldn't have to touch this
+        observed = map(transform, get_observed(sys))
+        initialization_eqs = map(transform, get_initialization_eqs(sys))
+        parameter_dependencies = map(transform, get_parameter_dependencies(sys))
+        defaults = Dict(transform(var) => transform(val) for (var, val) in get_defaults(sys))
+        guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
+        assertions = Dict(transform(condition) => msg for (condition, msg) in get_assertions(sys))
+        systems = get_systems(sys) # save before reconstructing system
+        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...
+        )
+        systems = map(transform, systems) # recurse through subsystems
+        sys = compose(sys, systems) # rebuild hierarchical system
+        if wascomplete
+            wasflat = isempty(systems)
+            sys = complete(sys; flatten = wasflat) # complete output if input was complete
+        end
+        return sys
+    end
+    return transform(sys)
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index cfb3bf5db6..4fed453390 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -62,25 +62,34 @@ end
     @independent_variables t
     D = Differential(t)
     @variables a(t) ȧ(t) Ω(t) ϕ(t)
+    a, ȧ = GlobalScope.([a, ȧ])
+    species(w; kw...) = ODESystem([D(Ω) ~ -3(1 + w) * D(a)/a * Ω], t, [Ω], []; kw...)
+    @named r = species(1//3)
+    @named m = species(0)
+    @named Λ = species(-1)
     eqs = [
-        Ω ~ 123
+        Ω ~ r.Ω + m.Ω + Λ.Ω
         D(a) ~ ȧ
         ȧ ~ √(Ω) * a^2
         D(D(ϕ)) ~ -3*D(a)/a*D(ϕ)
     ]
-    M1 = ODESystem(eqs, t; name = :M) |> complete
+    M1 = ODESystem(eqs, t, [Ω, a, ȧ, ϕ], []; name = :M)
+    M1 = compose(M1, r, m, Λ)
+    M1 = complete(M1; flatten = false)
 
     # 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)
-    @independent_variables a
-    @variables ȧ(a) Ω(a) ϕ(a) a_t(a) a_tt(a)
+    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, a_t, a_tt = M2.a, M2.ȧ, M2.Ω, M2.r.Ω, M2.m.Ω, M2.Λ.Ω, M2.ϕ, M2.a_t, M2.a_tt
     Da = Differential(a)
     @test Set(equations(M2)) == Set([
-        a_tt*Da(ϕ) + a_t^2*(Da^2)(ϕ) ~ -3*a_t^2/a*Da(ϕ)
-        ȧ ~ √(Ω) * a^2
-        Ω ~ 123
         a_t ~ ȧ # 1st order dummy equation
         a_tt ~ Da(ȧ) * a_t # 2nd order dummy equation
+        Ω ~ Ωr + Ωm + ΩΛ
+        ȧ ~ √(Ω) * a^2
+        a_tt*Da(ϕ) + 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
     ])
 
     @variables b(M2.a)
@@ -143,6 +152,4 @@ end
     @test_throws "not specified" change_independent_variable(M, v)
     @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 = compose(M, M)
-    @test_throws "hierarchical" change_independent_variable(M, M.x)
 end

From 2afacb5c696c2e760fd62597b427fae21e19a560 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 18:34:25 +0100
Subject: [PATCH 18/32] Forbid DDEs for now

---
 src/systems/diffeqs/basic_transformations.jl | 2 ++
 test/basic_transformations.jl                | 4 ++++
 2 files changed, 6 insertions(+)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index a8f3db11e4..0df7a5d51d 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -105,6 +105,8 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
 
     if !iscomplete(sys)
         error("System $(nameof(sys)) is incomplete. Complete it first!")
+    elseif is_dde(sys)
+        error("System $(nameof(sys)) contains delay differential equations (DDEs). This is currently not supported!")
     elseif isscheduled(sys)
         error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!")
     elseif !iscall(iv2_of_iv1) || !isequal(only(arguments(iv2_of_iv1)), iv1)
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 4fed453390..798548ff44 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -152,4 +152,8 @@ end
     @test_throws "not specified" change_independent_variable(M, v)
     @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) |> complete
+    @test_throws "DDE" change_independent_variable(M, M.x)
 end

From 2d5aa12109c8b494e5f4af975e09a63630e14fce Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 18:48:29 +0100
Subject: [PATCH 19/32] Use vars(ex; op = Nothing) instead of get_variables(ex)

---
 src/systems/diffeqs/basic_transformations.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 0df7a5d51d..98c50758c2 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -123,8 +123,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
 
     # 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)
-        vars = get_variables(ex)
-        for var_of_iv1 in vars # loop through e.g. f(t)
+        for var_of_iv1 in vars(ex; op = Nothing) # loop over all variables in expression, e.g. f(t) (op = Nothing so "D(f(t))" yields only "f(t)" and not "D(f(t))"; we want to replace *inside* derivative signs)
             if iscall(var_of_iv1) && isequal(only(arguments(var_of_iv1)), iv1) && !isequal(var_of_iv1, iv2_of_iv1) # handle e.g. f(t) -> f(u(t)), but not u(t) -> u(u(t))
                 var_of_iv2 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # e.g. f(t) -> f(u(t))
                 ex = substitute(ex, var_of_iv1 => var_of_iv2)

From 5780d91fd26f5fdcdc48980bc22b0a318284a2f1 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 19:04:34 +0100
Subject: [PATCH 20/32] Print detected transformation equations when there is
 not exactly one

---
 src/systems/diffeqs/basic_transformations.jl | 2 +-
 test/basic_transformations.jl                | 6 ++++--
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 98c50758c2..3d9d92b13f 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -137,7 +137,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     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!")
+        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)
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 798548ff44..b1d66fc9ee 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -148,8 +148,10 @@ end
     M = complete(M)
     @test_throws "singular" change_independent_variable(M, M.x)
     @test_throws "structurally simplified" change_independent_variable(structural_simplify(M), y)
-    @test_throws "not specified" change_independent_variable(M, w)
-    @test_throws "not specified" change_independent_variable(M, v)
+    @test_throws "Got 0 equations:" change_independent_variable(M, w)
+    @test_throws "Got 0 equations:" change_independent_variable(M, v)
+    M = ODESystem([D(x) ~ 1, v ~ 1], t; name = :M) |> complete
+    @test_throws "Got 2 equations:" change_independent_variable(M, 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)

From 5ca058d64b1918c40ec910acdf652e8263f537d6 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 19:38:33 +0100
Subject: [PATCH 21/32] Use default_toterm (with special underscore) for
 dummies

---
 src/systems/diffeqs/basic_transformations.jl |  5 ++--
 test/basic_transformations.jl                | 24 ++++++++++----------
 2 files changed, 14 insertions(+), 15 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 3d9d92b13f..7ed3c16049 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -149,9 +149,8 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     # 3) If requested, insert extra dummy equations for e.g. du/dt, d²u/dt², ...
     #    Otherwise, replace all these derivatives by their explicit expressions
     if dummies
-        div2name = Symbol(iv2name, :_, iv1name) # e.g. :u_t # TODO: customize
-        ddiv2name = Symbol(div2name, iv1name) # e.g. :u_tt # TODO: customize
-        div2, ddiv2 = @variables $div2name(iv2) $ddiv2name(iv2) # e.g. u_t(u), u_tt(u)
+        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
         @set! sys.unknowns = [get_unknowns(sys); [div2, ddiv2]] # add dummy variables
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index b1d66fc9ee..7a8c063283 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -79,17 +79,17 @@ end
 
     # 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)
-    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, a_t, a_tt = M2.a, M2.ȧ, M2.Ω, M2.r.Ω, M2.m.Ω, M2.Λ.Ω, M2.ϕ, M2.a_t, M2.a_tt
+    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, aˍt, aˍtt = M2.a, M2.ȧ, M2.Ω, M2.r.Ω, M2.m.Ω, M2.Λ.Ω, M2.ϕ, M2.aˍt, M2.aˍtt
     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 ~ ȧ # 1st order dummy equation
+        aˍtt ~ Da(ȧ) * aˍt # 2nd order dummy equation
         Ω ~ Ωr + Ωm + ΩΛ
         ȧ ~ √(Ω) * a^2
-        a_tt*Da(ϕ) + 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
+        aˍtt*Da(ϕ) + 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
     ])
 
     @variables b(M2.a)
@@ -105,7 +105,7 @@ end
     @variables x(t)
     Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
     Mx = change_independent_variable(Mt, Mt.x; dummies = true)
-    @test (@variables x x_t(x) x_tt(x); Set(equations(Mx)) == Set([x_t ~ 2*x, x_tt ~ 2*x_t]))
+    @test (@variables x xˍt(x) xˍtt(x); Set(equations(Mx)) == Set([xˍt ~ 2*x, xˍtt ~ 2*xˍt]))
 end
 
 @testset "Change independent variable (free fall)" begin
@@ -133,11 +133,11 @@ end
     # Compare to pen-and-paper result
     @independent_variables x
     Dx = Differential(x)
-    @variables x_t(x) x_tt(x) y(x) t(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ˍ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
     ])
 end
 

From 2397d9a22cda11008f0288d50e1b23a4d24b80b5 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 20:18:27 +0100
Subject: [PATCH 22/32] Change independent variable of incomplete systems

---
 .../tutorials/change_independent_variable.md  | 11 +++---
 src/systems/diffeqs/basic_transformations.jl  |  4 +--
 test/basic_transformations.jl                 | 36 +++++++++----------
 3 files changed, 23 insertions(+), 28 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index e49ededacd..f196c3cb64 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -24,11 +24,11 @@ D = Differential(t)
 M1 = ODESystem([
     D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
 ], t; defaults = [
-    y => 0.0
+    y => 0.0 # start on the ground
 ], initialization_eqs = [
     #x ~ 0.0, # TODO: handle? # hide
     D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
-], name = :M) |> complete
+], name = :M)
 M1s = structural_simplify(M1)
 ```
 This is the standard parametrization that arises naturally from kinematics and Newton's laws.
@@ -42,7 +42,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, M1.x; dummies = true)
+M2 = change_independent_variable(M1, x; dummies = true)
 @assert M2.x isa Num # hide
 M2s = structural_simplify(M2; allow_symbolic = true)
 ```
@@ -79,12 +79,11 @@ initialization_eqs = [
 ]
 M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M)
 M1 = compose(M1, r, m, Λ)
-M1 = complete(M1; flatten = false)
 M1s = structural_simplify(M1)
 ```
 Of course, we can solve this ODE as it is:
 ```@example changeivar
-prob = ODEProblem(M1s, [M1.a => 1.0, M1.r.Ω => 5e-5, M1.m.Ω => 0.3], (0.0, -3.3), [])
+prob = ODEProblem(M1s, [M1s.a => 1.0, M1s.r.Ω => 5e-5, M1s.m.Ω => 0.3], (0.0, -3.3), [])
 sol = solve(prob, Tsit5(); reltol = 1e-5)
 @assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide
 plot(sol, idxs = [M1.a, M1.r.Ω/M1.Ω, M1.m.Ω/M1.Ω, M1.Λ.Ω/M1.Ω])
@@ -111,7 +110,7 @@ M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)])
 We can now solve and plot the ODE in terms of $b$:
 ```@example changeivar
 M3s = structural_simplify(M3; allow_symbolic = true)
-prob = ODEProblem(M3s, [M3.r.Ω => 5e-5, M3.m.Ω => 0.3], (0, -15), [])
+prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), [])
 sol = solve(prob, Tsit5())
 @assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide
 plot(sol, idxs = [M3.r.Ω/M3.Ω, M3.m.Ω/M3.Ω, M3.Λ.Ω/M3.Ω])
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 7ed3c16049..a6ce1fcc14 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -103,9 +103,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
     iv2_of_iv1 = unwrap(iv) # e.g. u(t)
     iv1 = get_iv(sys) # e.g. t
 
-    if !iscomplete(sys)
-        error("System $(nameof(sys)) is incomplete. Complete it first!")
-    elseif is_dde(sys)
+    if is_dde(sys)
         error("System $(nameof(sys)) contains delay differential equations (DDEs). This is currently not supported!")
     elseif isscheduled(sys)
         error("System $(nameof(sys)) is structurally simplified. Change independent variable before structural simplification!")
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 7a8c063283..0aa679a817 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -28,8 +28,8 @@ end
 @testset "Change independent variable (trivial)" begin
     @variables x(t) y(t)
     eqs1 = [D(D(x)) ~ D(x) + x, D(y) ~ 1]
-    M1 = ODESystem(eqs1, t; name = :M) |> complete
-    M2 = change_independent_variable(M1, M1.y)
+    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
 end
@@ -43,8 +43,8 @@ end
         D(s) ~ 1 / (2*s)
     ]
     initialization_eqs = [x ~ 1.0, y ~ 1.0, D(y) ~ 0.0]
-    M1 = ODESystem(eqs, t; initialization_eqs, name = :M) |> complete
-    M2 = change_independent_variable(M1, M1.s; dummies = true)
+    M1 = ODESystem(eqs, t; initialization_eqs, name = :M)
+    M2 = change_independent_variable(M1, s; dummies = true)
 
     M1 = structural_simplify(M1; allow_symbolic = true)
     M2 = structural_simplify(M2; allow_symbolic = true)
@@ -75,11 +75,11 @@ end
     ]
     M1 = ODESystem(eqs, t, [Ω, a, ȧ, ϕ], []; name = :M)
     M1 = compose(M1, r, m, Λ)
-    M1 = complete(M1; flatten = false)
 
     # 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)
-    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, aˍt, aˍtt = M2.a, M2.ȧ, M2.Ω, M2.r.Ω, M2.m.Ω, M2.Λ.Ω, M2.ϕ, M2.aˍt, M2.aˍtt
+    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
     Da = Differential(a)
     @test Set(equations(M2)) == Set([
         aˍt ~ ȧ # 1st order dummy equation
@@ -103,16 +103,16 @@ end
 
 @testset "Change independent variable (simple)" begin
     @variables x(t)
-    Mt = ODESystem([D(x) ~ 2*x], t; name = :M) |> complete
-    Mx = change_independent_variable(Mt, Mt.x; dummies = true)
+    Mt = ODESystem([D(x) ~ 2*x], t; name = :M)
+    Mx = change_independent_variable(Mt, x; dummies = true)
     @test (@variables x xˍt(x) xˍtt(x); Set(equations(Mx)) == Set([xˍt ~ 2*x, xˍtt ~ 2*xˍt]))
 end
 
 @testset "Change independent variable (free fall)" 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) |> complete # gives (x, y) as function of t, ...
-    Mx = change_independent_variable(Mt, Mt.x; dummies = false) # ... but we want y as a function of x
+    Mt = ODESystem([D(D(y)) ~ -g, D(x) ~ v], t; name = :M) # gives (x, y) as function of t, ...
+    Mx = change_independent_variable(Mt, x; dummies = false) # ... 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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities
@@ -127,8 +127,8 @@ end
     M1 = ODESystem([ # crazy non-autonomous non-linear 2nd order ODE
         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) |> complete
-    M2 = change_independent_variable(M1, M1.x; dummies = true)
+    ], t; name = :M)
+    M2 = change_independent_variable(M1, x; dummies = true)
 
     # Compare to pen-and-paper result
     @independent_variables x
@@ -144,18 +144,16 @@ 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 "incomplete" change_independent_variable(M, M.x)
-    M = complete(M)
-    @test_throws "singular" change_independent_variable(M, M.x)
+    @test_throws "singular" change_independent_variable(M, x)
     @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([D(x) ~ 1, v ~ 1], t; name = :M) |> complete
-    @test_throws "Got 2 equations:" change_independent_variable(M, M.x, [D(x) ~ 2])
+    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) |> complete
-    @test_throws "DDE" change_independent_variable(M, M.x)
+    M = ODESystem([D(x(t)) ~ x(t-1)], t; name = :M)
+    @test_throws "DDE" change_independent_variable(M, x(t))
 end

From a27785418814fe917fdeb096943375e92349ac79 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 20:38:01 +0100
Subject: [PATCH 23/32] Warn user about subtle differences between variables
 after change of independent variable in documentation

---
 .../tutorials/change_independent_variable.md  | 22 +++++++++++++------
 1 file changed, 15 insertions(+), 7 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index f196c3cb64..b308bd3f82 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -23,11 +23,8 @@ D = Differential(t)
 @parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
 M1 = ODESystem([
     D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
-], t; defaults = [
-    y => 0.0 # start on the ground
-], initialization_eqs = [
-    #x ~ 0.0, # TODO: handle? # hide
-    D(x) ~ D(y) # equal initial horizontal and vertical velocity (45 °)
+], t; initialization_eqs = [
+    D(x) ~ D(y) # equal initial horizontal and vertical velocity (45°)
 ], name = :M)
 M1s = structural_simplify(M1)
 ```
@@ -43,15 +40,26 @@ We will demonstrate the last method by changing the independent variable from $t
 This transformation is well-defined for any non-zero horizontal velocity $v$.
 ```@example changeivar
 M2 = change_independent_variable(M1, x; dummies = true)
-@assert M2.x isa Num # hide
 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
+@assert allequal([M2.x, M2s.x]) # hide
+@assert allequal([y, M1s.y]) # hide
+@assert allunique([M1.x, M1.y, M2.y, M2s.y]) # hide
+M2s # display this # hide
 ```
 The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`.
+
+!!! warn
+    At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables.
+    Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables.
+    It can be instructive to inspect these yourself to see their subtle differences.
+
 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.
 It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
 ```@example changeivar
 using OrdinaryDiffEq, Plots
-prob = ODEProblem(M2s, [], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
+prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
 sol = solve(prob, Tsit5())
 plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
 ```

From aa29107b95d0fdf109df73b67cf49a07125eb471 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 20:51:00 +0100
Subject: [PATCH 24/32] Update change_independent_variable docstring

---
 src/systems/diffeqs/basic_transformations.jl | 28 ++++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index a6ce1fcc14..60159023f7 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -60,30 +60,30 @@ 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``).
 
-Keyword arguments
-=================
-If `dummies`, derivatives of the new independent variable with respect to the old one are expressed through dummy equations; otherwise they are explicitly inserted into the equations.
-If `simplify`, these dummy expressions are simplified and often give a tidier transformation.
-If `fold`, internal substitutions will evaluate numerical expressions.
+# 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.
+- `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
-======================================
+# Usage before structural simplification
+
 The variable change must take place before structural simplification.
 Subsequently, consider passing `allow_symbolic = true` to `structural_simplify(sys)` to reduce the number of unknowns, with the understanding that the transformation is well-defined.
 
-Usage with non-autonomous systems
-=================================
-If `sys` is 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))``).
+# Usage with non-autonomous systems
+
+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.
 
-Usage with hierarchical systems
-===============================
+# Usage with hierarchical systems
+
 It is recommended that `iv` is a non-namespaced variable in `sys`.
 This means it can belong to the top-level system or be a variable in a subsystem declared with `GlobalScope`.
 
-Example
-=======
+# Example
+
 Consider a free fall with constant horizontal velocity.
 Physics naturally describes position as a function of time.
 By changing the independent variable, it can be reformulated for vertical position as a function of horizontal position:

From ea5a0032165e56a2cf61985d664031e7b3409ce6 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Thu, 6 Mar 2025 20:57:47 +0100
Subject: [PATCH 25/32] Explicitly test that c*D(x) ~ something fails, but add
 TODO to fix it

---
 test/basic_transformations.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 0aa679a817..ea781536ef 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -148,6 +148,8 @@ end
     @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)

From 0e41e045616cc260ce61b430607ef33ccf44eedb Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Fri, 7 Mar 2025 19:30:05 +0100
Subject: [PATCH 26/32] Prepare test for array variables (until
 expand_derivatives bug is fixed)

---
 src/systems/diffeqs/basic_transformations.jl | 8 +++++---
 test/basic_transformations.jl                | 8 +++++---
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 60159023f7..0806a105bb 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -121,9 +121,11 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
 
     # 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)
-        for var_of_iv1 in vars(ex; op = Nothing) # loop over all variables in expression, e.g. f(t) (op = Nothing so "D(f(t))" yields only "f(t)" and not "D(f(t))"; we want to replace *inside* derivative signs)
-            if iscall(var_of_iv1) && isequal(only(arguments(var_of_iv1)), iv1) && !isequal(var_of_iv1, iv2_of_iv1) # handle e.g. f(t) -> f(u(t)), but not u(t) -> u(u(t))
-                var_of_iv2 = substitute(var_of_iv1, iv1 => iv2_of_iv1) # 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))
+                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)
             end
         end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index ea781536ef..cfbcea01b1 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -102,10 +102,12 @@ end
 end
 
 @testset "Change independent variable (simple)" begin
-    @variables x(t)
-    Mt = ODESystem([D(x) ~ 2*x], t; name = :M)
+    @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)
-    @test (@variables x xˍt(x) xˍtt(x); Set(equations(Mx)) == Set([xˍt ~ 2*x, xˍtt ~ 2*xˍt]))
+    @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]))
 end
 
 @testset "Change independent variable (free fall)" begin

From d66026947da0f471b8028348effce85e06bcc630 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Fri, 7 Mar 2025 20:22:49 +0100
Subject: [PATCH 27/32] Test change_independent_variable with registered
 functions and callable parameters

---
 test/basic_transformations.jl | 30 +++++++++++++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index cfbcea01b1..3853debcd4 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, Test
+using ModelingToolkit, OrdinaryDiffEq, DataInterpolations, Test
 
 @independent_variables t
 D = Differential(t)
@@ -143,6 +143,34 @@ end
     ])
 end
 
+@testset "Change independent variable (registered function / callable parameter)" begin
+    @independent_variables t
+    D = Differential(t)
+    @variables x(t) y(t)
+    @parameters f::LinearInterpolation (fc::LinearInterpolation)(..) # non-callable and callable
+    callme(interp::LinearInterpolation, input) = interp(input)
+    @register_symbolic callme(interp::LinearInterpolation, input)
+    M1 = ODESystem([
+        D(x) ~ 2*t
+        D(y) ~ 1*fc(t) + 2*fc(x) + 3*fc(y) + 1*callme(f, t) + 2*callme(f, x) + 3*callme(f, y)
+    ], t; name = :M)
+
+    # Ensure that interpolations are called with the same variables
+    M2 = change_independent_variable(M1, x, [t ~ √(x)])
+    @variables 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)
+    ])
+
+    _f = LinearInterpolation([1.0, 1.0], [-100.0, +100.0]) # constant value 1
+    M2s = structural_simplify(M2; allow_symbolic = true)
+    prob = ODEProblem(M2s, [M2s.y => 0.0], (1.0, 4.0), [fc => _f, f => _f])
+    sol = solve(prob, Tsit5(); abstol = 1e-5)
+    @test isapprox(sol(4.0, idxs=M2.y), 12.0; atol = 1e-5) # Anal solution is D(y) ~ 12 => y(t) ~ 12*t + C => y(x) ~ 12*√(x) + C. With y(x=1)=0 => 12*(√(x)-1), so y(x=4) ~ 12
+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)

From 10793d28f28b5ce8bad2b3f3ffe2dc439f41e845 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Fri, 7 Mar 2025 21:02:02 +0100
Subject: [PATCH 28/32] Optionally add differential equation for old
 independent variable

---
 .../tutorials/change_independent_variable.md  |  5 ++++
 src/systems/diffeqs/basic_transformations.jl  | 23 ++++++++++++++-----
 test/basic_transformations.jl                 |  6 ++---
 3 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index b308bd3f82..35a065e706 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -64,6 +64,11 @@ sol = solve(prob, Tsit5())
 plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
 ```
 
+!!! tip "Usage tips"
+    Look up the documentation of [`change_independent_variable`](@ref) for tips on how to use it.
+
+    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.
+
 ## 2. Alleviating stiffness by changing to logarithmic time
 
 In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe.
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 0806a105bb..5e3c09d1d4 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -51,7 +51,7 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
 end
 
 """
-    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, fold = false, kwargs...)
+    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, 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``).
@@ -63,6 +63,7 @@ This is satisfied if one is a strictly increasing function of the other (e.g. ``
 # 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.
 - `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`.
@@ -99,7 +100,7 @@ julia> unknowns(M′)
  y(x)
 ```
 """
-function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, simplify = true, fold = false, kwargs...)
+function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummies = false, 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
 
@@ -148,24 +149,34 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; dummi
 
     # 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
-        @set! sys.unknowns = [get_unknowns(sys); [div2, ddiv2]] # add dummy variables
+        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
 
-    # 4) Transform everything from old to new independent variable, e.g. t -> u.
+    # 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)
-        D1(iv2_of_iv1) => div2 # order 1, e.g. D(u(t)) -> u_t(u)
+        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)
     ]
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 3853debcd4..ac3720469a 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -114,12 +114,12 @@ end
     @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, ...
-    Mx = change_independent_variable(Mt, x; dummies = false) # ... but we want y as a function of x
+    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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities
+    prob = ODEProblem(Mx, [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0], (0.0, 20.0), [g => 9.81, v => 10.0]) # 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.x/v)^2/2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2)
+    @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

From 2b0998b74409f3cf04aebeae2652a113fd1ec5bd Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 9 Mar 2025 00:47:27 +0100
Subject: [PATCH 29/32] Rewrite change_independent_variable to handle any
 derivative order and nonlinear equations

---
 .../tutorials/change_independent_variable.md  |   4 +-
 src/systems/diffeqs/basic_transformations.jl  | 120 +++++++-----------
 test/basic_transformations.jl                 |  56 ++++----
 3 files changed, 77 insertions(+), 103 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index 35a065e706..8f9ee0cbd2 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/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$:
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index 5e3c09d1d4..c05e5e4394 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ 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
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index ac3720469a..8bf60ed68a 100644
--- a/test/basic_transformations.jl
+++ b/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))

From f2a1a5ae580487e3725ca85fbd9e1083719a1b48 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 9 Mar 2025 12:25:15 +0100
Subject: [PATCH 30/32] Test 3rd order nonlinear system

---
 test/basic_transformations.jl | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index 8bf60ed68a..dd6d677b60 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -134,24 +134,24 @@ 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 (crazy analytical example)" begin
+@testset "Change independent variable (crazy 3rd order nonlinear system)" begin
     @independent_variables t
     D = Differential(t)
     @variables x(t) y(t)
-    M1 = ODESystem([ # crazy non-autonomous non-linear 2nd order ODE
-        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
+    M1 = ODESystem([
+        (D^3)(y) ~ D(x)^2 + (D^2)(y^2) |> expand_derivatives
+        D(x)^2 + D(y)^2 ~ x^4 + y^5 + t^6
     ], t; name = :M)
-    M2 = change_independent_variable(M1, x)
+    M2 = change_independent_variable(M1, x; add_old_diff = true)
+    @test_nowarn structural_simplify(M2)
 
     # Compare to pen-and-paper result
-    @independent_variables x
+    @variables x xˍt(x) xˍt(x) y(x) t(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ˍ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
-    ])
+    areequivalent(eq1, eq2) = isequal(expand(eq1.lhs - eq2.lhs), 0) && isequal(expand(eq1.rhs - eq2.rhs), 0)
+    @test areequivalent(equations(M2)[1], xˍt^3*(Dx^3)(y) + xˍt^2*Dx(y)*(Dx^2)(xˍt) + xˍt*Dx(y)*(Dx(xˍt))^2 + 3*xˍt^2*(Dx^2)(y)*Dx(xˍt) ~ xˍt^2 + 2*xˍt^2*Dx(y)^2 + 2*xˍt^2*y*(Dx^2)(y) + 2*y*Dx(y)*Dx(xˍt)*xˍt)
+    @test areequivalent(equations(M2)[2], xˍt^2 + xˍt^2*Dx(y)^2 ~ x^4 + y^5 + t^6)
+    @test areequivalent(equations(M2)[3], Dx(t) ~ 1 / xˍt)
 end
 
 @testset "Change independent variable (registered function / callable parameter)" begin

From aaae59df56fa8a200317420fba65579d56e12ecc Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Sun, 9 Mar 2025 12:33:09 +0100
Subject: [PATCH 31/32] Clean up and format

---
 .../tutorials/change_independent_variable.md  |  78 +++++++-----
 src/ModelingToolkit.jl                        |   3 +-
 src/systems/diffeqs/basic_transformations.jl  |  49 ++++---
 test/basic_transformations.jl                 | 120 +++++++++++-------
 4 files changed, 150 insertions(+), 100 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index 8f9ee0cbd2..e1243fb2ee 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -2,11 +2,11 @@
 
 Ordinary differential equations describe the rate of change of some dependent variables with respect to one independent variable.
 For the modeler it is often most natural to write down the equations with a particular independent variable, say time $t$.
-In many cases there are good reasons for reparametrizing ODEs in terms of a different independent variable:
+However, in many cases there are good reasons for changing the independent variable:
 
-1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
-2. Some differential equations vary more nicely (e.g. less stiff or better behaved) with respect to one independent variable than another.
-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).
+    1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
+    2. Some differential equations vary more nicely (e.g. less stiff) with respect to one independent variable than another.
+    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).
 
 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.
 This is mechanical and error-prone.
@@ -14,30 +14,34 @@ ModelingToolkit provides the utility function [`change_independent_variable`](@r
 
 ## 1. Get one dependent variable as a function of another
 
-Consider a projectile shot with some initial velocity in a gravitational field.
+Consider a projectile shot with some initial velocity in a vertical gravitational field with constant horizontal velocity.
+
 ```@example changeivar
 using ModelingToolkit
 @independent_variables t
 D = Differential(t)
 @variables x(t) y(t)
-@parameters g = 9.81 v # gravitational acceleration and constant horizontal velocity
-M1 = ODESystem([
-    D(D(y)) ~ -g, D(x) ~ v # constant horizontal velocity
-], t; initialization_eqs = [
-    D(x) ~ D(y) # equal initial horizontal and vertical velocity (45°)
-], name = :M)
+@parameters g=9.81 v # gravitational acceleration and horizontal velocity
+eqs = [D(D(y)) ~ -g, D(x) ~ v]
+initialization_eqs = [D(x) ~ D(y)] # 45° initial angle
+M1 = ODESystem(eqs, t; initialization_eqs, name = :M)
 M1s = structural_simplify(M1)
+@assert length(equations(M1s)) == 3 # hide
+M1s # hide
 ```
+
 This is the standard parametrization that arises naturally from kinematics and Newton's laws.
 It expresses the position $(x(t), y(t))$ as a function of time $t$.
 But suppose we want to determine whether the projectile hits a target 10 meters away.
 There are at least three ways of answering this:
-* 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.
-* 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.
-* Solve the ODE for $y(x)$ and evaluate it at 10 meters.
+
+    - 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.
+    - 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.
+    - Solve the ODE for $y(x)$ and evaluate it at 10 meters.
 
 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$.
+This transformation is well-defined for any non-zero horizontal velocity $v$, so $x$ and $t$ are one-to-one.
+
 ```@example changeivar
 M2 = change_independent_variable(M1, x)
 M2s = structural_simplify(M2; allow_symbolic = true)
@@ -46,17 +50,21 @@ M2s = structural_simplify(M2; allow_symbolic = true)
 @assert allequal([M2.x, M2s.x]) # hide
 @assert allequal([y, M1s.y]) # hide
 @assert allunique([M1.x, M1.y, M2.y, M2s.y]) # hide
+@assert length(equations(M2s)) == 2 # hide
 M2s # display this # hide
 ```
+
 The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`.
 
 !!! warn
+    
     At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables.
     Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables.
     It can be instructive to inspect these yourself to see their subtle differences.
 
 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.
 It is straightforward to evolve the ODE for 10 meters and plot the resulting trajectory $y(x)$:
+
 ```@example changeivar
 using OrdinaryDiffEq, Plots
 prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
@@ -65,67 +73,75 @@ plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
 ```
 
 !!! tip "Usage tips"
+    
     Look up the documentation of [`change_independent_variable`](@ref) for tips on how to use it.
-
+    
     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.
 
 ## 2. Alleviating stiffness by changing to logarithmic time
 
 In cosmology, the [Friedmann equations](https://en.wikipedia.org/wiki/Friedmann_equations) describe the expansion of the universe.
 In terms of conformal time $t$, they can be written
+
 ```@example changeivar
 @variables a(t) Ω(t)
 a = GlobalScope(a) # global var needed by all species
 function species(w; kw...)
-    eqs = [D(Ω) ~ -3(1 + w) * D(a)/a * Ω]
+    eqs = [D(Ω) ~ -3(1 + w) * D(a) / a * Ω]
     return ODESystem(eqs, t, [Ω], []; kw...)
 end
-@named r = species(1//3) # radiation
+@named r = species(1 // 3) # radiation
 @named m = species(0) # matter
 @named Λ = species(-1) # dark energy / cosmological constant
-eqs = [
-    Ω ~ r.Ω + m.Ω + Λ.Ω # total energy density
-    D(a) ~ √(Ω) * a^2
-]
-initialization_eqs = [
-    Λ.Ω + r.Ω + m.Ω ~ 1
-]
+eqs = [Ω ~ r.Ω + m.Ω + Λ.Ω, D(a) ~ √(Ω) * a^2]
+initialization_eqs = [Λ.Ω + r.Ω + m.Ω ~ 1]
 M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M)
 M1 = compose(M1, r, m, Λ)
 M1s = structural_simplify(M1)
 ```
+
 Of course, we can solve this ODE as it is:
+
 ```@example changeivar
 prob = ODEProblem(M1s, [M1s.a => 1.0, M1s.r.Ω => 5e-5, M1s.m.Ω => 0.3], (0.0, -3.3), [])
 sol = solve(prob, Tsit5(); reltol = 1e-5)
 @assert Symbol(sol.retcode) == :Unstable # surrounding text assumes this was unstable # hide
-plot(sol, idxs = [M1.a, M1.r.Ω/M1.Ω, M1.m.Ω/M1.Ω, M1.Λ.Ω/M1.Ω])
+plot(sol, idxs = [M1.a, M1.r.Ω / M1.Ω, M1.m.Ω / M1.Ω, M1.Λ.Ω / M1.Ω])
 ```
-The solver becomes unstable due to stiffness.
+
+But the solver becomes unstable due to stiffness.
 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.
-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.
+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.
 
 It is therefore common to write these ODEs in terms of $b = \ln a$.
-To do this, we will change the independent variable in two stages; from $t$ to $a$ to $b$.
-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.
+To do this, we will change the independent variable in two stages; first from $t$ to $a$, and then from $a$ to $b$.
+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.
 First, we transform from $t$ to $a$:
+
 ```@example changeivar
 M2 = change_independent_variable(M1, M1.a)
+@assert !ModelingToolkit.isautonomous(M2) # hide
+M2 # hide
 ```
+
 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$:
+
 ```@example changeivar
 a = M2.a
 Da = Differential(a)
 @variables b(a)
 M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)])
 ```
+
 We can now solve and plot the ODE in terms of $b$:
+
 ```@example changeivar
 M3s = structural_simplify(M3; allow_symbolic = true)
 prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), [])
 sol = solve(prob, Tsit5())
 @assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide
-plot(sol, idxs = [M3.r.Ω/M3.Ω, M3.m.Ω/M3.Ω, M3.Λ.Ω/M3.Ω])
+plot(sol, idxs = [M3.r.Ω / M3.Ω, M3.m.Ω / M3.Ω, M3.Λ.Ω / M3.Ω])
 ```
+
 Notice that the variables vary "more nicely" with respect to $b$ than $t$, making the solver happier and avoiding numerical problems.
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 51156ba993..41f9292c73 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -255,7 +255,8 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        tunable_parameters, isirreducible, getdescription, hasdescription,
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
-export ode_order_lowering, dae_order_lowering, liouville_transform, change_independent_variable
+export ode_order_lowering, dae_order_lowering, liouville_transform,
+       change_independent_variable
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
diff --git a/src/systems/diffeqs/basic_transformations.jl b/src/systems/diffeqs/basic_transformations.jl
index c05e5e4394..a08c83ffb6 100644
--- a/src/systems/diffeqs/basic_transformations.jl
+++ b/src/systems/diffeqs/basic_transformations.jl
@@ -47,11 +47,17 @@ function liouville_transform(sys::AbstractODESystem; kwargs...)
     neweq = D(trJ) ~ trJ * -tr(calculate_jacobian(sys))
     neweqs = [equations(sys); neweq]
     vars = [unknowns(sys); trJ]
-    ODESystem(neweqs, t, vars, parameters(sys); checks = false, name = nameof(sys), kwargs...)
+    ODESystem(
+        neweqs, t, vars, parameters(sys);
+        checks = false, name = nameof(sys), kwargs...
+    )
 end
 
 """
-    change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_old_diff = false, simplify = true, fold = false, kwargs...)
+    change_independent_variable(
+        sys::AbstractODESystem, iv, eqs = [];
+        add_old_diff = false, simplify = true, fold = false
+    )
 
 Transform the independent variable (e.g. ``t``) of the ODE system `sys` to a dependent variable `iv` (e.g. ``u(t)``).
 The transformation is well-defined when the mapping between the new and old independent variables are one-to-one.
@@ -68,13 +74,13 @@ Any extra equations `eqs` involving the new and old independent variables will b
 # Usage before structural simplification
 
 The variable change must take place before structural simplification.
-Subsequently, consider passing `allow_symbolic = true` to `structural_simplify(sys)` to reduce the number of unknowns, with the understanding that the transformation is well-defined.
+In following calls to `structural_simplify`, consider passing `allow_symbolic = true` to avoid undesired constraint equations between between dummy variables.
 
 # Usage with non-autonomous systems
 
-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`.
+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))``).
+Otherwise the transformed system can be underdetermined.
+If an algebraic relation is not known, consider using `add_old_diff` instead.
 
 # Usage with hierarchical systems
 
@@ -102,7 +108,10 @@ julia> unknowns(M)
  yˍx(x)
 ```
 """
-function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; 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
+)
     iv2_of_iv1 = unwrap(iv) # e.g. u(t)
     iv1 = get_iv(sys) # e.g. t
 
@@ -114,6 +123,7 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_o
         error("Variable $iv is not a function of the independent variable $iv1 of the system $(nameof(sys))!")
     end
 
+    # Set up intermediate and final variables for the transformation
     iv1name = nameof(iv1) # e.g. :t
     iv2name = nameof(operation(iv2_of_iv1)) # e.g. :u
     iv2, = @independent_variables $iv2name # e.g. u
@@ -127,23 +137,22 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_o
     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
+    @set! sys.eqs = [get_eqs(sys); eqs] # add extra equations we derived
+    @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)
 
-    # 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!
+    # 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)
     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)?
+            is_function_of_iv1 = iscall(var) && isequal(only(arguments(var)), iv1) # of the form f(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_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
-        # 2) Expand with chain rule until nothing changes anymore
+        # 2) Repeatedly expand chain rule until nothing changes anymore
         orgex = nothing
         while !isequal(ex, orgex)
             orgex = ex # save original
@@ -156,22 +165,25 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_o
         return ex
     end
 
+    # Use the utility function to transform everything in the system!
     function transform(sys::AbstractODESystem)
         eqs = map(transform, get_eqs(sys))
         unknowns = map(transform, get_unknowns(sys))
         unknowns = filter(var -> !isequal(var, iv2), unknowns) # remove e.g. u
         ps = map(transform, get_ps(sys))
-        ps = filter(!isinitial, ps) # remove Initial(...) # # TODO: shouldn't have to touch this
+        ps = filter(!isinitial, ps) # remove Initial(...) # TODO: shouldn't have to touch this
         observed = map(transform, get_observed(sys))
         initialization_eqs = map(transform, get_initialization_eqs(sys))
         parameter_dependencies = map(transform, get_parameter_dependencies(sys))
-        defaults = Dict(transform(var) => transform(val) for (var, val) in get_defaults(sys))
+        defaults = Dict(transform(var) => transform(val)
+        for (var, val) in get_defaults(sys))
         guesses = Dict(transform(var) => transform(val) for (var, val) in get_guesses(sys))
-        assertions = Dict(transform(condition) => msg for (condition, msg) in get_assertions(sys))
+        assertions = Dict(transform(ass) => msg for (ass, msg) in get_assertions(sys))
         systems = get_systems(sys) # save before reconstructing system
         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,
+            eqs, iv2, unknowns, ps; observed, initialization_eqs,
+            parameter_dependencies, defaults, guesses,
             assertions, name = nameof(sys), description = description(sys)
         )
         systems = map(transform, systems) # recurse through subsystems
@@ -182,6 +194,5 @@ function change_independent_variable(sys::AbstractODESystem, iv, eqs = []; add_o
         end
         return sys
     end
-
     return transform(sys)
 end
diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl
index dd6d677b60..544a89cd29 100644
--- a/test/basic_transformations.jl
+++ b/test/basic_transformations.jl
@@ -32,16 +32,19 @@ end
     M2 = change_independent_variable(M1, y)
     @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])
+    @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
     @variables x(t) y(t) z(t) s(t)
     eqs = [
-        D(x) ~ y
-        D(D(y)) ~ 2 * x * D(y)
-        z ~ x + D(y)
-        D(s) ~ 1 / (2*s)
+        D(x) ~ y,
+        D(D(y)) ~ 2 * x * D(y),
+        z ~ x + D(y),
+        D(s) ~ 1 / (2 * s)
     ]
     initialization_eqs = [x ~ 1.0, y ~ 1.0, D(y) ~ 0.0]
     M1 = ODESystem(eqs, t; initialization_eqs, name = :M)
@@ -55,8 +58,8 @@ end
     sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10)
     ts = range(0.0, 1.0, length = 50)
     ss = .√(ts)
-    @test all(isapprox.(sol1(ts, idxs=M1.x), sol2(ss, idxs=M2.x); atol = 1e-7)) &&
-          all(isapprox.(sol1(ts, idxs=M1.y), sol2(ss, idxs=M2.y); atol = 1e-7))
+    @test all(isapprox.(sol1(ts, idxs = M1.x), sol2(ss, idxs = M2.x); atol = 1e-7)) &&
+          all(isapprox.(sol1(ts, idxs = M1.y), sol2(ss, idxs = M2.y); atol = 1e-7))
 end
 
 @testset "Change independent variable (Friedmann equation)" begin
@@ -64,15 +67,15 @@ end
     D = Differential(t)
     @variables a(t) ȧ(t) Ω(t) ϕ(t)
     a, ȧ = GlobalScope.([a, ȧ])
-    species(w; kw...) = ODESystem([D(Ω) ~ -3(1 + w) * D(a)/a * Ω], t, [Ω], []; kw...)
-    @named r = species(1//3)
+    species(w; kw...) = ODESystem([D(Ω) ~ -3(1 + w) * D(a) / a * Ω], t, [Ω], []; kw...)
+    @named r = species(1 // 3)
     @named m = species(0)
     @named Λ = species(-1)
     eqs = [
-        Ω ~ r.Ω + m.Ω + Λ.Ω
-        D(a) ~ ȧ
-        ȧ ~ √(Ω) * a^2
-        D(D(ϕ)) ~ -3*D(a)/a*D(ϕ)
+        Ω ~ r.Ω + m.Ω + Λ.Ω,
+        D(a) ~ ȧ,
+        ȧ ~ √(Ω) * a^2,
+        D(D(ϕ)) ~ -3 * D(a) / a * D(ϕ)
     ]
     M1 = ODESystem(eqs, t, [Ω, a, ȧ, ϕ], []; name = :M)
     M1 = compose(M1, r, m, Λ)
@@ -80,20 +83,22 @@ end
     # Apply in two steps, where derivatives are defined at each step: first t -> a, then a -> b
     M2 = change_independent_variable(M1, M1.a)
     M2c = complete(M2) # just for the following equation comparison (without namespacing)
-    a, ȧ, Ω, Ωr, Ωm, ΩΛ, ϕ, aˍt = M2c.a, M2c.ȧ, M2c.Ω, M2c.r.Ω, M2c.m.Ω, M2c.Λ.Ω, M2c.ϕ, M2c.aˍt
+    a, ȧ, Ω, ϕ, aˍt = M2c.a, M2c.ȧ, M2c.Ω, M2c.ϕ, M2c.aˍt
+    Ωr, Ωm, ΩΛ = M2c.r.Ω, M2c.m.Ω, M2c.Λ.Ω
     Da = Differential(a)
     @test Set(equations(M2)) == Set([
-        aˍt ~ ȧ # dummy equation
-        Ω ~ Ωr + Ωm + ΩΛ
-        ȧ ~ √(Ω) * a^2
-        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
+        aˍt ~ ȧ, # dummy equation
+        Ω ~ Ωr + Ωm + ΩΛ,
+        ȧ ~ √(Ω) * a^2,
+        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
     ])
 
     @variables b(M2.a)
-    M3 = change_independent_variable(M2, b, [Differential(M2.a)(b) ~ exp(-b), M2.a ~ exp(b)])
+    extraeqs = [Differential(M2.a)(b) ~ exp(-b), M2.a ~ exp(b)]
+    M3 = change_independent_variable(M2, b, extraeqs)
 
     M1 = structural_simplify(M1)
     M2 = structural_simplify(M2; allow_symbolic = true)
@@ -103,55 +108,69 @@ 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)
+    Mt = ODESystem([D(x) ~ 2 * x, D(y1) ~ y1], t; name = :M)
     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ˍ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 with 1st order horizontal equation)" begin
     @variables x(t) y(t)
-    @parameters g v # gravitational acceleration and constant horizontal velocity
+    @parameters g=9.81 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, ...
     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], (0.0, 20.0), [g => 9.81, v => 10.0]) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities
+    u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0]
+    p = [v => 10.0]
+    prob = ODEProblem(Mx, u0, (0.0, 20.0), p) # 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)
+    @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
+    @parameters g = 9.81 # 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
+    u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0, Mx.xˍt => 10.0]
+    prob = ODEProblem(Mx, u0, (0.0, 20.0), []) # 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)
+    @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 3rd order nonlinear system)" begin
     @independent_variables t
     D = Differential(t)
     @variables x(t) y(t)
-    M1 = ODESystem([
-        (D^3)(y) ~ D(x)^2 + (D^2)(y^2) |> expand_derivatives
+    eqs = [
+        (D^3)(y) ~ D(x)^2 + (D^2)(y^2) |> expand_derivatives,
         D(x)^2 + D(y)^2 ~ x^4 + y^5 + t^6
-    ], t; name = :M)
+    ]
+    M1 = ODESystem(eqs, t; name = :M)
     M2 = change_independent_variable(M1, x; add_old_diff = true)
     @test_nowarn structural_simplify(M2)
 
     # Compare to pen-and-paper result
     @variables x xˍt(x) xˍt(x) y(x) t(x)
     Dx = Differential(x)
-    areequivalent(eq1, eq2) = isequal(expand(eq1.lhs - eq2.lhs), 0) && isequal(expand(eq1.rhs - eq2.rhs), 0)
-    @test areequivalent(equations(M2)[1], xˍt^3*(Dx^3)(y) + xˍt^2*Dx(y)*(Dx^2)(xˍt) + xˍt*Dx(y)*(Dx(xˍt))^2 + 3*xˍt^2*(Dx^2)(y)*Dx(xˍt) ~ xˍt^2 + 2*xˍt^2*Dx(y)^2 + 2*xˍt^2*y*(Dx^2)(y) + 2*y*Dx(y)*Dx(xˍt)*xˍt)
-    @test areequivalent(equations(M2)[2], xˍt^2 + xˍt^2*Dx(y)^2 ~ x^4 + y^5 + t^6)
-    @test areequivalent(equations(M2)[3], Dx(t) ~ 1 / xˍt)
+    areequivalent(eq1, eq2) = isequal(expand(eq1.lhs - eq2.lhs), 0) &&
+                              isequal(expand(eq1.rhs - eq2.rhs), 0)
+    eq1lhs = xˍt^3 * (Dx^3)(y) + xˍt^2 * Dx(y) * (Dx^2)(xˍt) +
+             xˍt * Dx(y) * (Dx(xˍt))^2 +
+             3 * xˍt^2 * (Dx^2)(y) * Dx(xˍt)
+    eq1rhs = xˍt^2 + 2 * xˍt^2 * Dx(y)^2 +
+             2 * xˍt^2 * y * (Dx^2)(y) +
+             2 * y * Dx(y) * Dx(xˍt) * xˍt
+    eq1 = eq1lhs ~ eq1rhs
+    eq2 = xˍt^2 + xˍt^2 * Dx(y)^2 ~ x^4 + y^5 + t^6
+    eq3 = Dx(t) ~ 1 / xˍt
+    @test areequivalent(equations(M2)[1], eq1)
+    @test areequivalent(equations(M2)[2], eq2)
+    @test areequivalent(equations(M2)[3], eq3)
 end
 
 @testset "Change independent variable (registered function / callable parameter)" begin
@@ -161,35 +180,38 @@ end
     @parameters f::LinearInterpolation (fc::LinearInterpolation)(..) # non-callable and callable
     callme(interp::LinearInterpolation, input) = interp(input)
     @register_symbolic callme(interp::LinearInterpolation, input)
-    M1 = ODESystem([
-        D(x) ~ 2*t
-        D(y) ~ 1*fc(t) + 2*fc(x) + 3*fc(y) + 1*callme(f, t) + 2*callme(f, x) + 3*callme(f, y)
-    ], t; name = :M)
+    eqs = [
+        D(x) ~ 2t,
+        D(y) ~ 1fc(t) + 2fc(x) + 3fc(y) + 1callme(f, t) + 2callme(f, x) + 3callme(f, y)
+    ]
+    M1 = ODESystem(eqs, t; name = :M)
 
     # Ensure that interpolations are called with the same variables
     M2 = change_independent_variable(M1, x, [t ~ √(x)])
     @variables x xˍt(x) y(x) t(x)
     Dx = Differential(x)
     @test Set(equations(M2)) == Set([
-        t ~ √(x)
-        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)
+        t ~ √(x),
+        xˍt ~ 2t,
+        xˍt * Dx(y) ~ 1fc(t) + 2fc(x) + 3fc(y) +
+                      1callme(f, t) + 2callme(f, x) + 3callme(f, y)
     ])
 
     _f = LinearInterpolation([1.0, 1.0], [-100.0, +100.0]) # constant value 1
     M2s = structural_simplify(M2; allow_symbolic = true)
     prob = ODEProblem(M2s, [M2s.y => 0.0], (1.0, 4.0), [fc => _f, f => _f])
     sol = solve(prob, Tsit5(); abstol = 1e-5)
-    @test isapprox(sol(4.0, idxs=M2.y), 12.0; atol = 1e-5) # Anal solution is D(y) ~ 12 => y(t) ~ 12*t + C => y(x) ~ 12*√(x) + C. With y(x=1)=0 => 12*(√(x)-1), so y(x=4) ~ 12
+    @test isapprox(sol(4.0, idxs = M2.y), 12.0; atol = 1e-5) # Anal solution is D(y) ~ 12 => y(t) ~ 12*t + C => y(x) ~ 12*√(x) + C. With y(x=1)=0 => 12*(√(x)-1), so y(x=4) ~ 12
 end
 
 @testset "Change independent variable (errors)" begin
     @variables x(t) y z(y) w(t) v(t)
     M = ODESystem([D(x) ~ 1, v ~ x], t; name = :M)
-    @test_throws "structurally simplified" change_independent_variable(structural_simplify(M), y)
-    @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)
+    Ms = structural_simplify(M)
+    @test_throws "structurally simplified" change_independent_variable(Ms, y)
+    @test_throws "not a function of" change_independent_variable(M, y)
+    @test_throws "not a function of" change_independent_variable(M, z)
     @variables x(..) # require explicit argument
-    M = ODESystem([D(x(t)) ~ x(t-1)], t; name = :M)
+    M = ODESystem([D(x(t)) ~ x(t - 1)], t; name = :M)
     @test_throws "DDE" change_independent_variable(M, x(t))
 end

From 7f821e21fecdfdf857dc220fd4856d5c6388a9b3 Mon Sep 17 00:00:00 2001
From: Herman Sletmoen <hermansletmoen@gmail.com>
Date: Tue, 18 Mar 2025 20:09:10 +0100
Subject: [PATCH 32/32] Use t_nounits, D_nounits

---
 docs/src/tutorials/change_independent_variable.md | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md
index e1243fb2ee..4b7ba037fe 100644
--- a/docs/src/tutorials/change_independent_variable.md
+++ b/docs/src/tutorials/change_independent_variable.md
@@ -18,8 +18,7 @@ Consider a projectile shot with some initial velocity in a vertical gravitationa
 
 ```@example changeivar
 using ModelingToolkit
-@independent_variables t
-D = Differential(t)
+using ModelingToolkit: t_nounits as t, D_nounits as D
 @variables x(t) y(t)
 @parameters g=9.81 v # gravitational acceleration and horizontal velocity
 eqs = [D(D(y)) ~ -g, D(x) ~ v]
@@ -128,7 +127,7 @@ Unlike the original, notice that this system is *non-autonomous* because the ind
 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$:
 
 ```@example changeivar
-a = M2.a
+a = M2.a # get independent variable of M2
 Da = Differential(a)
 @variables b(a)
 M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)])