Back to tutorial index

Finite precision point arithmetic


Topics in this lab

Introduction

We all know that when we use the value of pi in Matlab, or compute cos(3.4), we are not getting the "exact" value of $\pi$ or the cosine function, but rather an approximation. In fact, we don't expect to be able to compute any irrational function to all of its digits, if only because we know that such values are nonterminating, non-repeating decimals values. In typical practical situations, our answers will not be significantly affected by the approximations that are made. What may be suprising, however, is that even typical "rational" numbers, such as 0.1 are only approximately represented on the computer.

In this lab, we will explore the number system represented by floating point arithmetic, and discuss some of the consequences for scientific computing. The ideas presented here extend to most modern computing systems, not just Matlab.

clear all
format long e

Back to the top

Examples : Floating point arithmetic

Example 1


Compute the following.

$$ x = 0.1 + 0.1 + 0.1 + ... + 0.1 \qquad \mbox{(10 times)} $$

x = 0;
for i = 1:10,
    x = x + 0.1;
end

What is the difference between $x$ and 1?

fprintf('%g\n',abs(x-1))
1.11022e-16

Example 2


Compute $x$.

$$ x = 2 - 3\left(\frac{4}{3} - 1\right) $$

% Algebra tells us that this should be 1
x = 2 - 3*(4/3 - 1);

What is the difference between $x$ and 1?

fprintf('%g\n',abs(x-1));
2.22045e-16

Example 3


Verify the following exact mathematical expression for the value $a = 0.3$.

$$ 1 + a + a^2 + a^3 + a^5 = \frac{1-a^6}{1-a} $$

a = 0.3;
S_left = sum(a.^(0:5));
S_right = (1-a^6)/(1-a);

How close are the left and right sides of this expression?

fprintf('%g\n',abs(S_left - S_right));
2.22045e-16

Example 4


For the function $f(x) = x$, we can compute the derivative $f'(x)$ exactly using the formula $$ \begin{eqnarray*} f'(x) & = & \frac{f(x+h) - f(x)}{h} = \frac{(x+h) - x}{h} = 1 \end{eqnarray*} $$

Verify this formula for $x = 1$ and $h = 10^{-3}$.

% Compute the derivative of f(x) = x using the secant method
x = 1;
h = 1e-3;
dfdx = ((x + h) - x)/h;    % This should be exactly 1

Do we get $f'(x) = 1$?

fprintf('%g\n',abs(dfdx-1));
1.10134e-13

Back to the top

More examples : Is floating point arithmetic commutative and associative?

You may recall from your first introduction to algebra that we can often re-arrange the order of operations like addition or multiplication. This communative property of these operations allows us to write $$ a + b + c = a + c + b = b + c + a $$

and so on. Unfortunately, this is not always true for floating point arithmetic.

The associative property that we learned in algebra, i.e. that $$ a + (b + c) = (a + c) + b $$

may not always hold either.

Example 5


The order in which we add numbers can matter. Consider these two expressions $$ x = 10^{16} + 1 - 10^{16} $$ $$ y = 10^{16} - 10^{16} + 1 $$ $$ z = 10^{16} - (10^{16} - 1) $$

x = 1e16 + 1 - 1e16;
y = 1e16 - 1e16 + 1;
z = 1e16 - (1e16 - 1);  % Test the associative property

Are $x$, $y$ and $y$ equal?

fprintf('x = %g\n',x);
fprintf('y = %g\n',y);
fprintf('z = %g\n',z);
x = 0
y = 1
z = 0

Example 6

The above can happen for much smaller values as well

x = 1 + 0.1 - 1;
y = 1 - 1 + 0.1;
z = 1 - (1 - 0.1);

Comparing the differences of these three values, we get

fprintf('x-y = %g\n',x-y);
fprintf('y-z = %g\n',y-z);
fprintf('x-z = %g\n',x-z);
x-y = 8.32667e-17
y-z = 2.77556e-17
x-z = 1.11022e-16

Example 7


In this example, we add up 100 random numbers in different orders.

rand('seed',1110);          % Get the same random numbers each time
x = rand(100,1);
xperm = x(randperm(100));   % Permute the values in array x
sum_orig = sum(x);
sum_perm = sum(xperm);

Compare the sum of the original array and the permuted array.

fprintf('%g\n',abs(sum_orig-sum_perm));
1.42109e-14

Back to the top

Examples (modified)

In some of the examples above, a slight change of the constants involved can fix the problems that we saw above. Below are modified versions of some of the above examples.

Example 1 (modified)


Compute the following.

$$ x = 0.125 + 0.125 + 0.125 + ... + 0.125 \qquad \mbox{(8 times)} $$

x = 0;
for i = 1:8,
    x = x + 0.125;
end

What is the difference between $x$ and 1?

fprintf('%g\n',abs(x-1))
0

Example 2 (modified)


Compute $x$.

$$ x = 2 - 2\left(\frac{3}{2} - 1\right) $$

% Algebra tells us that this should be 1
x = 2 - 2*(3/2 - 1);

What is the difference between $x$ and 1?

fprintf('%g\n',abs(x-1));
0

Example 3 (modified)


Verify the following exact mathematical expression for the value $a = 0.0625$.

$$ 1 + a + a^2 + a^3 + a^5 = \frac{1-a^6}{1-a} $$

a = 0.0625;
S_left = sum(a.^(0:5));
S_right = (1-a^6)/(1-a);

How close are the left and right sides of this expression?

fprintf('%g\n',abs(S_left - S_right));
0

Example 4 (modified)


For the function $f(x) = x$, verify the following formula for $x = 1$ and $h = 0.015625$. $$ \begin{eqnarray*} f'(x) & = & \frac{f(x+h) - f(x)}{h} = \frac{(x+h) - x}{h} = 1 \end{eqnarray*} $$

x = 1;
h = 0.015625;
dfdx = ((x + h) - x)/h;    % This should be exactly 1

Do we get $f'(x) = 1$?

fprintf('%g\n',abs(dfdx-1));
0

Back to the top

Machine epsilon

As as way to gauge how close two numbers are to each other, we can use the Matlab function eps(x). This measures the distance to the next representable number. For example,

Example 7

In this example, we want to see if the following mathematical expression holds for all $\x > 0$.

Compute $x$.

$$ 1 + x > 1, \qquad x > 0 $$

x = 1;
while (1 + x > 1)
    x = x/2;
end

Is $x$ equal to 0?

fprintf('x = %8.4e \n',x);
x = 1.1102e-16 

Back to the top

Lab exercises

  1. Can we modify examples 5,6 and 7 above to fix the apparent problems?
  2. Using a loop similar to the one above to compute the smallest value we could add to 1 and still have something greater than 1, try to find the smallest value we can add to $10^6$ and still get something larger than a million?

Compare your answers with the solutions.

Back to the top

Get the code

Do you want to try the above code fragments on your own? Download the Matlab script that produces this page here. (lab_19.m)

Powered by MathJax

Published with MATLAB® 8.5