所以关于这个有很多话要说，我们实际上会发表一篇论文来尝试总结一下，但让我把它缩小到可以放入一个快速 StackOverflow 帖子的内容。我将很早地发表一个声明并不断重复它：您无法将方法的效率与软件的效率分开。软件实现的细节是这个领域真正重要的。

本质上，最先进的方法完全取决于您正在研究的问题。DiffEqBenchmarks.jl有很多最新软件的基准测试（包括DifferentialEquations.jl、SUNDIALS、Hairer 的FORTRAN 东西、Shampine 的FORTRAN 东西、LSODA 等）。它基于Hairer II中的 Hairer 基准，但 DiffEqBenchmarks 与该软件更同步，并更彻底地探索了此处提到的许多方法（同时仍包括 Hairer 的 FORTRAN 方法）。可通过此进行基准测试的算法的完整列表可以在DifferentialEquations.jl 文档中找到，它涵盖了目前大多数已知的算法（并且在今年夏天之后变得更加完整）。并非文档中的所有方法都显示在每个图中，但您可以下载基准并添加/减去内容以及更改选项。

尽管您想将“软件”与“方法”分离，但许多细节的实际软件实现与实际使用的方法一样重要。因此，您可以在许多关于刚性问题的基准测试中看到，三种不同的 BDF 实现（Fortran 中的 ddebdf、LSODA、CVODE_BDF）考虑到它们如何选择时间步长、重用雅可比矩阵等，在不同问题上具有非常不同的效率。那些更多“工程”方面（尤其是 CVODE 的雅可比重用结构）是使软件高效或低效的原因。正如您从基准测试中看到的那样，改变选项也很重要。例如，SUNDIALS ARKODE 在不更改某些选项的情况下无法解决我们向它提出的大多数严格基准（这在 SUNDIALS 示例中也有说明，所以它 s不仅仅是我们的错误）！因此，在谈论方法以及阅读有关效率的论文时，请记住这一点。

鉴于所有这些警告，当您拥有足够小的系统时，高阶 Rosenbrock 方法似乎占主导地位，因为它们可以花费大量时间步长并且不需要牛顿迭代的稳定性（这将真正的稳定性限制在隐式方法上）。虽然以前的文献提到雅可比精度问题会阻止 Rosenbrock 方法在高阶上很好地收敛，但现代自动微分技术（如ForwardDiff.jl中的那些）用数值微分的雅可比精度规避了许多传统问题，这为该领域注入了新的活力（这也是 Julia 被用于该领域的原因之一）。BLAS 的并行化似乎只会阻碍这种大小（Jacobian <50x50，请参阅，工程细节很重要！），所以我怀疑在查看并行 Rosenbrock 方法（它允许多个 lu-factorizations通过创建类似 DIRK 的 Rosenbrock 并行完成），没有任何真正的软件实现来测试这个想法。我们将很快推出软件以开始基准测试，但文献只创建了双核并行的方法Rosenbrock 方法（我知道），因此专门用于 4 或 5 阶现代 4-8 核 CPU 的方法将是一个有趣的话题。

我要在这里提出的一个警告是，当您要求越来越低的容差时（比如1e-8），您总是会看到 Hairer 的 Radau 实现变得最有效。它也不能很好地扩展到更大的 ODE 大小，因为它使用比其他方法大得多的雅可比行列式作为 FIRK 方法，但是当需要高精度时，它的高阶确实成为导致令人印象深刻的效率的一个因素。

