An important consideration when performing a numerical calculation is the fact that all calculations are performed with finite precision.
There can be VERY different behavior when compared to infinite precision calculations (i.e. doing the calculation by hand).
Consider the numerical representation of $\sqrt{2} = 1.4142\ldots$. As we know, $\sqrt{2}$ is an irrational number and it has an infinite number of decimal digits that never repeat. To store this "simple" number on a computer using only information about its decimal digits we would need an infinite amount of storage.
import numpy as np
a = np.sqrt(2)
print('a = ', a)
print('|2 - a^2| = ', abs(2-a**2))
a = 1.41421356237 |2 - a^2| = 4.4408920985e-16
This is called round-off error.
It results from finite-precision arithmetic.
Some "symbolic" mathematical software (such as Maple and Mathematica) treat $\sqrt{2}$ using its mathematical definition that $(\sqrt{2})^2 = 2$. In this way, the full precision is retained without infinite storage. But we are going to use Python in this class, and the only information we will keep about numbers are their (binary) digits.
Computers store numbers using their binary form.
For our purposes, a binary number will be a sequence of 64 ones and zeros:
$$ 000000010101000001010001010101010101000010101010000001111111111 $$The digits in a binary number are partitioned for different purposes:
This binary sequence in the computer represents numbers as:
$$(-1)^s 2^{c-1023} (1 + f)$$'%.10f' % 0.1
'0.1000000000'
'%.15f' % 0.1
'0.100000000000000'
'%.20f' % 0.1
'0.10000000000000000555'
Consider the binary number
$$ 0 \,11000000001\, 110100001000\ldots $$$= 0.814453125$
Then the number represented by this binary number, in decimal notation is
$$(-1)^s 2^{c-1023} (1 + f) = 1 \cdot 2^{514} (1 + .814453125) \approx 9.7311 \times 10^{154}$$Smallest possible positive number is given by $s = 0, c = 1, f = 0$:
$$ N_{\min} = 2^{-1022} \cdot (1+ 0) $$Any numbers that occur in computation that are less (in absolute value) then $N_{\min}$ give rise to underflow, are treated as zero, which is represented in the machine by setting $c = 0$, $f = 0$ and $s = 0,1$ (two possible choices of zero):
$$ 10000000\ldots \quad \text{and} \quad 0000000\ldots$$np.exp(-800)
0.0
The largest positive number is $s = 0$, $c = 2046$ (this is one less than the maximum value for $c$) and $f = 1 - 2^{-52}$ (max value for $f$):
$$ N_{\max} = 2^{1023} \cdot (2 - 2^{-52})$$Numbers that are larger (in magnitude) than $N_\max$ give overflow, are treated as infinity, and are represented in the machine by putting $c=2047$ (its max value).
np.exp(800)
/Users/petermchale/anaconda/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: overflow encountered in exp """Entry point for launching an IPython kernel.
inf
To get a feeling for the implications of using binary approximations, let's suppose that we approximate numbers using $k$-digit machine numbers:
$$ \pm 0.d_1d_2\cdots d_k \times 10^n. $$And while the mechanism for rounding in a computer can be complicated, we take the convention of chopping to the closest $k$-digit machine number.
Suppose that $p^*$ is an approximation to $p\neq 0$. The absolute error of $p^*$ is given by $|p^* - p|$ and the relative error is given by
$$ \left|\frac{p^*-p}{p} \right|.$$Suppose we approximate $0.89 \times 10^3$ by the 1-digit machine number, $0.9 \times 10^3$. The absolute error is
$$ 0.89 \times 10^3 - 0.9 \times 10^3 = -0.1 \times 10^2 $$and the relative error is
$$ \left | \frac{ 0.89 \times 10^3 - 0.9 \times 10^3}{ 0.89 \times 10^3} \right| = 0.112 \times 10^{-1}. $$Because the numbers we are comparing are large (compared to one), the absolute error is large but the relative error is small.
Compute the roots of $\frac{1}{2} x^2 + 60x + 1 = 0$ to a relative accuracy of $10^{-4}$ using $5$-digit machine arithmetic
We know that
$$ x = {-60 \pm \sqrt{60^2-2}}. $$How accurate is this naive algorithm when implemented on a machine?
from computer_arithmetic import x_plus, x_minus, relative_error
print(relative_error(x_plus))
print(relative_error(x_minus))
print(relative_error(x_plus)/relative_error(x_minus))
0.019858313651939973934458797771105 0.000027762338710156256856180444719621 715.29685806604021135731070209644
The relative error for computing $x_{+}$ is orders of magnitude higher than that for $x_{-}$!
Recall that $x_+ = -b + \sqrt{b^2 - 4ac}$. When $b^2 \gg 4ac$, $x_+$ is the result of subtracting nearly equal numbers.
Subtracting nearly equal numbers can lead to large relative errors on a machine.
As an example, consider the task of computing
$p = 1.01 - 1.00 = 0.01$.
A 2-digit machine would approximate this computation as follows:
$p^* = 1.0 - 1.0 = 0.0$,
leading to a relative error of
$ \left| \frac{p - p^*}{p} \right| = 1. $
Returning to our quadratic-formula example, let's begin by computing what the value of $x_+$ ought to be:
from decimal import getcontext
from computer_arithmetic import minus_b, square_root
getcontext().prec = 32 # don't worry about this; consider it as a way to get an estimate of the true value of x_plus
print('actual value of x_plus = ', minus_b() + square_root())
actual value of x_plus = -0.016668982124708948887200263703
Let's now see what a machine with low precision would compute $x_+$ to be:
getcontext().prec = 5
print('5-digit machine number approximation of x_plus = ', minus_b() + square_root())
5-digit machine number approximation of x_plus = -0.017
actual value of $x_+$ = -0.016668982124708948887200263703
5-digit machine number approximation of $x_+$ = -0.017
What has happened here?
When all is said and done, the computer spits out the approximation above for $x_+$.
What if we used an even lower precision machine...
getcontext().prec = 3
print(minus_b() + square_root())
0
An approximate value of zero, corresponds to a relative error of one. Very bad!
We now understand why the relative accuracy of x_plus is so poor (compared to that of x_minus).
But is there anything we can do about it?
from computer_arithmetic import x_plus_better
print(relative_error(x_plus_better))
print(relative_error(x_minus))
print(relative_error(x_plus_better)/relative_error(x_minus))
0.0000010723684816132654996294022047794 0.000027762338710156256856180444719621 0.038626734325554585367248902589600
Another situation where round-off error comes into play is the evaulation of polynomials. Consider
$$ f(x) = x^3 - 4x^2 + \frac{1}{4} x - 1.$$
We want to evaulate $f(4.16)$ with $3$-digit machine arithmetic.
One approach is to code up an algorithm using this expression for $f(x)$. This leads to 8 operations.
Is there a way to rewrite the polynomial in a way that uses fewer computations?
Rewrite the polynomial as a nested expression:
$$ f(x) = x^3 - 4 x^2 + 1/4x -1 = x(x^2 - 4x + 1/4) - 1\\ = x(x(x-4) + 1/4) -1. $$Because this algorithm has fewer opportunities for making finite precision errors (only 5 operations), we expect it to reduce the relative error...
from computer_arithmetic import f, f_better, relative_error
print(relative_error(f))
print(relative_error(f_better))
0.00067499829114452160707599503891303 0.0000014240470267735154557312198069073