Lab 2: Errors

In Lab 1, we discussed the represenation of binary floating point numbers. One of the key results we discussed was the fact that there are only a finite number of bits, usually 64, used to store a floating point number.

Today, we're going to explore the often surprising consequences of this fact. First, we'll recall how digits are stored: in an IEEE single precision (32 bit) float, a number is represented as

where the bias $B = 127$. For example,

What decimal number does this represent? Work it out on paper!

Adding Floating point numbers

To begin to see how error creeps in just from our finite floating point representation, let's consider adding two numbers in regular, decimal scientific notation with bias $B= 0$. If we had say $7 + 3 \times 10^{-5}$, we can easily see the answer is $7.00003$. If we use double precision numbers (as python does by default), this poses no challenge to the computer either:

However, let's consider doing the same thing, but now we'll assume we have only 4 digits of precision in the mantissa of our numbers. That gives us $7.000$ and $3.000 \times 10^{-5}$. So far, so good. Floating point allows us to represent both of these quantities with the same number of digits of precision.

When we go to add them, the computer will do that by shifting the smaller number such that it has the same exponent of as the larger number, in this case $10^0$. Remember, we can represent the same number many ways--but we must keep only 4 digits of precision:

$3.000 \times 10^{-5} = 0.300 \times 10^{-4} = 0.030 \times 10^{-3} = 0.003 \times 10^{-2}$

You can see where this is going. We shift one more time, and we end up with $0.000 \times 10^{-1}$. This is not good, since this is zero, and we still have one more shift to do! We end up with

$7.000 \times 10^0 + 0.000 \times 10^0$ = $7.000 \times 10^0,$ which is clearly the wrong answer.

Machine Precision

This leads us to a very important definition: there exists a number $\epsilon_m$ such that $1_c + \epsilon_m = 1_c$, where that subscript $c$ means "as represented on the computer" rather than "the mathematical object 1".

Let's find an approximation for $\epsilon_m$ experimentally. In doing so, we're going to meet another kind of loop, called the while loop. A while loop tests an expression at the top, just like an if statement. If the expression is true, then the loop body executes. If it is false, the loop stops. Be careful! While loops may never stop! This is sometimes what you want, but often not. You must make sure that the expression at the top changes!

hint

If your while loop seems to never end, you can go to the Kernel menu in the Jupyter notebook and choose Interrupt Kernel. That should cause your code to stop and an error message to be printed. You can then fix your loop and restart. Any calculation like these ones should be more or less instantaneous. Later in the class, we'll have calculations that will take longer.

So what does this code do? How does it find epsilon?

That expression above was using standard python floating point numbers, which are double precision. We can use numpy to create single precision numbers. Here's how:

Note we have to be careful! Every number in python is double precision unless you say otherwise. Also, python is smart enough to automatically "promote" a variable like eps to double precision if you add, subtract, multiply, or divide by something of double precision. That's why I've carefully set all numbers to be np.float32s in the above cell.

Questions

What is the difference in $\epsilon_m$ between single and double precision? Try to explain what is going on in terms of what you've learned about floating point numbers, comparing them to numbers written in scientific notation with a limited precision.

Subtractive cancellation

By hand, subtract $5.324$ from $5.321$. How many digits of precision did you start with? How many did you end with? What does that tell you about your new answer? Think about these questions for a minute before proceeding.

This is called subtractive cancellation: you are cancelling off all of your more precise digits, and leaving only the least precise one (in this case). The subtraction you compute may be exact itself, but if there is any error in the starting values ($5.324$ and $5.321$) it becomes amplified. The technical term for this is "bad."

Example 1

Calculate $y = \sqrt{x + 1} - \sqrt{x}$ for $x = 10^i$ with $i = 1,2,3,\ldots,20$. Be careful: you want $x$ to be a float, not an int. To do that make sure you use 10.**i Store the answers in a numpy array called y1.

Now do the same thing, but get rid of the subtraction. Multiply $y$ by

You should get a new expression for $y$ with no subtractions. Now, using your new formulation, calculate $y$ again, and store it in a new numpy array y2.

Which one of these is more accurate? If it's not obvious, you probably made a coding mistake.

The less accurate one is called unstable. The errors caused by subtractive cancellation grow as you move to certain values of $x$, in this case, large values.

Now using pyplot, plot the answer with subtraction y1 and the answer without subtraction y2 as a function of x. Code to do this is in the cell below. Look at the code and see the sytle of the plot.

What's wrong with this plot? one thig to note is that the $x$ axis clearly varies by a lot. Using a logarithmic plot is helpful here. Matplotlib allows log plots on the x, y, or both axes easily:

plt.semilogx(x,y)
plt.semilogy(x,y)
plt.loglog(x,y)

But even if you try that you might still not be able to see a clear difference between them. How could you fix this? hint: it's something you learned in 107, if not before that

Fix it, and show a plot that displays how much worse one solution is than the other. Make sure to label all axes on the plot!

Example 2

You need to compute the following

for $\theta \to \pi/2$.

Your first thought is to simply solve the equation for a set of $\theta$, something like $\theta = [\pi/2 - n \epsilon, \pi/2 - (n-1) \epsilon,\pi/2 - (n-2) \epsilon\ldots,\pi/2-\epsilon]$, where $\epsilon < 1$ and $n \simeq 200$. Don't use $\theta = \pi/2$! Ask yourself why not before continuing.

First, set up your $\theta$ variable as a numpy array called theta. Make sure the values get closer to $\theta$. Try experimenting with different values of $\epsilon < 1$.

Now that you have theta, compute $y$ using the above formula.

Now, use your knowledge of trigonometry to find an equivalent form for $y$. Hint: it's a single trig function

numpy has that function built-in. Run it on your theta array, and print the difference between your computed y and the answer given by the numpy function. This is called the "absolute error", $E_{abs} \equiv |x_{true} - x_{calc}|$.

Now compute the relative error (you may have called this the "percent error" before). $E_{rel} \equiv |x_{true} - x_{calc}|/x_{true}$.

Finally, plot the relative error as a function of $\theta$.

Very important note: The true solution is often unknown, and in any real problem is totally unknown. If you knew the true solution, you wouldn't be bothering to run a computer simulation! It's very important for you to remember that the built in numpy function you just calculated is also an approximation, but in this case, it's a much, much better approximation.

Every computational answer is wrong, but some wrong answers are useful!