随着系统规模增加到 >100 ODE，其他因素开始逐渐变得更加重要。本质上，对于一些，计算 Jacobian 矩阵和分解 W 矩阵$W = I - \gamma f'$$\gamma \in \mathbb{R}$成为任何半隐式或隐式方法中的主导因素（因此任何隐式 RK、隐式多步、Rosenbrock）。当这种情况发生时，重用您已经计算的雅可比矩阵或重用您已经创建的分解雅可比矩阵的能力变得非常重要。因此，由于标准 Rosenbrock 方法的准确度取决于其 Jacobian 行列式的准确度，因此 Rosenbrock 方法开始落后，因为它们每一步都需要计算和反转一个新的 Jacobian 行列式。隐式方法没有这个问题，因为雅可比仅用作求解隐式方程的线搜索，因此实际上不会影响准确性。“坏雅可比行列式”确实导致必须进行更多线性求解才能使牛顿法收敛，但是与雅可比计算或反演相比，具有预分解矩阵的线性求解是渐近微不足道的。反向替换是$O(n^2)$，雅可比计算是，而反演是。这一事实使人们通常只认为后者是一个问题，实际时间和分析证明这是一种过度简化，因为计算雅可比需要执行中的计算，即 ODE 导数函数，而反演是纯线性代数（矩阵着色和稀疏雅可比矩阵的解析形式降低了复杂性，但仍需要以某种形式因此，在这个阶段，您必须根据计算的成本对问题进行更多细分，因为即使您也可能会遇到非常慢的问题$O(n^2)$$O(n^3)$$f$$f$$f$$N=100$$f$计算，因此调用可能会在矩阵求逆中占主导地位（这通常是 MATLAB 或 Python ODE 程序中的情况）。$f$

所以让我们假设你有一个主导的半大型 ODE 系统，我们将半大型定义为一个足够大的大小以关心这些更渐近的因素，但又足够小以至于密集/稀疏分解仍然适合内存。我们发现在这些情况下效果很好的是 SDIRK 方法之类的方法，因为它们可以以比 BDF 高得多的阶数稳定地步进（L 稳定性和 B 稳定性），同时在一定程度上重用雅可比行列式。发生的情况是这最小化计算（尽管目前确实需要更多的反演，这是我们正在调查的事情，再次是工程“数学后”挑战），因为 BDF 通过进行大量雅可比重用和抨击非常小的步骤变得高效（哈，打赌你没看到）同样倒置$f$$f$$W$表示连续的多个步骤。我们需要将其形式化为 DiffEqBenchmarks.jl 更多，但对此的一个很好的来源是来自 PR 的图，该图正在查看具有昂贵 f evals 的僵硬 1000 ODE 模型。在这种情况下，KenCarp4似乎是更好的方法之一，尽管我们正在研究替代的 (E)SDIRK 画面，但还没有找到更好的方法（5 阶似乎没有一个也能做到这一点，肯尼迪和卡彭特在他们的论文中指出了这一事实在 ESDIRK 方法上也是如此）。在这里可能证明有用的另一种策略是 RosW 方法，它们是 Rosenbrock 方法，当雅可比行列不准确时不会失去顺序。今年夏天将在该领域对这些方法进行全面测试。但是，很多效率将取决于开发有效的雅可比重用策略。

当最终到达实际上由矩阵求逆支配的领域时，BDF 似乎做得很好，原因比人们想象的更奇怪。不是 BDF，而是具体的 VODE 和 CVODE 实现线。原因是因为与其他刚性求解器相比，它们采取了非常小的步骤，从而限制了它们的dt变化，这有助于保持$W$矩阵相同 ($\gamma$总是与dt) 成正比，并减少所需的反转次数。再加上他们特殊的固定前导系数 (FLC) 形式，他们可以尽可能多地重复使用已经倒置的矩阵。然而，这里的替代策略，例如宽松的牛顿迭代，可能会被证明是成功的。这完全取决于工程，取决于如何处理自适应性并使其平滑。

下一个阶段是当你得到足够大的 ODE 时，稀疏分解不再适合内存。此时，您几乎必须更改非线性求解器策略。有两种方法正在研究中。一是做安德森加速。另一种是 Newton-Krylov 方法。这两种方法都不能仅仅将某些因素分解以使同一步骤中的后续阶段变得微不足道，这意味着多阶段方法（Rosenbrock、ESDIRK 等）现在比单阶段方法具有更多的非线性求解时间，因为它们有在每个阶段从头开始执行整个非线性求解。也就是说，还有一些工程因素进入了那里。例如，当您为诸如 ESDIRK 方法之类的东西实现非线性求解器时，您可以使用前面的阶段对下一阶段的非线性迭代的起始值进行外推。在许多情况下，$c_i$画面不是单调的，这反过来意味着它不是外推而是（低阶）插值，这意味着某些阶段的非线性迭代会在一两次收敛。因此，如果不提及阶段预测器的准确性（这是用于此过程的术语），在这里说多阶段方法处于不利地位并不是那么简单，因此“方法之外”的某些东西和部分软件过程再次影响实际结果。但是，在这个领域，BDF 似乎确实大放异彩，因为它只有一个阶段。然而，它使用如此小的步长这一事实有点令人担忧（缺乏高于 2 阶的 L 稳定性，并且缺乏优化的前导截断误差系数给出了这个属性）并且确实导致了相当多的线性求解，所以我不会

