CODEWITHSUNDEEP
matrixjulia
LysandrosAn
LysandrosAn

Reputation: 41

Mass matrix for ODE with Julia (Jacobian sparsity)

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: enter image description here 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: enter image description here 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

Upvotes: 0

Views: 659

Answers (1)

Chris Rackauckas
Chris Rackauckas

Reputation: 19142

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.

enter image description here

Upvotes: 1

Related Questions