Mass matrix for ODE with Julia (Jacobian sparsity)

Question

while implementing a stiff matrix ODE in Julia, I notice that the inverse Mass matrix required for the state-space representation (2nd order ODE), leads to a very dense Jacobian matrix: for dx.=[-invM*C -invM*K; eye(Float64, 406, 406) zeros(Float64, 406, 406)]*x

If instead we remove the inv(Ms) expression, the "Jacobian" of the right-hand side of the ODE is very sparse, which is no suprise: for "Mass"*dx.=[-C -K; eye(Float64, 406, 406) zeros(Float64, 406, 406)]*x

Is it possible to improve performance by providing the mass on the left-hand side of the equation, in the form:

[Ms zeros(Float64, 406, 406); zeros(Float64, 406,406) eye(Float, 406,406)]

I guess this option is available in DifferentialEquations.jl ?

Thank you

Mass matrices are described in the stiff ODE tutorial under the title Handling Mass Matrices. Is this what you were looking for? — Chris Rackauckas

Chris Rackauckas Chris Rackauckas · Accepted Answer · 2020-01-07T13:36:48

Mass matrices are described in the stiff ODE tutorial under the title Handling Mass Matrices. You define the mass matrix by adding it to the ODEFunction as a keyword argument. For example:

using DifferentialEquations
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  y₁ + y₂ + y₃ - 1
  nothing
end
M = [1. 0  0
     0  1. 0
     0  0  0]
f = ODEFunction(rober,mass_matrix=M)
prob_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob_mm,Rodas5(),reltol=1e-8,abstol=1e-8)

plot(sol, xscale=:log10, tspan=(1e-6, 1e5), layout=(3,1))

is the mass matrix DAE form of the ROBER equation.

Mass matrix for ODE with Julia (Jacobian sparsity)

1 Answers