关于数值积分器，“辛”是什么意思，SciPy 的 odeint 是否使用它们？

计算科学 颂 时间积分 一体化 微分方程
2021-12-10 19:53:45

在这篇评论中，我写道：

...默认 SciPy 积分器，我假设它只使用辛方法。

其中我指的是 SciPy's odeint，它使用“非刚性（Adams）方法”或“刚性（BDF）方法”。根据消息来源：

def odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0,
           ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0,
           hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
           mxords=5, printmessg=0):
    """
    Integrate a system of ordinary differential equations.

    Solve a system of ordinary differential equations using lsoda from the
    FORTRAN library odepack.

    Solves the initial value problem for stiff or non-stiff systems
    of first order ode-s::
        dy/dt = func(y, t0, ...)
    where y can be a vector.
    """

这是一个示例，我将卫星绕地球的轨道传播了三个月，只是为了表明它按预期进动。

我相信非辛积分器具有不希望的特性，即它们倾向于不保存能量（或其他量），因此例如在轨道力学中是不希望的。但是我不确定是什么使辛积分器辛。

是否有可能以一种直截了当且（相当）容易理解但并非不准确的方式来解释该属性是什么（使辛积分器成为辛积分器）？我是从积分器内部功能的角度来询问的，而不是它在测试中的表现。

我的怀疑是否正确，odeint只使用辛积分器？

2个回答

让我从更正开始。不，odeint没有任何辛积分器。不，辛积分并不意味着能量守恒。

辛是什么意思，什么时候应该使用它？

首先，辛是什么意思？辛意味着解存在于辛流形上。辛流形是由 2-形式定义的解集。辛流形的细节可能听起来像数学废话，所以它的要点是在这样的流形上的两组变量之间存在直接关系。这对物理学很重要的原因是因为哈密顿方程自然具有解驻留在相空间中的辛流形上，自然分裂是位置和动量分量。对于真正的哈密顿解，该相空间路径是恒定能量。

辛积分器是其解存在于辛流形上的积分器。由于离散化误差，当它求解哈密顿系统时，它并不能准确地得到流形上的正确轨迹。相反，该轨迹本身受到干扰O(Δtn)为了订单n从真实轨迹。然后由于该轨迹随时间的数值误差而存在线性漂移。正常的积分器往往具有二次（或更多）漂移，并且对此相空间路径（仅局部）没有任何良好的全局保证。

这往往意味着辛积分器倾向于比普通积分器更好地捕获长时间模式，因为没有漂移并且几乎可以保证周期性。这个笔记本很好地显示了开普勒问题的这些属性。第一张图片显示了我所说的解决方案的周期性。

在此处输入图像描述

使用来自 Kahan 的 6 阶辛积分器和来自 DifferentialEquations.jl 的 Li解决了这个问题。你可以看到能量不是完全守恒的，但它的变化取决于扰动解流形与真实流形的距离。但是由于数值解本身存在于一个辛流形上，它往往几乎是完全周期性的（你可以看到一些线性数值漂移），这使得它非常适合长期积分。如果你对 RK4 做同样的事情，你可能会遇到灾难：

在此处输入图像描述

您可以看到问题在于数值解中没有真正的周期性，因此随着时间的推移它往往会漂移。

这突出了选择辛积分器的真正原因：辛积分器擅长长期积分具有辛性质的问题（哈密顿系统）。所以让我们来看看一些事情。请注意，即使在辛问题上，您也并不总是需要辛积分器。对于这种情况，自适应 5 阶 Runge-Kutta 方法可以做得很好。这是Tsit5：

在此处输入图像描述

注意两件事。一，它获得了足够好的精度，以至于您无法在相空间图中看到实际的漂移。但是，在右侧，您可以看到存在这种能量漂移，因此，如果您进行足够长的积分，则此方法的效果不如具有周期性属性的求解方法。但这提出了一个问题，它在效率方面与仅极其准确地集成有何不同？嗯，这有点不太确定。在SciMLBenchmarks.jl 中，您可以找到一些调查此问题的基准。例如，这个笔记本查看来自四重玻色子模型的哈密顿方程系统上的能量误差与运行时间，并表明如果您想要非常高的精度，那么即使对于相当长的积分时间，仅使用高阶 RK 或龙格-库塔尼斯特罗姆（ RKN) 方法。这是有道理的，因为为了满足辛属性，积分器放弃了一些效率，并且几乎必须是固定的时间步长（有一些研究在后者方面取得了进展，但距离并不远）。

