Examples of Pitfalls in Mathematical Computation

Introduction

The following examples shows that doing mathematical computations is not simply a matter of programming the standard definitions from a mathematics textbook. The examples are from George Forsythe's, Pitfalls in Computation, or Why a Math Book isn't Enough, American Math Monthly, 1970.

Checking equation solutions is difficult

It is difficult to check if an answer is correct. Consider the following pair of equations. 

   Equation 1   0.780*x + 0.563*y − 0.217 = 0
   Equation 2   0.913*x + 0.659*y − 0.254 = 0

You are given the following solutions and must determine which of them is the best solution. 

   Solution 1  x1 = 0.999, y1 = -1.001
   Solution 2  x2 = 0.341, y2 = −0.087

Of course you substitute the solutions into the equations and get the following values for the left hand sides.

Substitute in Equation 1
Substitute in Equation 2

Looks like the second solution is better but the true answer is (1, -1) so the first solution is actually closer to the correct answer.

Simply substituting into equations is not enough to identify which of multiple solutions is the better one!

Commutativity does not hold for computer arithmetic

First add and subtract alternate reciprocals of the integers from 1 to 10.
x = +1.0 − 1.0/2 + 1.0/3 − 1.0/4 + 1.0/5 − 1.0/6 + 1.0/7 − 1.0/8 + 1.0/9 − 1.0/10 =

Now compute the same as above but in the reverse order.
y = −1.0/10 + 1.0/9 − 1.0/8 + 1.0/7 − 1.0/6 + 1.0/5 − 1.0/4 + 1.0/3 − 1.0/2 + 1.0 =

By mathematics x and y are equal but note the diffence in the last digit. To emphasze the difference lets do a bit of identical arithmetic to x and y.
x *= 1.0e16; x /= 3.0; and y *= 1.0e16; y /= 3.0;

By mathematics we expect x and y should be the same so the following values for a and b should be zero.
a = x − y =
b = y − x =

But note that a and b are of opposite sign. From here on any arithmetic done on a and b could diverge drastically with an if statement comparing a (or b) to 0, with the then and else phrases doing different arithmetic.

Computing a series is not straight forward — ex

ex is defined by the infinite series 1 + x + x2/2! + x3/3! + x4/4! ...

We cannot compute all the terms in an infinite series. Instead we compute until the next term is within epsilon (say 10−16) of the current result.

The series works well with positive x but with negative x the computed answer is more and more incorrect the more negative x becomes. For example consider computing e−100 with epsilon = 10−16 =

Using Javascript Math.E**−100 the result =

When x is negative the correct algorithm to use is 1/e−x =

Computing the roots of quadratic equations

A quadratic equation is described by the expression: a*x**2 + b*x + c = 0

The standard method of computing the roots is given by: x = (−b +/− sqrt(b**2 − 4*a*c)) / (2*a)

An improved method of computing the roots is given by the following. 

        if (b > 0) { root1 = (2*c)/(−b − sqrt(b**2 − 4*a*c));
        root2 = (−b − sqrt(b**2 − 4*a*c))/(2*a); }
        else { root1 = (−b + sqrt(b**2 − 4*a*c))/(2*a);
        root2 = (2*c)/(−b + sqrt(b**2 − 4*a*c));
*** Example 1: Given the quadratic equation coefficients: a = 1.0; b = −1.0e10; c = 1.0
The roots are the following
Using the standard algorithm root 1 is root 2 is
Using the improved algorithm root 1 is root 2 is

The improved algorithm gives the correct result for root 2 but no digits are correct in the standard algorithm.

*** Example 2: Given the quadratic equation coefficients: a = 1e-30; b = -1e30; c = 1e30
The roots are the following
Using the standard algorithm root 1 is root 2 is
Using the improved algorithm root 1 is root 2 is

The improved algorithm gives the correct result for root 2 but no digits are correct in the standard algorithm.