Mass matrix DAEs

Specific Domains Modelling & Simulations
1

Hi, I am trying to run some DAE simulations using DifferentialEquations.jl, and am having some trouble. I am doing a very simple test function:

function f!(t, u, du)
    du[1] = u[2]
    du[2] = u[2] - 1.
    return
end

u0 = [0. 1.]
tspan = (0.0, 1.0)

M = zeros(2,2)
M[1,1] = 1.

m_ode_prob = ODEProblem(f!, u0, tspan, mass_matrix=M)
sol = solve(m_ode_prob, Rosenbrock23())

This produces a Base.LinAlg.SingularException(2)

Stacktrace 
 [1] A_ldiv_B! at ./linalg/lu.jl:238 [inlined]
 [2] A_ldiv_B!(::Array{Float64,2}, ::Base.LinAlg.LU{Float64,Array{Float64,2}}, ::Array{Float64,2}) at ./linalg/factorization.jl:55
 [3] (::DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!})(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}, ::Bool) at /home/frank/.julia/v0.6/DiffEqBase/src/linear_nonlinear.jl:10
 [4] ode_determine_initdt(::Array{Float64,2}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::Int64, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/initdt.jl:15
 [5] #init#1831(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:339
 [6] init(::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:51
 [7] #solve#1830(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [8] solve(::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [9] include_string(::String, ::String) at ./loading.jl:515

If instead I use the diagonal matrix:

M = zeros(2,2)
M[1,1] = 1.
M[2,2] = 1.

I get a DimensionMismatch(“B has leading dimension 1, but needs 2”).

Stacktrace 
 [1] getrs!(::Char, ::Array{Float64,2}, ::Array{Int64,1}, ::Array{Float64,2}) at ./linalg/lapack.jl:926
 [2] A_ldiv_B! at ./linalg/lu.jl:238 [inlined]
 [3] A_ldiv_B!(::Array{Float64,2}, ::Base.LinAlg.LU{Float64,Array{Float64,2}}, ::Array{Float64,2}) at ./linalg/factorization.jl:55
 [4] (::DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!})(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}, ::Bool) at /home/frank/.julia/v0.6/DiffEqBase/src/linear_nonlinear.jl:10
 [5] ode_determine_initdt(::Array{Float64,2}, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::Int64, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/initdt.jl:15
 [6] #init#1831(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:339
 [7] init(::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:51
 [8] #solve#1830(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [9] solve(::DiffEqBase.ODEProblem{Array{Float64,2},Float64,true,#f_3!,Void,Array{Float64,2},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Rosenbrock23{0,true,DiffEqBase.LinSolveFactorize{Base.LinAlg.#lufact!},DataType}) at /home/frank/.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [10] include_string(::String, ::String) at ./loading.jl:515

Am I doing something wrong here, is the mass matrix somehow specified differently than I thought (and is there documentation/tutorials/examples somewhere that would explain this, that I overlooked?)? Or should I fill an issue on github?

2

Your initial condition, u0, should be a Vector rather than a Matrix, so changing u0 to u0 = [0., 1.] fixes your problem.

3

Thanks for the quick reply! That does indeed fix the second error, with a non-degenerate mass matrix, but not the first one unfortunately.

4

Nope, that was me being dumb. I thought I came up with a great scheme for handling initial timestep calculations like in the standard ODE case the other day, but allowing for mass matrices. Turns out that scheme requires that the mass matrix isn’t singular… I’ll get that fixed up.

In the meantime, you can work around it yourself by setting an initial timestep.

function f!(t, u, du)
    du[1] = u[2]
    du[2] = u[2] - 1.
    return
end

u0 = [0.,1.]
tspan = (0.0, 1.0)

M = zeros(2,2)
M[1,1] = 1.

m_ode_prob = ODEProblem(f!, u0, tspan, mass_matrix=M)
sol = solve(m_ode_prob, Rosenbrock23(), dt = 0.1)

for example works just fine. If you want to match what other popular software does, Hairer says in his book:

H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
     THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
     ADAPTS ITS STEP SIZE (IF H=0.D0, THE CODE PUTS H=1.D-6).

so

sol = solve(m_ode_prob, Rosenbrock23(), dt = 1e-6)

matches the default you’d get from RADAU/RODAS/etc. I’ll get this fixed up.

5

Fixed on master and will be tagged soon.

6

Awesome, thanks!

7

Someone told me they found this old post, so let me update it. This is all documented here: DAE Solvers · DifferentialEquations.jl

All of the things mentioned in the OP post have been fixed since 2019.