13
$\begingroup$

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.

$\endgroup$
9
$\begingroup$

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

$\endgroup$
10
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ -1: Where is the cancellation that you attribute to the Verlet formula? Typically $\delta t$ is small, which makes $r(\t-\delta t)\approx r(t)$, with no cancellation resulting. Try $r(t)=1$! $\endgroup$ – Arnold Neumaier Sep 3 '12 at 17:15
  • $\begingroup$ @ArnoldNeumaier: Yes, Allen and Tildesley's example doesn't seem to make much sense, I only wanted to provide some reference for problems arising when "a small term [..] is added to a difference of large terms", which is what the OP asked, not whether it is a problem in the given case. $\endgroup$ – Pedro Sep 5 '12 at 7:14
  • $\begingroup$ But adding a small term to a large term is just a rounding error, nothing dangerous at all. Cancellation is when two nearly equal large terms are subtracted to get a tiny term. This becomes a problem only when either the subtracted intermediates are far bigger than the final result of a computation, or when the small intermediate result affected by cancellation is divided by another small element. $\endgroup$ – Arnold Neumaier Sep 5 '12 at 8:11
  • $\begingroup$ @ArnoldNeumaier: As, I think, is quite obvious from my answer, I was referring to the problem of computing the difference, not the sum. $\endgroup$ – Pedro Sep 5 '12 at 13:08
  • 1
    $\begingroup$ @ArnoldNeumaier: Point taken, but I hope you understand that I consider that quite petty for a "-1". $\endgroup$ – Pedro Sep 5 '12 at 13:25
5
$\begingroup$

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$,

$\endgroup$
  • $\begingroup$ But the effect is that large only in 3-digit decimal arithmetic. $\endgroup$ – Arnold Neumaier Sep 3 '12 at 18:49
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ cancellation in double precision arithmetic can be as big a problem as in single precision. A case in point is Gaussian elimination without pivoting, which often gives very poor results due to cancellation, even in double precision. $\endgroup$ – Arnold Neumaier Sep 3 '12 at 9:21
  • $\begingroup$ -1: The Verlet formula typically retains all digits of accuracy, not just 1 or 2 of 8 in single precision. $\endgroup$ – Arnold Neumaier Sep 3 '12 at 17:17
  • $\begingroup$ @ArnoldNeumaier: Sure, you can get the same kind of problems in double precision. All I said is that one doesn't encounter them quite as frequently. $\endgroup$ – Wolfgang Bangerth Sep 4 '12 at 1:36
  • $\begingroup$ If you lose 6 digits three times in a chain of computations you have lost all digits even in double precision. Algorithms suffering from cancellation will usually be poor even in double precision. Verlet's algorithm is different since there is no cancellation but a mild linear growth of errors. Thus the loss of accuracy cannot multiply, making it suitable for much longer integration times. This was surely known to Allen & Tildesley. $\endgroup$ – Arnold Neumaier Sep 4 '12 at 14:10
  • $\begingroup$ Right. But what I mean is that if you have an algorithm without cancellation you still incur an error on the order of 1e-8 in single precision, and if you do 1e8 time steps then you may have a problem even if everything else is exact. 1e8 time steps is an order of magnitude you may have for ODE. On the other hand, in double precision, your inaccuracy in each step is 1e-16 and it would require 1e16 time steps to get a complete loss of accuracy. That is a number of steps that you will not encounter in practice. $\endgroup$ – Wolfgang Bangerth Sep 5 '12 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.