BDF 与隐式龙格库塔时间步进

计算科学 时间积分 龙格库塔 平流扩散
2021-12-11 21:40:29

为什么要选择高阶隐式龙格库塔(IMRK)而不是 BDF 时间步长,有什么理由吗?BDF 对我来说似乎更容易,因为q阶段 IMRK 需求q每个时间步的线性求解。BDF 和 IMRK 的稳定性似乎是一个有争议的问题。我找不到任何比较/对比隐式时间步进器的资源。

如果有帮助,我的最终目标是为平流扩散 PDE 选择一个高阶隐式时间步进器。

1个回答

是的,由于某种原因,没有太多资源。很长一段时间以来,标准 goto 是“只使用 BDF 方法”。出于几个历史原因,这个口头禅是一成不变的:一方面,Gear 的代码是第一个广泛使用的刚性求解器,另一方面,MATLAB 套件没有/不包含隐式 RK 方法。然而,这种启发式方法并不总是正确的,我会说从测试中它通常是错误的。让我详细解释一下。

BDF 方法具有较高的固定成本

BDF 方法中的自适应时间步长和自适应顺序的成本非常高。必须重新计算系数,或者需要将值内插到不同的时间。为了使当前的 BDF 代码在这里做得更好,已经做了很多工作(有两种众所周知的实现“形式”试图以不同的方式处理这个问题),但这实际上只是一个非常困难的软件工程问题。正因为如此,实际上大部分 BDF 代码都是 Gear 原始代码的后代:Gear 的、vode、lsoda、Sundial 的 CVODE,甚至 DAE 求解器 DASKR 和 DASSL 都是相同代码的后代。

这意味着,如果您的问题“太小”,那么高固定成本确实很重要,而隐式 RK 方法会做得更好。

高阶 BDF 方法对于复杂的特征值非常不稳定

BDF 方法允许您控制最大阶数并使其更加保守,原因是:高阶 BDF 方法甚至不能很好地处理“中等大小”的复杂特征值。所以在这些情况下,为了稳定,他们不得不放弃他们的订单。这就是为什么 6 阶 BDF 方法虽然在技术上稳定,但经常被忽略,因为即使是 5 阶也已经存在问题(并且 6 阶更不稳定)。只有二​​阶是 A 稳定的,所以它总是可以回退到那里,但是步进是错误主导的。

因此,当在非平凡问题上使用 BDF 代码时,您不会一直得到 5 阶。震荡导致其订单下降。

BDF 方法的起始成本很高

自动启动的 BDF 方法的启动成本很高。它们从一个小的 Euler 步长开始,然后是一个小的 BDF-2 步长,等等。对于更短的时间积分,这具有非常不可忽略的影响。如果你经常停下来,比如事件处理,这将严重影响代码的效率。这就是为什么您很少看到在繁重的事件处理或延迟方程情况中使用多步方法的原因之一(延迟方程在每个足够高的阶不连续处重新开始,它在每个iτ对于每个延迟时间τ)

BDF 方法是多步方法,具有最佳的缩放比例

但这就是 BDF 方法对他们的作用:它们的扩展性最好。这是因为它们的函数调用最少。这错误地导致了您应该“始终默认使用 BDF 方法”的想法(出于上述原因,出于稳定性和效率的原因,这并不总是正确的),但这确实意味着,如果函数调用f足够昂贵,那么 BDF 方法将是最好的。通常这意味着对于足够大和足够刚性的 PDE 离散化,BDF 方法将是最有效的。也就是说,这个“足够大”非常依赖于问题,我发现 PDE 离散化在使用radau多达 1,000 个空间变量的方法时效果更好,所以 YMMV。唯一知道的方法是测试。

有哪些可用的基准?

Hairer 的书和 DiffEqBenchmarks(解释如下)在易于获得的工作精度图方面可能是最好的。Hairer 的求解常微分方程 II在第 154 页和第 155 页上有一堆工作精度图。他选择的问题的结果与我上面所说的一致,原因是我上面所说的:如果问题不是,隐式 RK 会更有效“足够大”。另一个需要注意的有趣的事情是,高阶 Rosenbrock 方法在他的许多测试(如Rodas)中在误差较高的情况下是最有效的,而隐式 RKradau5在误差较低的情况下是最有效的。但他的测试问题大多不是大型 PDE 的离散化,因此请考虑上述几点。

你如何测试/基准测试?

我喜欢在 Julia 中使用DifferentialEquations.jl对此进行测试(免责声明:我是开发人员之一)。这是在朱莉娅。编程语言真的应该在这里得到一个注释。请记住,随着函数调用成本的增加,BDF 方法会更好。在 R/MATLAB/Python 中,用户函数是优化求解器必须实际使用的唯一 R/MATLAB/Python 代码:如果您使用 SciPy 或 Sundials 包装器,除了您传递的函数之外,它都是 C/Fortran 代码. 这意味着,在动态语言(不是 Julia)中,BDF 方法比通常做得更好,因为函数调用非常未优化(这可能是 Shampine 包含ode15s在 MATLAB 套件中,但没有隐式 RK 方法的原因) .

