Change independent variable of 2nd-order systems

Author: hersle

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
 

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