peel off Constant in mass matrix

Author: shashi

src/simplify.jl

@@ -43,6 +43,7 @@ function simplify_constants(expr)
 end
 
 to_mtk(x) = x
+to_mtk(x::Number) = Constant(x)
 to_mtk(v::SymbolicUtils.Variable{T}) where {T} = Variable{T}(nameof(v))
 to_mtk(v::SymbolicUtils.Variable{FnType{X,Y}}) where {X,Y} = Variable{Y}(nameof(v))
 function to_mtk(expr::SymbolicUtils.Term)

src/systems/diffeqs/abstractodesystem.jl

@@ -114,6 +114,8 @@ function calculate_massmatrix(sys::AbstractODESystem, simplify=true)
         end
     end
     M = simplify ? simplify_constants.(M) : M
+    # M should only contain concrete numbers
+    M = map(x->x isa Constant ? x.value : x, M)
     M == I ? I : M
 end
 

test/derivatives.jl

@@ -6,6 +6,8 @@ using Test
 @variables x y z
 @derivatives D'~t D2''~t Dx'~x
 
+test_equal(a, b) = @test isequal(simplify_constants(a), simplify_constants(b))
+
 @test @macroexpand(@derivatives D'~t D2''~t) == @macroexpand(@derivatives (D'~t), (D2''~t))
 
 @test isequal(expand_derivatives(D(t)), 1)
@@ -38,15 +40,15 @@ eqs = [0 ~ σ*(y-x),
        0 ~ x*y - β*z]
 sys = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(sys)
-@test isequal(jac[1,1], -1σ)
-@test isequal(jac[1,2], σ)
-@test isequal(jac[1,3], 0)
-@test isequal(jac[2,1], -1z + ρ) # FIXME
-@test isequal(jac[2,2], -1)
-@test isequal(jac[2,3], -1x)
-@test isequal(jac[3,1], y)
-@test isequal(jac[3,2], x)
-@test isequal(jac[3,3], -1*β)
+test_equal(jac[1,1], -1σ)
+test_equal(jac[1,2], σ)
+test_equal(jac[1,3], 0)
+test_equal(jac[2,1], -1z + ρ) # FIXME
+test_equal(jac[2,2], -1)
+test_equal(jac[2,3], -1x)
+test_equal(jac[3,1], y)
+test_equal(jac[3,2], x)
+test_equal(jac[3,3], -1β)
 
 # Variable dependence checking in differentiation
 @variables a(t) b(a)

test/mass_matrix.jl

@@ -10,9 +10,9 @@ eqs = [D(y[1]) ~ -k[1]*y[1] + k[3]*y[2]*y[3],
 
 sys = ODESystem(eqs,t,y,k)
 M = calculate_massmatrix(sys)
-M == [1 0 0
-      0 1 0
-      0 0 0]
+@test M == [1 0 0
+            0 1 0
+            0 0 0]
 
 f = ODEFunction(sys)
 prob_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))

test/simplify.jl

@@ -28,4 +28,4 @@ d2 = Differential(x)(d1)
 d3 = Differential(x)(d2)
 
 @test simplified_expr(expand_derivatives(d3)) == :(0)
-#@test simplify_constants(x^0) == 1
+@test simplified_expr(simplify_constants(x^0)) == :(1)