How to define ode without mass_matrix

General Usage
1

Dear All,
I am not so familary with Julia ode.
Here as to the ode definition using ODEfunction:
also online:CMSB-2023/UKF.jl at main · H-Sax/CMSB-2023 · GitHub
the code as below:
function NIK!(du, u, p, t)
pLV, psa, psv, Vlv, Qav, Qmv, Qs = u
τₑₛ, τₑₚ, Rmv, Zao, Rs, Csa, Csv, Eₘₐₓ, Eₘᵢₙ = p

# 1) Left Ventricle
du[1] = (Qmv - Qav) * ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) + pLV / ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) * DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
# 2) Systemic arteries 
du[2] = (Qav - Qs ) / Csa     
# 3) Venous
du[3] = (Qs - Qmv) / Csv 
# 4) Left Ventricular Volume
du[4] = Qmv - Qav 
# 5) Aortic Valve flow
du[5]    = Valve(Zao, (pLV - psa)) - Qav
# 6) Mitral Valve flow
du[6]   = Valve(Rmv, (psv - pLV)) - Qmv 
# 7) Systemic flow
du[7]     = (du[2] - du[3]) / Rs
nothing

end

M = [1. 0 0 0 0 0 0
0 1. 0 0 0 0 0
0 0 1. 0 0 0 0
0 0 0 1. 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 1. ]

Nik_ODE = ODEFunction(NIK!,mass_matrix=M)

u0 = [8.0, 8.0, 8.0, 265.0, 0.0, 0.0, 0.0]

p = [0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]

tspan = (0, τ * 30)

prob = ODEProblem(Nik_ODE, u0, tspan, p)

how to change the ode equation definition without mass matrix :Nik_ODE = ODEFunction(NIK!,mass_matrix=M)

thanks.

2

why don’t you want to use a mass matrix? It does change what the equations mean (specifically this is defining a Differential Algebraic equation in the semi-explicit form).

3

thanks.
as my next need is to estimate the parameter of p using Optimization.OptimizationProblem with the cost_function = build_loss_objective(prob, Tsit5(), L2Loss(t, data),Optimization.AutoForwardDiff(),maxiters = 10000, verbose = false,save_idxs = 1)

but the solver is not able to use mass matrices.
so…

the example is from:Getting Started with Optimization-Based ODE Parameter Estimation · DiffEqParamEstim.jl

4

yeah you just need to use a solver that can use mass matrices. I would probably suggest Rodas5P, but you might want to try FBDF as well. While you’re at it, you might want to rewrite

using LinearAlgebra
M = Diagonal([1., 1., 1., 0., 0. 1.])

which is a little easier to read and write.

5

thanks very much, I used Rodas5P according to your suggestion.