In floating point arithmetic, why does numerical imprecision result from adding a small term to a difference of large terms?

Question

I have been reading the book Computer Simulation of Liquids by Allen and Tildesley. Starting on page 71, the authors discuss the various algorithms that are used to integrate Newton's equations of motion in molecular dynamics (MD) simulations. Starting on page 78, the authors discuss the Verlet algorithm, which is perhaps the canonical integration algorithm in MD. They state:

Perhaps the most widely used method of integrating the equations of motion is that initially adopted by Verlet (1967) and attributed to Stormer (Gear 1971). This method is a direct solution of the second-order equation $m_i \ddot{\textbf{r}}_i = \textbf{f}_i$. The method is based on postions $\textbf{r}(t)$, accelerations $\textbf{a}(t)$, and the positions $\textbf{r}(t - \delta t)$ from the previous step. The equation for advancing the positions reads as follows:

$$\tag{3.14}\textbf{r}(t + \delta t) = 2\textbf{r}(t) - \textbf{r}(t - \delta t) + \delta t^2 \textbf{a}(t).$$

There are several points to note about eqn (3.14). It will be seen that the velocities do not appear at all. They have been eliminated by addition of the equations obtained by Taylor expansion about $\textbf{r}(t)$:

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

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

Then, later (on page 80), the authors state:

Against the Verlet algorithm, ... the form of the algorithm may needlessly introduce some numerical imprecision. This arises because, in eqn (3.14), a small term ($\mathcal{O}(\delta t^2)$) is added to a difference of large terms ($\mathcal{O}(\delta t^0)$), in order to generate the trajectory.

I guess that the "small term" is $\delta t^2 \textbf{a}(t)$, and the "difference of large terms" is $2\textbf{r}(t) - \textbf{r}(t - \delta t)$.

My question is, why does numerical imprecision result from adding a small term to a difference of large terms?

I am interested in a rather basic, conceptual reason, since I am not familiar at all with details of floating point arithmetic. Also, do you know of any "overview-type" references (books, articles, or websites) that would introduce me to basic ideas of floating point arithmetic related to this question? Thanks for your time.

Answer 1

Their observation ''the form of the algorithm may needlessly introduce some numerical imprecision'' is correct. But their explanation ''This arises because, in eqn (3.14), a small term ($O(δt^2)$) is added to a difference of large terms ($O(δt^0)$), in order to generate the trajectory.'' is spurious.

The true reason for the slight numerical instability of Verlet's algorithm is that it is only marginally stable, because the difference equation $x_{k+1}=2x_k-x_{k-1}$ (essentially the case where you neglect $a$ in Verlet) has a parasitic solution proportional to $k$, which causes errors introduced to grow linearly in $k$ whereas for a fully stable multistep method applied to a dissipative differential equation, error growth is bounded.

Edit: Note that the book is about numerical simulation of molecular dynamics, and for getting a reasonable accuracy of the resulting expectations one needs a huge number$N$ of steps, as the accuracy scales with $O(N^{-1/2})$ only. (Often the time step is in the picoseconds to follow the intrinsic oscillation scale. But biologically relevant time scales are in the milliseconds or larger ($N\sim 10^9$), though usually one doesn't compute that far.)

For more details, see http://en.wikipedia.org/wiki/Linear_multistep_method#Stability_and_convergence

Answer 2

If you're looking for a good introduction, I would suggest David Goldberg's What Every Computer Scientist Should Know About Floating Point Arithmetic. It may be a bit too detailed, but it's available online for free.

If you've got a good library, I'd suggest Michael Overton's Numerical Computing with IEEE Floating Point Arithmetic, or the first few chapters of Nick Higham's Accuracy and Stability of Numerical Algorithms.

What Allen and Tildesley are referring to specifically is numerical Cancellation. The short of it is that if you have, say, only three digits and you subtract 100 from 101, you get 1.00 (in three digits). The number looks like it is accurate to three digits, but actually, only the first digit is true and the trailing .00 are garbage. Why? Well, 100 and 101 are only inexact representations of, say 100.12345 and 101.4321, but you could only store them as three-digit numbers.

Answer 3

To apply Pedro's example to the equation $(3.14)$, assume your variables are stored with the following values:

$$ \mathbf{r}(t) = 101 $$ $$ \mathbf{r}(t - \delta t) = 100 $$ $$ \delta t^2 \mathbf{a}(t) = 1.49 $$

From $(3.14)$ it should follow that

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

but, since we can only use three digits, the result becomes truncated to

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

This error will propagate, so that after 20 steps, assuming $\mathbf{a}(t)$ remains unchanged, you get $\mathbf{r}(t + 20\delta t) = 331$ instead of $433.90$,

Answer 4

Pedro already gives the important fact, namely cancellation. The point is that every number you compute with has an associated accuracy; for example, a single precision floating point number can only represent things up to approximately 8 digits of accuracy. If you have two numbers that are almost exactly the same but differ in the 7th digit, then the difference will again be an 8-digit single precision floating point number and it looks like it is accurate to 8 digits, but in reality only the first 1 or 2 digits are accurate because the quantities you computed it from are not accurate beyond this first 1 or 2 digits of the difference.

Now, the book you cite is from 1989. Back then, computations were most often done in single precision and round-off and cancellation were serious problems. Today, most computations are done using double precision with 16 digits of accuracy, and this is far less a problem today than it used to be. I think it is worthwhile reading the paragraphs you cite with a grain of salt and take them in the context of their time.