Change independent variable of simple 1st order system

Author: hersle

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

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