odeint来自 SciPy 库的默认值为此处描述lsoda的积分器。然而，任何对渐近计算时间的简单描述都是不可能的。原因很多。

首先，让我描述一下算法。非刚性方程的常见多步算法是Adams-Moulton 方法。虽然这些是隐含的，但Adams-Bashforth 方法通常太不稳定而无法使用。为了解决隐式方程，一个简单的 Picard 迭代（或者如果你想花哨并让它更好一点，Anderson 加速，但 LSODA 使用 Picard），因为它很便宜。在这里使用 Picard 迭代会降低稳定性区域（如Hairer II中所讨论的），但这被认为是可以的，因为该算法适用于非刚性问题。这是亚当斯积分器。正如我在一篇博文中描述的那样，比较了多年来开发的不同 ODE 套件, 几乎所有的多步方法（除了一些来自 Shampine 的方法）都基于从 GEAR 开始的相同实现，因此lsode, lsoda, vode, CVode(SUNDIALS) 对于非刚性问题都有几乎相同的 Adams 实现。

对于刚性问题，首选多步方法是反向微分公式 (BDF) 方法。Wiki 页面非常好，显示稳定性如何变差 >5，因此通常没有人使用 6 阶。每一步都必须求解一个隐式方程，并且这些使用准牛顿方法来这样做，因为这不会像定点迭代那样增加理论稳定性。然而，为了降低成本，它是准牛顿而不是牛顿，因为在许多情况下，雅可比矩阵在两个步骤之间可能是相似的，因此重新使用相同的雅可比矩阵将允许它仍然收敛到正确的隐式解决方案而不必分解$(I - \Delta t J)$再次。

多步方法需要以前的时间点来计算它们的解决方案（查看公式）。因此，他们需要一些其他解决方案才能开始。因此，一个非常简单的方法是创建一个自适应排序方法。这些从一阶方法开始（亚当斯的欧拉，BDF 的反向欧拉），然后在计算一个值之后可以移动到二阶（基于某些标准），依此类推。随着阶数的增加，每步的误差减小，因此步长可以改变。因此，标准非刚性多步求解器是具有自适应步长的自适应 1-12 阶 Adams 方法，标准刚性多步求解器是具有自适应步长的自适应 1-5 BDF 方法。我什至不会在这里描述下一个顺序或时间步长的选择是如何发生的，因为这会很快变得相当复杂。

现在在 90 年代，像 Shampine 和 Petzold 这样的少数人认为可以将这些方法结合起来。Petzold 使用 Shampine 的想法创建了 LSODA，它与 LSODE 类似，只是它可以自适应地从 Adams 方法更改为 BDF 方法，反之亦然。基本上，它会进行特征值估计来计算方程的刚度，如果它太僵硬/太不僵硬，它会切换行为。

所以 LSODA 既是 Adams 方法又是 BDF 方法，它根据内部误差估计改变其步长和顺序，并且在许多情况下它必须为 BDF 分解雅可比行列式，但如果它意识到它不会在每一步都分解雅可比行列式：隐式方程将足够快地求解。完全不可能将任何渐近算法复杂性归因于此！

一个切线虽然：

我也被告知，使用这种方法时，无法估计错误

不，ODE 求解器确实有估计（局部）误差的内部方法，这实际上是自适应时间步进方法（如 LSODA）中使用的方法。其中大多数使用嵌入式方法，即单独的低阶方法，它是相同阶段的不同线性组合。由于是相同的阶段，它几乎可以免费计算。但是，没有真正使用 RK4 的嵌入式方法。RK4 有一个非常好的错误估计器，但在实践中非常罕见。事实上，它只用于依赖于状态的延迟求解器，Shampine 在他的论文中描述了它ddesd。DifferentialEquations.jl 实现已经实现了这个误差估计器（主要用于延迟微分方程），所以如果你好奇这里有一些开源代码，你可以看看了解它是如何完成的（以及更有效的可变数组版本）。不过，那是一些小众的东西。

那么好吧，那么您如何在实践中找到像 ODE 求解器的“算法复杂性”之类的东西呢？最好和最常用的方法是工作精度图。你要做的是你拿一些测试方程，以非常高的精度求解它们以获得参考解，然后用一堆不同的自适应容差和步长来求解它，以找到时间和误差之间的关系。这很关键：如果没有与错误的关系，时间就没有任何意义。

有三个大的基准集。第一个是Hairer I和Hairer II，它们都是求解算法的很好参考。第二个是这个 IVP 测试套件。不过这些都已经过时了，而且只有非常古老的算法。他们也不倾向于使用以多线程方式求解线性系统的方法，这在当今非常普遍并且对于刚性求解器效率非常必要。他们也没有LSODA。

第三个基准集是DiffEqBenchmarks.jl。它利用了DifferentialEquations.jl软件，该软件直接链接到Julia 编程语言中的大量方法，然后有一堆专门用于基准求解器的方法。这些基准包含lsoda大量其他选项。在非刚性问题上，lsoda 往往做得很好，但不如高阶 Runge-Kutta 方法（我将链接许多突出显示总体摘要的测试方程之一，请随意查看更多）。对于刚性问题，它往往在高容差方面做得很好，但在需要高精度时非常好（注意：Rodas5和lsoda很接近，但您可以通过与下一组进行比较来查看哪个是Rodas5...需要修复）。在某些情况下，在高精度下，它的时间可能会像这里看到的那样爆炸，这可能是因为自适应算法开始错误地切换到显式模式，时间步长足够小，只需要切换回来。

一个比较点是有序的。这是SciPy 使用的 Fortran 代码lsoda的 C 翻译的包装器，因此它具有相似的特征。SciPy 中的使用会非常慢，因为 Python 会像 Julia 那样链接编译后的导数函数，但除此之外，算法的行为不会有任何不同。f

概括

最后，让我总结一下该算法的一些特征。这是非常多范式的。它在几乎所有标准基准测试案例中都表现良好，尽管从来都不是最好的。但是它可以切换模式以处理非刚性问题，并且可以很好地处理刚性问题（但当容差变低时，它会出现刚性问题）。但是，如果您想要的只是一种在任何情况下都可以使用的算法，那么这将是一个不错的选择（这很可能是 SciPy 的默认设置）。

请注意，我说的是“标准基准案例”。在许多情况下，您不想使用此算法。例如，我在一篇博客文章中描述，如果您有很多事件，使用多步方法的成本会大大增加（因为它们不能再使用过去的时间点，需要“回到欧拉”）。此外，BDF 方法仅在二阶之前是 L 稳定的，因此它们在处理非常僵硬的 ODE 时存在问题（这就是为什么例如 MATLAB 建议您使用除ode15s如果问题太僵硬）。一般来说，多步方法对大型系统比对小型系统更有利，因为它们每个时间步只有一步，尽管其他方法（如龙格-库塔方法）可能会采取更大的步数，并且由于优化常数中的常数而得到相同的误差扩展（在此处描述）。

这可能是全面描述 LSODA 的算法复杂性（或求解 ODE 的效率）的最简单方法，LSODA 是odeintSciPy 1.0 中的默认选择。

编辑

为了将来参考，这是关于lsoda. lsoda是对lsodeAdams 和 BDF 方法的修改，所以这里没有提到它的名字，而是被称为“LSODE 包的修改版本”（这使得它很难找到搜索哈哈）。