Merge pull request #714 from SciML/myb/order_lowering

Fix order lowering

Author: YingboMa

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "4.4.2"
+version = "4.4.3"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

src/systems/diffeqs/first_order_transform.jl

@@ -15,15 +15,12 @@ function ode_order_lowering(eqs, iv, states)
     diff_eqs = Equation[]
     diff_vars = []
     alge_eqs = Equation[]
-    alge_vars = []
 
-    for (i, (eq, ss)) ∈ enumerate(zip(eqs, states))
-        if _iszero(eq.lhs)
-            push!(alge_vars, ss)
+    for (i, eq) ∈ enumerate(eqs)
+        if !isdiffeq(eq)
             push!(alge_eqs, eq)
         else
             var, maxorder = var_from_nested_derivative(eq.lhs)
-            # only save to the dict when we need to lower the order to save memory
             maxorder > get(var_order, var, 1) && (var_order[var] = maxorder)
             var′ = lower_varname(var, iv, maxorder - 1)
             rhs′ = diff2term(eq.rhs)
@@ -45,5 +42,5 @@ function ode_order_lowering(eqs, iv, states)
     end
 
     # we want to order the equations and variables to be `(diff, alge)`
-    return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, alge_vars))
+    return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, setdiff(states, diff_vars)))
 end

test/odesystem.jl

@@ -145,7 +145,7 @@ eqs = [D(x) ~ -A*x,
        D(y) ~ A*x - B*_x]
 de = ODESystem(eqs)
 @test begin
-    local f
+    local f, du
     f = eval(generate_function(de, [x,y], [A,B,C])[2])
     du = [0.0,0.0]
     f(du, [1.0,2.0], [1,2,3], 0.0)
@@ -225,3 +225,12 @@ using ModelingToolkit
 @derivatives D'~t
 sys = ODESystem([D(x) ~ a])
 @test sys.eqs[1].rhs isa Sym
+
+# issue 708
+@parameters t a
+@variables x(t) y(t) z(t)
+@derivatives D'~t
+sys = ODESystem([D(x) ~ y, 0 ~ x + z, 0 ~ x - y], t, [z, y, x], [])
+sys2 = ode_order_lowering(sys)
+M = ModelingToolkit.calculate_massmatrix(sys2)
+@test M == Diagonal([1, 0, 0])