matrix - 具有 Julia 的 ODE 的质量矩阵（雅可比稀疏性）

Question

在 Julia 中实现刚性矩阵 ODE 时，我注意到状态空间表示（二阶 ODE）所需的反质量矩阵会导致非常密集的雅可比矩阵 ：dx.=[-invM*C -invM*K; eye(Float64, 406, 406) zeros(Float64, 406, 406)]*x

相反，如果我们删除 inv(Ms) 表达式，则 ODE 右侧的“Jacobian”非常稀疏，这不足为奇： 对于"Mass"*dx.=[-C -K; eye(Float64, 406, 406) zeros(Float64, 406, 406)]*x

是否可以通过在等式左侧提供质量来提高性能，形式为：

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

我猜这个选项在 DifferentialEquations.jl 中可用？

谢谢

标签： matrixjulia

Answer 1

质量矩阵在标题为处理质量矩阵的僵硬 ODE 教程中进行了描述。您可以通过将质量矩阵添加到ODEFunction作为关键字参数来定义质量矩阵。例如：

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))

是 ROBER 方程的质量矩阵 DAE 形式。