这集中在基准测试最多的方法上，但或者 EPRIK 和指数 Rosenbrock 方法是一个完全可供探索的替代分支。文献中报道了一些好的结果，DifferentialEquations.jl 确实在 ExponentialUtilities.jl 中有完整的 Krylov-expmv + Krylov 自适应设置，并以5 阶 EPIRK 方法实现，所以如果你想加入我们的基准测试团队，我们很想看看这些结果如何。当然，由于工程如此重要，这些都是相当新的，可能需要几轮分析。一些初步结果令人兴奋（有趣的是，虽然 ARKODE 在很多 ODE 基准测试中都失败了，但在离散化 PDE 基准测试中却相当出色）。

这些基准还强调的是，在 2019 年改进刚性求解器的方法之一是不将其作为单个求解f，而是作为一个拆分f = f1 + f2并通过 IMEX 积分器仅隐式处理一个部分。IMEX BDF（称为 SBDF）、IMEX ESDIRK、IMEX Rosenbrock 等都存在并且是需要更多测试的有趣路线。

此外，如果您只有“半刚性”，那么稳定显式方法是具有足够自适应稳定性的显式方法来处理一些刚性问题。实现存在但需要基准。我们目前正在进行基准大修以更新使用 Wea​​ve.jl 文件（用于自动运行）并将所有基准更新到 Julia v1.0，一旦完成，稳定的 RK、RosW、指数积分器和其他人将加入基准测试。当 GPU 和 TPU 集成到软件中时，我们的基准测试也缺少效率：原生 Julia 方法接近能够使用 GPU，因此测试我们使用它时会发生什么会很有趣。这给记忆带来了更高的压力，但对线性代数例程的成本压力却更低，这可能会使事情更趋向于$f$- 主导域。

因为这是一个巨大的话题，所以我也有可能错过了您特别想要更多讨论的无限选项之一。糟糕，我很抱歉，但这需要$O(n^5)$是时候审查所有方法了。5 年后再问，我们会有更多的基准 :)。

所以总的来说，你想要一个关于方法而不是软件的简短摘要，但是既然你问到什么是最先进的，我不得不用所有不同的方式来回应，结果实际上与软件相关联. 然而，这些是一堆不同的方向，似乎在不同的领域取得了丰硕的成果，JuliaDiffEq 团队正在继续实施、优化和基准测试。所以现在任何时候这些结果中的一些可能会发生变化，这就是它的样子，因为渐近因子和纯数学约束只在无穷大时才重要。如果您想了解有关该领域的更多信息，这样的摘要可能会有所帮助，但在您动手并开始进行基准测试和分析之前，您永远不会得到一个好主意。我们'

PS 要查找有关方法的论文，请参阅我们的引用页面。这不是最新的，因此在某些情况下您可能需要深入研究我们的问题。当然，我们总是欢迎PR 来改进我们的引文页面。

编辑：Yingbo Ma 分享了这个漂亮的小脚本，突出了 BDF 与其他方法的步进行为差异：

using OrdinaryDiffEq, Plots, ParameterizedFunctions, Sundials, ODEInterfaceDiffEq

hires = @ode_def Hires begin
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(hires,u0,(0.0,321.8122))
function plotdts(sols::Tuple, title; kwargs...)
    ts = map(s->s.second.t, sols)
    dts = map(t->[0; diff(t)], ts)
    xlim = max(map(maximum, ts)...)
    plt = plot(title=title, xlabel="t", ylabel="dt", ylims=(0, max(map(maximum, dts)...)*1.3), xlims=(0, xlim*1.5); kwargs...)
    vline!(plt, [xlim], lab="t1")
    for i in 1:length(ts)
        plot!(plt, ts[i], dts[i], lab=sols[i].first, marker=(2, ))
    end
    plt
end

plotdts(("CVODE_BDF"=>solve(prob, CVODE_BDF()),
         "KenCarp4"=>solve(prob, KenCarp4()),
         "radau5"=>solve(prob, radau5())), "Time vs Step Size", dpi=200)
savefig("$(homedir())/tvsdt.png")