
47 Numerical stability of computations *************************************** You may hear sometimes the following naive argument: Why spend time on all those small floating point errors? Even REAL*4 gives us at least 7 significant digits, and surely that is enough for most applications, and there is the REAL*8, and that must surely suffice. One counterargument is that errors may propagate and grow rapidly, and since "small" roundoff errors are inevitable, there is always the danger that our computations will produce garbage. Introduction  Problems that have an inherent tendency for "error amplification" are called 'illconditioned', they require using multiprecision arithmetic packages and tight error control. A classical example of an illconditioned problem is a certain polynomial of degree 20, upon multiplying one of its coefficients (the one that multiplies the 19th power) by a factor of 1.0000000005 the location of most of the roots is completely changed. In the above example, simply entering the input data (the coefficients of the polynomial) to the computer may be enought to trash the whole calculation, as the conversion from the external decimal representation to the internal binary one may introduce too large errors. A 'wellconditioned' problem may have an error amplifying algorithm, such an algorithm is said to be 'unstable'. In simple terms, an unstable algorithm transforms a "small" change in the input quantities into a "large" change in the output quantities. Common mechanisms for error amplification are: I) Evaluating a rapidly changing quantity (large derivative), the errors are amplified by the magnitude of the derivative II) Mixing/alternating of different solutions to the same problem, e.g. illconditioned polynomials, ordinary differential equations with another solution that increases rapidly III) Iterative evaluation, may increase an initial error repeatedly via other mechanisms Sometimes an unstable algorithm may be modified to produce a stable one. ++  COMPUTATIONAL ERRORS MAKES YOU SOLVE A DIFFERENT PROBLEM  ++ A very simple example for mechanism I  square root  This problem is simple to the point of being artificial, but don't get it wrong, such situations do occur in real programs. We compute the square root of a small number DX, in two ways: after adding and subtracting 1.0 and directly. C For DEC machines use: SQUARE1 = 9.0E08, SQUARE2 = 5.76E08 C For IEEE machines use: SQUARE1 = 9.0E08, SQUARE2 = 5.76E08 C PROGRAM SQROOT C  REAL * X, * SQRT, * DIRECT_SQRT, * ERROR, * SQUARE1, SQUARE2 C  PARAMETER( * SQUARE1 = 9.0E08, * SQUARE2 = 5.76E08) C  WRITE(*,*) X = 1.0 + SQUARE1 WRITE(*,*) ' We choose DX = ', SQUARE1 WRITE(*,*) ' X = 1.0 + DX = ', X WRITE(*,*) ' X  1.0 = ', (X  1.0) C  SQRT = (X  1.0) ** 0.5 WRITE(*,*) ' (X  1.0) ** 0.5 = ', SQRT DIRECT_SQRT = (9.0E08 ** 0.5) WRITE(*,*) ' Direct root of 9.0E08 = ', DIRECT_SQRT ERROR = (SQRT  DIRECT_SQRT) / DIRECT_SQRT WRITE(*,*) ' Relative error = ', ERROR WRITE(*,*) C  WRITE(*,*) X = 1.0 + SQUARE2 WRITE(*,*) ' We choose DX = ', SQUARE2 WRITE(*,*) ' X = 1.0 + DX = ', X WRITE(*,*) ' X  1.0 = ', (X  1.0) C  SQRT = (X  1.0) ** 0.5 WRITE(*,*) ' (X  1.0) ** 0.5 = ', SQRT DIRECT_SQRT = (5.76E08 ** 0.5) WRITE(*,*) ' Direct root of 5.76E08 = ', DIRECT_SQRT ERROR = (SQRT  DIRECT_SQRT) / DIRECT_SQRT WRITE(*,*) ' Relative error = ', ERROR WRITE(*,*) C  END In the first case we get a relative error of 15%, and in the second we have an error of 100%, unbelievable, isn't it? The large relative errors are a consequence of an unstable algorithm. In the example program we computed the square root of a small number, we know that: y = x ** 0.5 y' = 0.5 / (x ** 0.5) The derivative goes infinite when we approach x = 0.0, so clearly a small change there in x will produce a LARGE change in y. Of course, in our finite arithmetic we can't get infinitely close to x = 0.0, near that point we have x = n * Xmin Let's stabilize our simple algorithm, in our example it is enough to replace: X = 1.0 + 9.0E08 SQRT = (X  1.0) ** 0.5 By the more reasonable: SQRT = 9.0E08 ** 0.5 Now we get the correct result, why? In the first case we introduced a large roundoff error in computing 1.0 + 9.0E08, because a small number was added to 1.0. Remember that the "density of points" near 1.0 is much lower than near 0.0, or said differently the "granularity" is higher. Subtracting the 1.0 again didn't help, on the contrary, and the large roundoff error stayed and produced a wrong square root. Of course you can't always transform an unstable algorithm to a stable one by simple arithmetic juggling, sometimes the problem is deeper, and you have to use a better algorithm. A more realistic example for mechanism I  the quadratic equation  An example program for computing error amplification  The following procedure (STDRVR) can track error propagation through a routine with a certain arrangement of arguments. The tested routine should have REAL arguments, the first is the "result", and 1 to 3 "variables" follow it. Such "variable argument list" routines usually work quite nicely, if however they don't please notify this guide. In this example (see the call to STDRVR) the routine COMPUTE is run 20 times, the arguments X1, X2 are varied by random increments on the order of the common REAL machine precision, and the resulting values are monitored. Try giving small positive values to X1 and X2, and note the large error amplifications that are generated. PROGRAM DRIVER C  REAL * X1, X2 C  EXTERNAL * COMPUTE C ================================================================== WRITE (*,'(1X,A,$)') 'ENTER X1, X2: ' READ (*,*) X1, X2 C  CALL STDRVR(COMPUTE, 2.0E7, 20, 2, X1, X2) C  END SUBROUTINE COMPUTE(Y, X1, X2) C  REAL * Y, X1, X2 C  Y = 1.0E+03 * (X1 ** 0.5) + 1.0E+03 / X2 C  RETURN END SUBROUTINE STDRVR(FUNC, ERRABS, NITER, NARGS, X1, X2, X3) C  INTEGER * NITER, NARGS, * SEED, COUNT, I C  REAL * ERRABS, X1, X2, X3, Y, * EPS, MINEPS, XIN(3), MIN, MAX, SUM, * RELERR, ABSERR C  PARAMETER( * MINEPS = 2.0E07) C  SAVE * COUNT C  DATA * COUNT /0/ C  EXTERNAL * FUNC C ================================================================== IF ((NARGS .LT. 1) .OR. (NARGS .GT. 3)) * STOP 'STDRVR: Wrong number of args ' C  COUNT = COUNT + 1 WRITE (*,*) WRITE (*,*) ' Data collection run # ', COUNT C  SEED = 1234567 C  IF (NARGS .GE. 1) XIN(1) = X1 IF (NARGS .GE. 2) XIN(2) = X2 IF (NARGS .GE. 3) XIN(3) = X3 WRITE (*,*) ' Base values = ', (XIN(I), I = 1, NARGS) C  IF (NARGS .EQ. 1) CALL FUNC(Y, XIN(1)) IF (NARGS .EQ. 2) CALL FUNC(Y, XIN(1), XIN(2)) IF (NARGS .EQ. 3) CALL FUNC(Y, XIN(1), XIN(2), XIN(3)) C  MIN = Y MAX = Y SUM = Y C  IF (ERRABS .GE. MINEPS) THEN EPS = ERRABS WRITE (*,*) ' Initial error = ', EPS ELSE EPS = MINEPS WRITE (*,*) ' Initial error raised to ', MINEPS ENDIF C  DO I = 2, NITER IF (NARGS .GE. 1) XIN(1) = X1 + EPS * 2.0 * (RAN(SEED)  0.5) IF (NARGS .GE. 2) XIN(2) = X2 + EPS * 2.0 * (RAN(SEED)  0.5) IF (NARGS .GE. 3) XIN(3) = X3 + EPS * 2.0 * (RAN(SEED)  0.5) C  IF (NARGS .EQ. 1) CALL FUNC(Y, XIN(1)) IF (NARGS .EQ. 2) CALL FUNC(Y, XIN(1), XIN(2)) IF (NARGS .EQ. 3) CALL FUNC(Y, XIN(1), XIN(2), XIN(3)) C  IF (Y .LT. MIN) MIN = Y IF (Y .GT. MAX) MAX = Y SUM = SUM + Y ENDDO C  ABSERR = ABS(MAX  MIN) RELERR = NITER * ABSERR / SUM C  WRITE (*,'(2X,70(1H))') WRITE (*,*) ' Number of points = ', NITER WRITE (*,*) ' Max = ', MAX, ' Min = ', MIN WRITE (*,*) ' Relative error = ', 100.0 * RELERR, ' percent ' WRITE (*,*) ' Error amplification = ', ABSERR / EPS WRITE (*,'(2X,70(1H=))') C  RETURN END
