# Mass matrix DAEs

**URL:** <https://discourse.julialang.org/t/mass-matrix-daes/6805>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [October 31, 2017, 2:55pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805 "2017-10-31T14:55:32Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![FHell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fhell/32/2421_2.png) [@FHell](https://discourse.julialang.org/u/FHell)\
**Post date:** [October 31, 2017, 2:55pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/1 "2017-10-31T14:55:32Z")

</div>

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:

```julia
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**
>
> ```julia
> [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:

```julia
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**
>
> ```julia
> [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?

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [October 31, 2017, 4:01pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/2 "2017-10-31T16:01:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![FHell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fhell/32/2421_2.png) [@FHell](https://discourse.julialang.org/u/FHell)\
**Post date:** [October 31, 2017, 4:46pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/3 "2017-10-31T16:46:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 31, 2017, 5:05pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/4 "2017-10-31T17:05:11Z")

</div>

> [@FHell](#):
>
> but not the first one unfortunately.

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.

```julia
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:

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 31, 2017, 10:48pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/5 "2017-10-31T22:48:24Z")

</div>

Fixed on master and will be tagged soon.

---

<div class="post-metadata">

**Author:** ![FHell](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fhell/32/2421_2.png) [@FHell](https://discourse.julialang.org/u/FHell)\
**Post date:** [November 1, 2017, 2:56pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/6 "2017-11-01T14:56:28Z")

</div>

Awesome, thanks!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 23, 2022, 8:59pm UTC](https://discourse.julialang.org/t/mass-matrix-daes/6805/7 "2022-12-23T20:59:16Z")

</div>

Someone told me they found this old post, so let me update it. This is all documented here: [DAE Solvers · DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryDiffEq.jl-(Mass-Matrix))

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