python的ODE求解器尊重周期性

Posted 2023-03-24

【中文标题】python的ODE求解器尊重周期性

【英文标题】：ODE solver for python respecting periodicity

【发布时间】：2021-04-27 02:23:38

【问题描述】：

我想整合一个控制问题，即 dx/dt = A.f(t) 形式的 ODE，其中 x(t) 是 R^3 中的函数，f(t) 是 R^4 中的函数，A 是 3x4 矩阵。在我的特殊情况下，f(t) = F'(t)，即函数 F 的时间导数。此外，F 是 1 周期的。因此，在区间 [0, 1] 上积分 ODE 应该会再次产生起始位置。然而，scipy.integrate 中的solve_ivp 之类的方法根本不尊重这种周期性（我已经尝试了所有可能的方法，例如RK45, Radau, DOP853, LSODA）。

是否有一种特殊的 ODE 求解器可以高度精确地尊重这种周期性？

【问题讨论】：

您能否举例说明如何违反周期性？您是否设置atol, rtol 并根据这些公差比较值？ // 你检查了导数计算中的错误吗？ // 你能在你的问题中添加一个代码示例吗（这是这样的）？

【参考方案1】：

无论如何，所有算法都会遭受精度损失，因此您很可能永远无法实现精确的周期性。尽管如此，您也可以尝试通过使用参数“atol”和“rtol”来提高积分的精度，这将大致将每个时间步的积分误差保持在 (atol + rtol*y) 以下。您通常可以低至 atol=rtol=1e-14。

【讨论】：

虽然辛积分器靠近能量表面，但它们在保持表面内的相位方面不太精确。检查圆周运动或如您所引用的行星运动。期间仅在方法的顺序内是精确的。 // 在给定任务中，你在哪里识别出哈密顿系统？

@LutzLehmann 确实，我弄错了。我相应地编辑了我的答案，感谢您的评论！

谢谢，atol和rtol确实要选的越低越好。然后，以机器精度获得结果。否则，就曲线 F 的范数而言，有趣的是，错误在 RK45 的 10E-6 和 Radau 的 10E-8 处达到了一个平台。

@mathology 可能会尝试最高阶方法“DOP853”。它是 8 阶的，而 Radau 和 RK45 是 5 阶的。这大致意味着，当时间步长足够小时，DOP853 的误差将比前两种方法低dt^3 倍。当使用非常低的公差 (~1e-14) 时，通常最好使用可能的最高阶。你可能还会在 Github 上找到更高阶的方法。另外，我想知道是否使用恒定时间步长（没有基于容错的自适应时间步长）这样您的周期是整数时间步长可能会提供更好的结果。

对于显式算法（不是 Radau），您可以通过设置 atol=rtol=1e+99（是）和“max_step”=dt 来“破解”您的方式，通过 solve_ivp 施加恒定时间步长.对于隐式算法，这不是一个好主意，因为 atol 和 rtol 用于牛顿终止标准，设置如此高的值可能会导致不合理的行为。