计算科学 - 在浮点算术中，为什么将一个小项加到一个大项的差上会导致数值不精确？ - 吾爱随笔录

在浮点算术中，为什么将一个小项加到一个大项的差上会导致数值不精确？

计算科学浮点计算机算术

2021-11-28 01:15:38

我一直在阅读Allen 和 Tildesley的计算机模拟液体一书。从第 71 页开始，作者讨论了用于在分子动力学 (MD) 模拟中整合牛顿运动方程的各种算法。从第 78 页开始，作者讨论了 Verlet 算法，它可能是 MD 中的规范积分算法。他们说：

也许最广泛使用的整合运动方程的方法是 Verlet (1967) 最初采用的方法，并归因于 Stormer (Gear 1971)。该方法是二阶方程的直接解 $m_i \ddot{\textbf{r}}_i = \textbf{f}_i$ . 该方法基于职位 $\textbf{r}(t)$ , 加速度 $\textbf{a}(t)$ , 和位置 $\textbf{r}(t - \delta t)$ 从上一步开始。推进职位的方程式如下：

$\begin{matrix} (3.14) & r (t + δ t) = 2 r (t) - r (t - δ t) + δ t^{2} a (t) . \end{matrix}$ $\tag{3.14}\textbf{r}(t + \delta t) = 2\textbf{r}(t) - \textbf{r}(t - \delta t) + \delta t^2 \textbf{a}(t).$

关于 eqn (3.14) 有几点需要注意。可以看出，速度根本没有出现。它们已通过添加由泰勒展开得到的方程来消除 $\textbf{r}(t)$ ：

$r (t + δ t) = r (t) + δ t v (t) + (1 / 2) δ t^{2} a (t) + . . .$ $\textbf{r}(t + \delta t) = \textbf{r}(t) + \delta t \textbf{v}(t) + (1/2) \delta t^2 \textbf{a}(t) + ...$

$\begin{matrix} (3.15) & r (t - δ t) = r (t) - δ t v (t) + (1 / 2) δ t^{2} a (t) - . . . . \end{matrix}$ $\tag{3.15}\textbf{r}(t - \delta t) = \textbf{r}(t) - \delta t \textbf{v}(t) + (1/2) \delta t^2 \textbf{a}(t) - ... .$

然后，后来（在第 80 页），作者说：

与 Verlet 算法相比，……算法的形式可能会不必要地引入一些数值不精确性。这是因为，在 eqn (3.14) 中，一个小项 ( $\mathcal{O}(\delta t^2)$ ) 加到大项的差上 ( $\mathcal{O}(\delta t^0)$ )，以生成轨迹。

我想“小项”是 $\delta t^2 \textbf{a}(t)$ ，并且“大项的差”是 $2\textbf{r}(t) - \textbf{r}(t - \delta t)$ .

我的问题是，为什么将小项添加到大项的差异会导致数值不精确？

我对一个相当基本的概念性原因感兴趣，因为我根本不熟悉浮点运算的细节。另外，您是否知道任何“概述类型”参考资料（书籍、文章或网站）可以向我介绍与此问题相关的浮点运算的基本概念？谢谢你的时间。

4个回答

如果您正在寻找一个好的介绍，我建议您阅读 David Goldberg 的What Every Computer Scientist Should Know About Floating Point Arithmetic。它可能有点过于详细，但它可以免费在线获得。

如果你有一个好的库，我建议 Michael Overton 的Numerical Computing with IEEE Floating Point Arithmetic或 Nick Higham 的Accuracy and Stability of Numerical Algorithms的前几章。

Allen 和 Tildesley 具体指的是数值取消。简而言之，如果你只有三位数，然后100从减去101，你会得到1.00（三位数）。这个数字看起来精确到三位数，但实际上，只有第一位是真实的，尾随.00是垃圾。为什么？好吧，100and101只是和的不精确表示100.12345，101.4321但您只能将它们存储为三位数字。

他们的观察“算法的形式可能会不必要地引入一些数值不精确”是正确的。但是他们的解释“这是因为，在 eqn (3.14) 中，一个小项 ( ) 被添加到一个大项的差异 ( ) 中，以生成轨迹。 ''是虚假的。 $O(δt^2)$ $O(δt^0)$

Verlet 算法的轻微数值不稳定的真正原因是它只是勉强稳定，因为差分方程（基本上是在 Verlet中忽略成比例的寄生解，这导致引入的误差在中线性增长，而对于应用于耗散微分方程的完全稳定的多步方法，误差增长是有界的。 $x_{k+1}=2x_k-x_{k-1}$ $a$ $k$ $k$

编辑：请注意，这本书是关于分子动力学的数值模拟，为了获得合理的预期精度，需要大量步，因为精度仅成比例. （通常时间步长以皮秒为单位，以遵循固有的振荡尺度。但生物学相关的时间尺度以毫秒或更大（）为单位，尽管通常不会计算那么远。） $N$ $O(N^{-1/2})$ $N\sim 10^9$

有关更多详细信息，请参阅 http://en.wikipedia.org/wiki/Linear_multistep_method#Stability_and_convergence

要将 Pedro 的示例应用于等式，假设您的变量存储有以下值： $(3.14)$

r (t) = 101

$\mathbf{r}(t) = 101$

r (t - δ t) = 100

$\mathbf{r}(t - \delta t) = 100$

δ t^{2} a (t) = 1.49

$\delta t^2 \mathbf{a}(t) = 1.49$

从应该遵循 $(3.14)$

r (t + δ t) = 103.49

$\mathbf{r}(t + \delta t) = 103.49$

但是，由于我们只能使用三位数字，因此结果会被截断为

r (t + δ t) = 103

$\mathbf{r}(t + \delta t) = 103$

这个错误会传播，所以在 20 步之后，假设保持不变，你会得到而不是， $\mathbf{a}(t)$ $\mathbf{r}(t + 20\delta t) = 331$ $433.90$

佩德罗已经给出了重要的事实，即取消。关键是您计算的每个数字都有相关的准确性。例如，单精度浮点数只能表示大约 8 位精度的事物。如果您有两个几乎完全相同但第 7 位不同的数字，那么差异将再次是 8 位单精度浮点数，看起来它精确到 8 位，但实际上只有第一个1 位或 2 位数字是准确的，因为您计算它的数量超出了前 1 位或 2 位数字的差异。

现在，您引用的书是 1989 年的。当时，计算通常以单精度进行，四舍五入和取消是严重的问题。今天，大多数计算都是使用具有 16 位精度的双精度来完成的，这在今天已经比过去少了很多问题。我认为值得一读你引用的段落，并把它们放在他们所处的时代背景下。

其它你可能感兴趣的问题

上一篇了解迭代方法的“收敛速度” 下一篇在频谱规范枪战中，C 是否比 Fortran 慢（使用 gcc、intel 和其他编译器）？