154
6 Computing Integrals and Testing Code
Rounding Errors When we compute with real numbers, these numbers are
inaccurately represented on the computer, and arithmetic operations with inaccurate
numbers lead to small rounding errors in the final results. Depending on the type of
numerical algorithm, the rounding errors may or may not accumulate.
Testing with a Tolerance If we cannot make tests like 0.1 + 0.2 == 0.3, what
should we then do? The answer is that we must accept some small inaccuracy and
make a test with a tolerance. Here is the recipe:
In [1]: a = 0.1; b = 0.2; expected = 0.3
In [2]: computed = a + b
In [3]: diff = abs(expected - computed)
In [4]: tol = 1E-15
In [5]: diff < tol
Out[5]: True
Here we have set the tolerance for comparison to 10 −15 , but calculating 0.3 -
(0.1 + 0.2) shows that it equals -5.55e-17, so a lower tolerance could be used
in this particular example. However, in other calculations we have little idea about
how accurate the answer is (there could be accumulation of rounding errors in more
complicated algorithms), so 10 −15 or 10 −14 are robust values. As we demonstrate
below, these tolerances depend on the magnitude of the numbers in the calculations.
Absolute and Relative Differences Doing an experiment with 10 k + 0.3 − (10 k +
0.1 + 0.2) for k = 1, . . . , 10 shows that the answer (which should be zero) is
around 10 16−k . This means that the tolerance must be larger if we compute with
larger numbers. Setting a proper tolerance therefore requires some experiments to
see what level of accuracy one can expect. A way out of this difficulty is to work
with relative instead of absolute differences. In a relative difference we divide by
one of the operands, e.g.,
a = 10
k
+ 0.3, b = (10
k
+ 0.1 + 0.2), c =
a − b
a
.
Computing this c for various k shows a value around 10 −16 . A safer procedure is
thus to use relative differences.
We may exemplify this in a quick session, using k = 10,
In [1]: a = 10**10 + 0.3
In [2]: b = 10**10 + 0.1 + 0.2
In [3]: diff = a-b
In [4]: diff
Out[4]: -1.9073486328125e-06
In [5]: rel_diff = (a-b)/a
6 Computing Integrals and Testing Code
Rounding Errors When we compute with real numbers, these numbers are
inaccurately represented on the computer, and arithmetic operations with inaccurate
numbers lead to small rounding errors in the final results. Depending on the type of
numerical algorithm, the rounding errors may or may not accumulate.
Testing with a Tolerance If we cannot make tests like 0.1 + 0.2 == 0.3, what
should we then do? The answer is that we must accept some small inaccuracy and
make a test with a tolerance. Here is the recipe:
In [1]: a = 0.1; b = 0.2; expected = 0.3
In [2]: computed = a + b
In [3]: diff = abs(expected - computed)
In [4]: tol = 1E-15
In [5]: diff < tol
Out[5]: True
Here we have set the tolerance for comparison to 10 −15 , but calculating 0.3 -
(0.1 + 0.2) shows that it equals -5.55e-17, so a lower tolerance could be used
in this particular example. However, in other calculations we have little idea about
how accurate the answer is (there could be accumulation of rounding errors in more
complicated algorithms), so 10 −15 or 10 −14 are robust values. As we demonstrate
below, these tolerances depend on the magnitude of the numbers in the calculations.
Absolute and Relative Differences Doing an experiment with 10 k + 0.3 − (10 k +
0.1 + 0.2) for k = 1, . . . , 10 shows that the answer (which should be zero) is
around 10 16−k . This means that the tolerance must be larger if we compute with
larger numbers. Setting a proper tolerance therefore requires some experiments to
see what level of accuracy one can expect. A way out of this difficulty is to work
with relative instead of absolute differences. In a relative difference we divide by
one of the operands, e.g.,
a = 10
k
+ 0.3, b = (10
k
+ 0.1 + 0.2), c =
a − b
a
.
Computing this c for various k shows a value around 10 −16 . A safer procedure is
thus to use relative differences.
We may exemplify this in a quick session, using k = 10,
In [1]: a = 10**10 + 0.3
In [2]: b = 10**10 + 0.1 + 0.2
In [3]: diff = a-b
In [4]: diff
Out[4]: -1.9073486328125e-06
In [5]: rel_diff = (a-b)/a