如何获得更准确的取消

计算科学 Python scipy 龙格库塔 微分方程 准确性
2021-12-20 02:44:44

我会尽量直截了当,所以如果还有什么问题,请告诉我,您需要更多详细信息。

我正在求解几个没有显式耦合的方程,但它们对应的未知变量,比如必须满足微分方程:xy

x˙=x+y,

其中点表示相对于自变量的导数,例如t

的方程是二阶的,因此可以从中得到,并且可以检查上面的方程是否始终满足。然而(见附图),事实证明,无论我使用 SciPy 中的哪种集成方法(已经实现的方法),上述等式在某些时候都不再满足。这是因为xxx˙xy以非常高的精度相互抵消,这似乎无法通过 SciPy 提供的任何方法来实现(我已经通过采用每种方法并尽可能降低绝对和相对公差来检查这一点。在所附的图中,采用的方法是 DOP853,当需要非常低的公差时,它应该非常有用)。

在此处输入图像描述

我的问题是您是否知道任何提高准确性的方法,以便取消变得更加精确(我希望在整个计算过程中都能方便地满足方程)。到目前为止,我更改的唯一参数是相对和绝对公差(当然还有我们可以使用的不同方法)。是否有任何我遗漏的参数可能对此有用?

1个回答

我不确定 Python 库是否可行,因为它们在后台使用 Fortran 并且不容易重新编译,但是 Julia DifferentialEquations.jl JIT 编译器根据您提供的数字类型专门化了求解器。这是一些奇怪类型的演示,例如有理数、MPFR BigFloats 和 ArbFloats(基于 Arb 库)

您可以在 Feagin 收敛图中看到这一点,该图中通过 BigFloats演示了 14 阶方法到1050的准确性。在使用 BigFloats 或 ArbFloats 的 Julia 中,您可以setprecision更改数字类型的精度以获得所需的精度。

虽然Julia 方法与 SciPy 相比非常快(数量级半数),并且即使它们编译为专门处理输入类型,因此会为高精度情况生成特殊的优化代码,但高精度算术仍然相当昂贵你应该记住这一点。专门针对这种高精度范围的积分器将是相当重要的。请注意,如果您想这样做,我可能会推荐Vern9一种多线程外推方法,例如ExtrapolationMidpointDeuflhard(它们将在f调用之间进行多线程,这将随着容差的降低而变得更加重要),或者可能是新的16th order symplectic integratorIRKGL16

此外,如果您需要经过验证的算术,您可以将TaylorIntegration.jl用于具有解决方案浮点精度界限的高阶泰勒方法。