Errors in Computing

The sources of errors in numerical computing are many. They can creep in during your measurement of the physical system being modeled, in the mathematical model created, in the algorithm to compute quantities associated with the math model, and from numerically representing floating point numbers on a computer. Most often they come from all of the above.

When computing something, often there is more than one algorithm available. Although those algorithms may be mathematically equivalent, typically they differ greatly numerically and in efficiency. Studying the first is the realm of traditional numerical analysis, while the second is traditionally a computer science field. We'll primarily handle the second. First, a small matter of endless confusion needs to be discussed:

Absolute versus Relative Errors

If a scalar a is approximated by á (for example, a computed version of a), the absolute error is |a-á| and the relative error is |a-á|/|a|, with the usual disclaimer. You should always bear in mind whether a person is talking about absolute or relative error! Here is a classical example from Charles Van Loan of contrasting those two types of error. Stirling's formula gives an approximation to the factorial function
n! = 1*2*3* ... *n = factorial(n)
where Stirling's formula is

Sn = sqrt(2*pi*n)*(n/e)n ~ n!

No doubt you were told at some time that "as n increases, Stirling's formula asymptotically approaches n!". Check it with a Matlab script 

% Stirling
%
% Prints a table showing error in Stirling's formula for n!
% Probably you should limit the value of n to 15 or less,
% because of the formatting used.
%
% Modification of a routine from Charles Van Loan
%
   close all
   nmax = input('Largest value of n to test: ');
   clc
   home
   disp('  ')
   disp('                                 Stirling          Absolute    Relative')
   disp('   n                n!         Approximation         Error       Error')
   disp('----------------------------------------------------------------------')
   e=exp(1);
   nfact=1;
   for n=1:nmax
      nfact = n*nfact;
      s = sqrt(2*pi*n)*((n/e)^n);
      abserror = abs(nfact - s);
      relerror = abserror/nfact;
      s1 = sprintf('  %2.0f   %14.0f   %17.2f',n,nfact,s);
      s2 = sprintf('   %13.2f    %5.2e',abserror,relerror);
      disp([s1 s2])
   end
Matlab lets you use "sprintf", which works like C's sprintf. It creates a formatted string and assigns it to a variable. The %2.0f indicates a floating point number with two digits and zero digits after the decimal point. The %5.2e means an exponential notation number with two digits after the decimal point and one before. It generates: 
                                 Stirling          Absolute    Relative
   n                n!         Approximation         Error       Error
----------------------------------------------------------------------
   1                1                0.92            0.08    7.79e-02
   2                2                1.92            0.08    4.05e-02
   3                6                5.84            0.16    2.73e-02
   4               24               23.51            0.49    2.06e-02
   5              120              118.02            1.98    1.65e-02
   6              720              710.08            9.92    1.38e-02
   7             5040             4980.40           59.60    1.18e-02
   8            40320            39902.40          417.60    1.04e-02
   9           362880           359536.87         3343.13    9.21e-03
  10          3628800          3598695.62        30104.38    8.30e-03
  11         39916800         39615625.05       301174.95    7.55e-03
  12        479001600        475687486.47      3314113.53    6.92e-03
  13       6227020800       6187239475.19     39781324.81    6.39e-03
  14      87178291200      86661001740.60    517289459.40    5.93e-03
  15    1307674368000    1300430722199.47   7243645800.53    5.54e-03
The absolute error in the formula grows with n, but the relative error decreases. So it is not correct to say "Sn converges to n!".

By the way, can Matlab's computation of 15! be trusted? Let's see: it has 7 even numbers in its product expansion, and 3 multiples of 5, so it should have three zeros on the end, which it does. Since the sum of its digits is 45, it is divisible by 9. Since the alternating sum of its digits is -11, the number is divisible by 11. The number comprised of the last seven digits (4368000) is divisible by 27 = 128. Finally, it agrees with Stirling's formula in order of magnitude and to the first three significant digits. OK, maybe it is right. This is something you should always do: try to validate your numerical results with an independent method, or at least get estimates. Often it is easy to get an idea of what ballpark the numbers should be playing in, e.g., Avagadro's number looking like 3.749853x10-29, or if the computed speed of light in a vacuum turns out to be slower than a campus bus, then you should immediately know to start debugging.


Next: Errors in Floating Point Arithmetic
or toddle on home to: P573