此外，从这两个笔记本中注意到，您也可以只采用标准方法，并在每个步骤（或每隔几个步骤）将其投影回解决方案流形。这就是使用 DifferentialEquations.jl ManifoldProjection 回调的示例正在做的事情。您会看到保证守恒定律得到维护，但每一步都需要额外的解决隐式系统的成本。您还可以使用完全隐式 ODE 求解器或奇异质量矩阵来添加守恒方程，但最终结果是这些方法的计算成本更高。

总而言之，您想要解决辛积分器的问题类别是那些在辛流形（哈密顿系统）上有解决方案的问题，您不想投入计算资源以获得非常精确的（容差<1e-12）解决方案，不需要精确的能量/等。保护。这突出表明，这一切都与长期集成特性有关，因此您不应该像某些文献所暗示的那样，不顾一切地涌向它们。但在天体物理学等许多领域中，它们仍然是一个非常重要的工具，在这些领域中，您确实需要长时间积分，您需要足够快地解决这些积分，而不会产生荒谬的准确性。

我在哪里可以找到辛积分器？存在什么样的辛积分器？

通常有两类辛积分器。有辛 Runge-Kutta 积分器（即上面示例中显示的积分器）和隐式 Runge-Kutta 方法，它们具有辛属性。正如@origimbo 所提到的，辛龙格-库塔积分器要求您为它们提供分区结构，以便它们可以分别处理位置和动量部分。然而，与评论相反，隐式 Runge-Kutta 方法是辛的，不需要这样做，而是需要求解非线性系统。这还不错，因为如果系统是非刚性的，这个非线性系统可以通过函数迭代或安德森加速来解决，但是辛 RK 方法仍然应该是效率的首选（它'

也就是说，odeint 没有来自这两个系列的方法，因此如果您正在寻找辛积分器，它不是一个好的选择。在 Fortran 中，Hairer 的站点有一小部分可供您使用。Mathematica 有一些内置的. GSL ODE 求解器具有隐式 RK 高斯点积分器，IIRC 是辛积分器，但这是使用 GSL 方法的唯一原因。

但是最全面的辛积分器集可以在 Julia 的 DifferentialEquations.jl 中找到（回想一下，上面的笔记本使用了它）。在此页面上可以找到可用的辛 Runge-Kutta 方法列表，您会注意到隐式中点方法也是辛的（隐式 Runge-Kutta Trapezoid 方法被认为是“几乎是辛的”，因为它是可逆的）。它不仅拥有最大的方法集，而且是开源的（您可以看到高级语言的代码及其测试）并且有很多基准测试。使用它来解决物理问题的一个很好的入门笔记本是这个教程笔记本. 但当然建议您通过第一个 ODE 教程开始使用该软件包。

一般来说，您可以在此博客文章中找到数值微分方程组的详细分析。它非常详细，但是由于它必须涵盖很多主题，因此每个主题的细节都比这要少，因此请随时要求以任何方式对其进行扩展。

为了补充 Chris Rackauckas 的回答，陈述一些数学上的废话以及你几乎肯定知道的一些东西，如果有坐标描述，动力系统就是哈密​​顿量p和q和一个功能，H(p,q)这样

dqdt=+∂H∂p
和
dpdt=−∂H∂q.
 这个动作保存了H沿着轨迹，但它还有一个额外的属性，即如果我们定义一个映射
p(t),q(t)=ϕt(p(t0),q(t0))
 那么这个映射保留了两种形式dp∧dq. 对于一个问题p和q是一维的，您可以将其视为相空间上闭合曲线内的区域是守恒的。这确保了各种良好的稳定性属性，因为轨迹的“球”必须彼此“靠近”。

在数字方面，辛积分器以相同的方式起作用，也保留了这个区域/两个形式。反过来，这意味着存在一个保守的“数值哈密顿量”（可能与确切的不一样）。请注意，稳定性与准确性不同，因此辛方法的大部分优点是在长时间积分时出现的（例如，您的方法可能会迅速将卫星放在地球的错误一侧，但绝不允许它衰减到它）。

其它你可能感兴趣的问题
上一篇发表强调软件的论文的场所 下一篇如果我希望我的代码将来在千万亿级机器上运行，我应该投资哪些编程范式？