[More Rounding-off Error Visualizations]

Jean-François COLONNA

[Contact me]

france telecom, France Telecom R&D

[

[

[

(CMAP28 WWW site: this page was created on 10/21/1996 and last updated on 11/14/2023 17:54:46 -CET-)

(published in

For chaotic processes, the sensivity to initial conditions is a well known phenomenon which is said to be the cause of the impossibility of long term forecast. In my point of view, there is a phenomenon of greater importance: the sensivity to numerical accuracy as already described in [01][02]. The rounding-off errors are inherent to the computer architecture and if for daily use 32 or 64 bit accuracy is enough, fixed-size accuracy will not be sufficient when running non linear programs. Moreover in this context the associative property of the multiplicative operator is lost and unfortunately, for most values A, B and C:

(this is also true for the addition, but the effect on the multiplication is worse). Let's recall this can be shown very easily with the Verhulst dynamics that models the variation of a population of individuals by the following iteration:(AxB)xC # Ax(BxC)

2 X = (R+1)X - RX n n-1 n-1('R' denotes the growing rate. When R > 2.57 this shows a so-called "chaotic" behaviour.) A very simple short C program in double precision can be written in order to compute Xn for any value of n; unfortunately there are many ways to get the results as shown in the following two tables obtained on an IBM ES9000 and an IBM RS6000 (with R = 3.0):

IBM ES9000:

(R+1)X-R(XX) (R+1)X-(RX)X ((R+1)-(RX))X RX+(1-(RX))X X+R(X-(XX))

X(00) = 0.500000 0.500000 0.500000 0.500000 0.500000 X(10) = 0.384631 0.384631 0.384631 0.384631 0.384631 X(20) = 0.418895 0.418895 0.418895 0.418895 0.418895 X(30) = 0.046399 0.046399 0.046399 0.046399 0.046399 X(40) = 0.320185 0.320183 0.320188 0.320182 0.320189 X(50) = 0.063406 0.064521 0.061895 0.064941 0.061244 X(60) = 1.040381 0.846041 0.529794 1.319900 1.214070 X(70) = 0.004104 1.199452 0.873553 0.573637 0.000009 X(80) = 0.108044 0.121414 1.260726 0.395871 0.280590 X(90) = 0.096374 0.089244 0.582157 0.344503 1.023735

IBM RS6000:

(R+1)X-R(XX) (R+1)X-(RX)X ((R+1)-(RX))X RX+(1-(RX))X X+R(X-(XX))

X(00) = 0.500000 0.500000 0.500000 0.500000 0.500000 X(10) = 0.384631 0.384631 0.384631 0.384631 0.384631 X(20) = 0.418895 0.418895 0.418895 0.418895 0.418895 X(30) = 0.046399 0.046399 0.046399 0.046399 0.046399 X(40) = 0.320177 0.320184 0.320188 0.320190 0.320189 X(50) = 0.067567 0.063747 0.061859 0.060822 0.061486 X(60) = 0.001145 0.271115 0.616781 0.298613 1.307350 X(70) = 1.296775 1.328462 0.486629 0.938605 1.054669 X(80) = 0.553038 0.817163 1.277151 1.325437 0.617058 X(90) = 0.094852 0.154184 1.174162 0.148151 0.237355

The consequence is that the results obtained may depend on many extra-parameters:

- the way parentheses are used,
- the system release (in particular the compiler's one),
- the compiling options,
- the instruction scheduling (see for example new out of order microprocessors where under certain circumstances, some results obtained could not be reproduced five minutes later and this, without changing anything in the environment...),
- the floating point processor and the internal format of numbers,
- the way floating point constants are translated into binary values,
- etc...

It is noteworthy to recall that the infinite Real number field is mapped on many different finite floating point number sets and that these sets are not themselves fields. One of the consequences is that one cannot study stricto sensu problems depending "heavily" on Real numbers (in the mathematical sense) with a computer: digital Real numbers do not exist! Today, where both fundamental and applied research rely heavily on computers [03], and where numerical techniques and devices are omnipresent, it is more than essential to recall to computer users those problems. To do so, I have chosen to exhibit the N-body problem. It is important for the reader to remember that my purpose here is not to talk about neither deterministic chaos, nor celestial mechanics, nor applied mathematics (this problem involving numerical integration methods); my only purpose here is only to show how dangerous it could be to do multiplications with a computer...

The N-body problem is the following one: N gravitationally interacting bodies are specified with their initial conditions (at time t=0) and we should like to know their individual positions at any time t. Since Kepler [04] and Newton [05] it is known that for N=2 the trajectories are elliptic. Unfortunately, as stated by Henri Poincare [06], for N strictly greater than 2, this problem cannot be studied analytically. Fortunately, in the forties, John Von Neumann and many others gave us the computer that can be used today to study anything... The N-body problem can then be numerically studied using the Newton law:

--> ------> F = m. Gammafor each of the studied body (k):

---> -------> F = m . Gamma k k kwith:

2 ----> d OC -------> k Gamma = --------- k 2 dtand:

2 ----> i=N d OC ----- M k \ i -----> --------- = / G ----------- C C 2 ----- 3 k i dt i=1 |----->| (i#k) | C C | k i(where O and Ck denote respectively the origin of the coordinate system and the position of the number k body). Then, it "suffices" to solve a second order differential equation system in order to compute the instaneous coordinates of each bodies studied. I wrote such a program; the following inputs are needed for each body Ck:

- the initial position {Xk(0),Yk(0),Zk(0)},
- the initial velocity {VXk(0),VYk(0),VZk(0)},
- the mass mk,

Figure displays the case where there are only two bodies (N=2). As expected the trajectory of the second one is an ellipse whose focus is the first one (chosen to be still). Figure displays our actual Solar System whose evolution is computed for about 18 years; its shows non chaotical trajectories (each of them are very close to an ellipse whose focus is the sun). Figure exhibits a very interesting case with four bodies (N=4); as a matter of fact a very heavy one (five times heavier than Jupiter) is located near the Sun in order to induce strong perturbations of the trajectories of the two light bodies. In this case, their trajectories are far from being elliptic or even periodic. Finally the same double precision C program with exactly the same initial conditions is run on three "compatible" computers:

- a HP 9000/755 workstation under
HP-UX A.09.05 with a cc compiler
(dubbed the
*Red*computer), - a VAX9000 mainframe under
ULTRIX 4.4 with a gcc compiler (dubbed
the
*Green*computer), - a SGI Power Challenge M server
under IRIX 6.0.1 with a cc compiler
(dubbed the
*Blue*computer).

Picture synthesis allows us to watch this phenomenon in a very dramatic way. The N-body problem is no exception; all non linear models I have numerically studied exhibit the same behaviour (please, see [01] where many examples are available as mpeg animations). So, what can be done ? To change the numerical method used to solve differential equations is of no help (by the way and paradoxically in these circumstances, higher order methods could be less efficient since they make use of a higher number of multiplications); anyway, the real problem lies at a lower level (the arithmetic used). To increase the accuracy of computers (going from to 64-bits to 128-bits or more), which is far from being economical, due to the nature of the problem, will just shift it some time steps later (after all, accuracy should be a function of speed...). As described in [01][02], I am suggesting a method to help telling which models are subject to this phenomenon. I should like to encourage reader involvement; in classroom settings students must be informed of this problem. Research must be conducted at all levels (hardware and software) to help preventing and solving it (if possible...). Finally a major question will remain unanswered maybe forever: which computer does God use to do the huge amount of computations needed to run our Universe in real time?

- [01] Please, visit the "The Subjectivity of Computers" WWW site.
- [02] "The subjectivity of computers", Jean-François COLONNA, Technical Correspondence, Communications of the ACM, volume 36, numero 8, page 15-18, 08/1993.
- [03] "Scientific Display: a Means of Reconciling Artists and Scientists", Jean-François COLONNA, "Frontiers of Scientific Visualization", Clifford A. Pickover and Stuart K. Tewksbury editors, John Wiley and Sons, New York, 1994.
- [04] "Astronomia Nova", Johannes Kepler, 1609.
- [05] "Philosophiae Naturalis Principia Mathematica", sir Isaac Newton, 1687.
- [06] "Les méthodes nouvelles de la mécanique céleste", Henri Poincaré, 1892-1899, Paris.

Copyright © France Telecom R&D and CMAP (Centre de Mathématiques APpliquées) UMR CNRS 7641 / École polytechnique, Institut Polytechnique de Paris, 1996-2023.