精确Runge-Kutta系数

Question

当使用数值方法(如龙格-库塔)时，浮子在计算机上的有限精度会影响解(布鲁沃尔定律)。

在本论文中，它建议将精确的Runge系数模拟为A=B+C，其中B是精确的机器数，C是一些小的修正。

有人能在实践中解释一下这是怎么回事吗？如果A= 3/10，那么如何确定B和C？

谢谢你的帮助。

Answer 1

在本文中，他们建议对分母为1024的A使用有理逼近。(这意味着A最多有10个重要的非零位)。你有(3/10)*1024 = 307.2，所以B是

B=307/1024 = 0.2998046875，C=A= 0.0001953125

C不能完全表示为IEEE Binary64，最近的浮点数将为

C=1.9531249999999988887775374843.E-4。

在公式中插入这些值(3.1f)

Answer 2

这个技巧可能在2007年提交论文时起了作用，但我认为它不太可能在现代平台上起作用。

在现代x86 ( 32位和64位)处理器上，有两个独立的浮点计算指令集：

• 旧的x87指令(可以追溯到原始的8087协处理器)，它有80位寄存器。
• 最近的SSE指令，它使用与格式相同的寄存器(32位用于float，64位用于double)。

较新的SSE指令通常被现代编译器所青睐，因为它们往往更快，因为它们可以完全流水线化，并且支持诸如SIMD操作之类的花哨操作。然而，在2007年，一些编译器仍然默认只使用x87指令，因为二进制文件随后可以在旧机器上使用(在32位计算机上尤其如此)。

80位寄存器支持高达64位的意义，比64位double的53位意义多11位.其思想是，您可以潜在地减少中间舍入错误，在这种情况下，您可以利用该错误。

考虑他们问题的一个更简单的版本:计算

Y = A*X

按照他们的建议，通过将A分解为B+C，B只有10个重要位。然后手术

B*X

不会产生任何舍入错误，因为它最多有63位重要位。全计算

Y = B*X + C*X

将给你的结果几乎完全64位的准确性。

如果没有扩展的精度，B*X通常会产生与直接计算A*X大小大致相同的舍入误差(除非X本身存储的精度降低)。

现在这听起来很棒:您可能想知道SSE指令为什么要去掉这个吗？不幸的是，它是不可预测的:在某些情况下，编译器会安排它以使其工作，但在另一些情况下，它需要将寄存器“溢出”到内存中，在这种情况下，您将失去这种额外的精度。这反过来又会给出一些奇怪的结果，比如将x+y == x+y等操作计算为false，这取决于各个操作的评估时间。

然而，并不是所有的东西都丢失了！如果您有一台最近的机器，您可能可以利用融合乘积(fma)操作来获得更高的精确度。在这种情况下，它看起来就像

Y = fma(B,X,C*X)