针对刚性方程的 ODE 求解器的建议

计算科学 C++ 微分方程 C 朱莉娅
2021-12-26 07:36:44

我正在继续前博士的研究。我小组中的学生需要求解 ODE 系统。在技​​术说明中,他们写道:

玻尔兹曼方程组的行为在数值上是刚性的。这意味着建议对其数值解使用隐式方法以实现可接受的步长(以及因此可接受的执行时间和数值误差)。这里使用具有牛顿迭代的后向微分公式的 CVODE 作为 ODE 求解器。在每个外部步骤中分析计算完整的雅可比行列式。

我无权访问他们的代码,而且我是 C/C++ 的新手,所以我想知道是否

  1. 你能推荐一个更简单的 ODE 求解器吗?
  2. 是否有教程/书籍/示例来学习如何使用 CVODE?我已经查看了官方文档,但恐怕它对我来说太高级了。

如果需要,可以提供更多详细信息。谢谢!

2个回答

这是一个巨大的开放式问题,但我将从当前版本的DifferentialEquations.jl中复制推荐部分:

僵硬的问题

对于高公差 (>1e-2?) 的刚性问题,建议您使用 Rosenbrock23 或 TRBDF2。这些对振荡和巨大的刚度具有鲁棒性,但仅在需要低精度时才有效。Rosenbrock23 对于重新评估和重新分解雅可比行列式成本不太高的小型系统更有效,而对于足够大的系统,TRBDF2 将更有效。

在中等公差 (>1e-8?) 时,建议您使用 Rodas5、Rodas4(前者效率更高,但后者更可靠)、Kvaerno5 或 KenCarp4。作为原生的 DifferentialEquations.jl 求解器,许多 Julia 数值类型(例如 BigFloats、ArbFloats 或 DecFP)都可以使用。

为了在低容差 (<1e-9) 下更快地求解,但在使用 Vector{Float64} 时,请使用 radau。

对于渐近大型 ODE 系统(N>1000?),其中 f 非常昂贵且复特征值最小(低振荡),在这种情况下 QNDF 将是最有效的,但需要 Vector{Float64}。如果解是平滑的,QNDF 的表现也会出奇的好。但是,这种方法可以处理比其他方法更小的刚度,并且它的牛顿迭代可能会在精度较低的情况下失败。在此方案中要考虑的其他选择是 CVODE_BDF 和 lsoda。

刚性积分器的特殊性质

ImplicitMidpoint 是一个对称的辛积分器。梯形是具有自适应时间步长的对称(几乎是辛)积分器。ImplicitEuler 是通用算法的扩展,具有自适应时间步长和有效的准牛顿雅可比重用,它是双曲 PDE 的完全强稳定性保持 (SSP)。

请注意,Rodas4 在非线性抛物 PDE 的离散化上失去了准确性,因此建议您在那些是 3 阶的情况下将其替换为 Rodas4P。ROS3P 只是三阶并且在此类问题上实现了三阶,因此在这种情况下可以更有效。

(特别是这个页面

这些建议都得到SciMLBenchmarks 的支持, SciMLBenchmarks是一个不断更新的跨语言基准测试环境,用于科学机器学习活动,包括微分方程求解和参数估计。

此外,还有其他属性需要专门研究。u'=A(t)u例如,如果系统是这样的形式,可以使用 Magnus 方法,这将是便宜且稳定的。有关特殊形式的非自治线性 ODE 的方法列表,请参阅此页面。同样,在某些情况下,只有部分方程是刚性的,或者部分方程是刚性和线性的,在这种情况下,像 IMEX 或指数积分器这样的拆分方法是一个好主意。有关 IMEX 或指数积分器的说明,请参阅此页面

简而言之,总结是:

  1. Rosenbrock 方法对于较小的系统更有效,但需要更准确的雅可比行列式。如果您使用自动微分或对雅可比矩阵使用符号微分(如果它是 W 方法,但不多),并且您的系统足够小,您真的只想使用它们。如果这些情况成立,那么它们是迄今为止最快的。

  2. L-stable 方法通常是最有效的,高阶 SDIRK 和自适应时间自适应阶 BDF 方法是不错的选择。对于 BDF 部分,纯 Julia QNDF 在许多情况下(与一些细节和线性求解器有关)往往优于 CVODE_BDF,尝试 TRBDF2 或 KenCarp 方法通常是个好主意。

  3. 当您想要更低的容差时,请使用更高阶的方法。刚性 ODE 的最高阶方法是自适应阶 Radau 方法。


使用刚性 ODE 求解器

这与推荐几乎是不同的问题,但是DifferentialEquaitons.jl文档中有很多教程可以通过复制和粘贴来入门。有一个专门针对刚性 ODE 的完整教程最后,这将使用 CVODE_BDF 求解 Robertson 方程:

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] =  k₂*y₂^2
  nothing
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
sol = solve(prob,CVODE_BDF()) # Don't pass CVODE_BDF() and it will auto-select a method
plot(sol,tspan=(1e-2,1e5),xscale=:log10)

这可能是启动和运行 CVODE 的最简单方法,因为安装过程会自动构建所有相关库并包括 LAPACK/BLAS 的构建。

CVODE 是 SUNDIALS 软件包的一部分,它是开源的。它附带有数百页和许多示例代码的手册。很不错的软件,值得你用!