因此,对于更优化的情况,Julia 是一个很好的测试场,因为如果函数调用是类型稳定的,那么您的f与任何 C/Fortran 函数一样快速/高效。DifferentialEquations.jl 的好处在于,您只需切换一行代码即可在大量算法之间进行测试。对于ODEProblem问题,您可以使用以下方法在两者之间切换:

@time sol = solve(prob,CVODE_BDF())
@time sol = solve(prob,radau())

其中第一个使用日晷CVODE(BDF 方法),第二个使用 Hairer 的radau(隐式 RK)。

对于任何ODEProblem情况,您都可以使用基准测试工具来查看不同算法如何针对不同的自适应容差进行缩放。一些结果发布在 DiffEqBenchmarks.jl 上例如,在 ROBER 问题(3 个刚性 ODE 的系统)中,您可以看到日晷实际上是不稳定的,并且以足够高的容差发散(而其他方法收敛得很好),显示了上面关于稳定性问题的说明。在范德波尔问题上,你可以看到它更像是一个洗牌。我有很多我没有发布的东西,但在我完成一些更高阶的 Rosenbrock 方法(Rodas 是这些方法的 Fortran 版本)之前可能不会得到它。

(注意:那些僵硬的基准需要更新。一方面,文本需要更新,因为出于某种原因 ODE.jl 的方法不同......)

外推方法和 RKC

还有一些外推方法,例如seulex用于解决僵硬问题的方法。这些是“无限自适应顺序”,但这仅意味着当您正在寻找非常低的错误时它们是渐近良好的(即希望以低于1e-10或更低的速度解决棘手的问题,在这种情况下您可能可以使用显式方法) . 但是,在大多数情况下,它们效率不高,应避免使用。

Runge-Kutta Chebyschev 是一种显式方法,它也适用于您应该考虑的僵硬问题。我还没有把它包裹在DifferentialEquations.jl 中,所以我自己没有任何确凿的证据证明它是如何公平的。

其他需要考虑的方法:刚性 PDE 的专用方法

应该注意的是,当您可以轻松计算雅可比行列式时,高阶 Rosenbrock 方法对于中小型刚性问题的表现非常好,甚至好很多倍。但是,对于某些 PDE,我相信对流扩散问题属于这一类,Rosenbrock 可能会失去一些精度。此外,他们需要非常准确的雅可比矩阵来保持其准确性。在 Julia 中,这很容易,因为求解器带有符号和自微分,这对于机器 epsilon 来说是正确的。然而,像 MATLAB 的东西ode23s可能会出现问题,因为这些实现使用有限差分。对于 BDF 和隐式 RK 方法,Jacobian 中的错误会导致收敛速度较慢,而对于 Rosenbrock,由于这些不是隐式方程,而是其中包含 Jacobian 反演的 RK 方法,这些只是失去了准确性的顺序。

其他要研究的方法是指数 RK 方法、指数时间差分 (ETD)、隐式积分因子 (IIF) 和指数 Rosenbrock 方法。前三个利用了在许多 PDE 离散化中,

ut=Au+f(u)

在哪里A是一个线性算子,这种线性可以让你准确地求解Au操作员部分。缺点是,生成的方法必须使用eA即使是密集的A是稀疏的,所以在 Krylov 方法中有很多工作。尽管如此,这些都是避免隐式方程的 A 稳定方法,因此它们对于 PDE 离散化非常有效。

指数 Rosenbrock 方法不需要A, 而是由Ju+g(u)在哪里J是雅可比行列式,并且f=Ju+g,但这是相同的想法。

还有其他方法:最新创作

完全隐式的方法显然适用于刚性方程。如果 PDE 不够大,则无法有效地“在空间上并行化”以充分利用 HPC。相反,您可以创建完全隐式的并行时间离散化,因此适用于刚性方程,但并行以充分利用硬件。XBraid是使用这种技术的求解器,主要方法是 PFFAST 和 parareal 方法。DifferentialEquations.jl 正在开发一种类似的神经网络方法。

同样,当您没有足够大的空间离散化来有效地并行化,但有可用的并行计算资源时,这是最佳选择。

结论:谨慎考虑

在数值差分历史中很长一段时间内,人们认为外推法是最好的,因为它们在技术上可以实现任何阶。这意味着,对于一个足够小的Δt,它们将是最有效的(因此,如果您正在寻找一个非常精确的解决方案,是的!)。然而,几乎没有人真正在这些方法有效的情况下工作,这就是为什么人们通常不再使用外推法的原因。这是一个谨慎的词:对渐近考虑持保留态度

BDF 方法是渐近最好的,但在大多数情况下,您可能不在该状态下工作。但是如果空间离散化有足够多的点,BDF 方法可以在空间中有效地并行化(通过并行化线性求解),并且函数调用最少,因此会做得最好。但是如果你的 PDE 离散化不够大,可以好好看看隐式 RK、Rosenbrock、指数 RK 等方法。使用像 DifferentialEquations.jl 这样的软件套件可以轻松地在不同方法之间进行交换,这对您了解哪种方法最适合您的问题域非常有帮助,因为在许多情况下,它无法提前知道。

[如果您有任何要添加到测试套件中的示例问题,请随时提供帮助。我想保留一个非常全面的资源。我希望“尽快”在可运行的笔记本形式中拥有海尔的所有基准,并希望其他代表科学领域的问题。任何帮助表示赞赏!]