Sources of Error

Error can come from many sources when using applying a numerical method:

• Model/Data Error
• Truncation Error
• Floating Point Error

Goal: Categorize and understand each type of error and explore some simple approaches to analyzing error.

Model and Data Error

Errors in fundamental formulation

• Lotka-Volterra - fractional rabbits, no extinctions, etc.
• Data Error - Inaccuracy in measurement or uncertainties in parameters

Unfortunatley we cannot control model and data error directly but we can use methods that may be more robust in the presense of these types of errors.

Truncation Error

Errors arising from approximating a function with a simpler function (e.g. $sin(x) \approx x$ for $|x| \approx 0$.

Floating Point Error

Errors arising from approximating real numbers with finite-precision numbers and arithmetic.

Basic Definitions

Given a true value of a function $f$ and an approximate solution $\hat{f}$ define:

Absolute Error: $e = |f - \hat{f}|$

Relative Error: $r = \frac{e}{|f|} = \frac{|f - \hat{f}|}{|f|}$

Decimal precision p (number of significant digits) - $p$ is the precision if and only if$$

is true where

Truncation Error and Taylor's Theorem

Taylor's Theorem: Let $f(x) \in C^{m+1}[a,b]$ and $x_0 \in [a,b]$, then for all $x \in (a,b)$ there exists a number $c = c(x)$ that lies between $x_0$ and $x$ such that

where $T_N(x)$ is the Taylor polynomial approximation

and $R_N(x)$ is the residual (the part of the series we left off)

Another way to think about these results involves replacing $x - x_0$ with $\Delta x$. The primary idea here is that the residual $R_N(x)$ becomes smaller as $\Delta x \rightarrow 0$.

and $R_N(x)$ is the residual (the part of the series we left off)

Example 1

$f(x) = e^x$ with $x_0 = 0$

Using this we can find expressions for the relative and absolute error as a function of $x$ assuming $N=2$.

We can also use the package sympy which has the ability to calculate Taylor polynomials built-in!

Lets plot this numerically for a section of $x$.

Example 2

$f(x) = \frac{1}{x} ~~~~~~ x_0 = 1$, approximate with $\hat{f}(x) = T_2(x)$

Symbols and Definitions

Big-O notation: $f(x) = \text{O}(g(x))$ as $x \rightarrow a$ if and only if $|f(x)| \leq M |g(x)|$ as $|x - a| < \delta$ for $M$ and $a$ positive.

In practice we use Big-O notation to say something about how the terms we may have left out of a series might behave. We saw an example earlier of this with the Taylor's series approximations:

Example:

$f(x) = \sin x$ with $x_0 = 0$ then

We can actually write $f(x)$ then as

This becomes more useful when we look at this as we did before with $\Delta x$:

We can also develop rules for error propagation based on Big-O notation:

In general, there are two theorems that do not need proof and hold when the value of x is large:

Let $$f(x) = p(x) + O(x^n)$$ $$g(x) = q(x) + O(x^m)$$ $$k = \max(n, m)$$ then $$f+g = p + q + O(x^k)$$ $$f \cdot g = p \cdot q + O(x^{n\cdot m})$$

On the other hand, if we are interested in small values of x, say Δx, the above expressions can be modified as follows:

then

Note 1: In this case we suppose that at least the polynomial with $k = \max(n, m)$ has the following form: $$p(\Delta x) = 1 + p_1x + p_2x^2 + ...$$ or $$q(\Delta x) = 1 + q_1x + q_2x^2 + ...$$ so that there is an O(1) term that guarantees the existence of $O(\Delta x^r)$ in the final product.

To get a sense of why we care most about the power on $\Delta x$ when considering convergence the following figure shows how different powers on the convergence rate can effect how quickly we converge to our solution. Note that here we are plotting the same data two different ways. Plotting the error as a function of $\Delta x$ is a common way to show that a numerical method is doing what we expect and exhbits the correct convergence behavior. Since errors can get small quickly it is very common to plot these sorts of plots on a log-log scale to easily visualize the results. Note that if a method was truly of the order $n$ that they will be a linear function in log-log space with slope $n$.

Horner's Method for Evaluating Polynomials

Given

want to find best way to evaluate $P_N(x)$. First consider two ways to write $P_3$:

and using nested multiplication:

Consider how many operations it takes for each...

Adding up all the operations we can in general think of this as a pyramid

We can estimate this way that the algorithm written this way will take approximately $O(N^2 / 2)$ operations to complete.

Looking at our other means of evaluation:

Here we find that the method is $O(N)$ (the 2 is usually ignored in these cases). The important thing is that the first evaluation is $O(N^2)$ and the second $O(N)$!

Algorithm

Fill in the function and implement Horner's method:

def eval_poly(p, x):
    """Evaluates polynomial given coefficients p at x

    Function to evaluate a polynomial in order N operations.  The polynomial is defined as

    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]

    The value x should be a float.
    """
    pass

or an alternative version that allows x to be a vector of values:

def eval_poly(p, x):
    """Evaluates polynomial given coefficients p at x

    Function to evaluate a polynomial in order N operations.  The polynomial is defined as

    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]

    The value x can by a NumPy ndarray.
    """

    y = numpy.ones(x.shape) * p[0]
    for coefficient in p[1:]:
        y = y * x + coefficient

    return y

This version calculates each y value simultaneously making for much faster code!

Truncation Error vs. Floating Point Error

Truncation error: Errors arising from approximation of a function, truncation of a series...

Floating-point Error: Errors arising from approximating real numbers with finite-precision numbers

or $\frac{1}{3} \approx 0.333333333$ in decimal, results form finitely number of registers to represent each number.

Floating Point Systems

Numbers in floating point systems are represented as a series of bits that represent different pieces of a number. In nomalized floating point systems there are some standar conventions for what these bits are used for. In general the numbers are stored by breaking them down into the form

where

1. $\pm$ is a single bit and of course represents the sign of the number
2. $d_1 . d_2 d_3 d_4 \ldots d_p$ is called the mantissa. Note that technically the decimal could be moved but generally, using scientific notation, the decimal can always be placed at this location. The digits $d_2 d_3 d_4 \ldots d_p$ are called the fraction with $p$ digits of precision. Normalized systems specifically put the decimal point in the front like we have and assume $d_1 \neq 0$ unless the number is exactly $0$.
3. $\beta$ is the base. For binary $\beta = 2$, for decimal $\beta = 10$, etc.
4. $E$ is the exponent, an integer in the range $[E_{\min}, E_{\max}]$

The important points on any floating point system is that

1. There exist a discrete and finite set of representable numbers
2. These representable numbers are not evenly distributed on the real line
3. Airthmetic in floating point systems yield different results from infinite precision arithmetic (i.e. "real" math)

Example: Toy System

Consider the toy 2-digit precision decimal system (normalized) $$f = \pm d_1 . d_2 \times 10^E$$ with $E \in [-2, 0]$.

Number and distribution of numbers

1. How many numbers can we represent with this system?

2. What is the distribution on the real line?

3. What is the underflow and overflow limits?

How many numbers can we represent with this system?

$f = \pm d_1 . d_2 \times 10^E$ with $E \in [-2, 0]$.

What is the distribution on the real line?

What is the underflow and overflow limits?

Smallest number that can be represented is the underflow: $1.0 \times 10^{-2} = 0.01$ Largest number that can be represented is the overflow: $9.9 \times 10^0 = 9.9$

Properties of Floating Point Systems

All floating-point systems are characterized by several important numbers

• Smalled normalized number (underflow if below)
• Largest normalized number (overflow if above)
• Zero
• Machine $\epsilon$ or $\epsilon_{\text{machine}}$
• inf and nan, infinity and Not a Number respectively
• Subnormal numbers

Binary Systems

Consider the 2-digit precision base 2 system:

Number and distribution of numbers

1. How many numbers can we represent with this system?

2. What is the distribution on the real line?

3. What is the underflow and overflow limits?

How many numbers can we represent with this system?

What is the distribution on the real line?

Smallest number that can be represented is the underflow: $1.0 \times 2^{-1} = 0.25$ Largest number that can be represented is the overflow: $1.1 \times 2^1 = 2.2$

Real Systems - IEEE 754 Binary Floating Point Systems

Single Precision

• Total storage alloted is 32 bits
• Exponent is 8 bits $\Rightarrow E \in [-126, 127]$
• Fraction 23 bits ($p = 24$)

s EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF
0 1      8 9                     31

Overflow $= 2^{127} \approx 3.4 \times 10^{38}$

Underflow $= 2^{-126} \approx 1.2 \times 10^{-38}$

$\epsilon_{\text{machine}} = 2^{-23} \approx 1.2 \times 10^7$

Double Precision

• Total storage alloted is 64 bits
• Exponent is 11 bits $\Rightarrow E \in [-1022, 1024]$
• Fraction 52 bits ($p = 53$)

s EEEEEEEEEE FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FF
0 1       11 12                                                      63

Overflow $= 2^{1024} \approx 1.8 \times 10^{308}$

Underflow $= 2^{-1022} \approx 2.2 \times 10^{-308}$

$\epsilon_{\text{machine}} = 2^{-52} \approx 2.2 \times 10^{-16}$

Python Access to IEEE Numbers

Access many important parameters, such as machine epsilon:

import numpy
numpy.finfo(float).eps

Why should we care about this?

• Floating point arithmetic is not commutative or associative
• Floating point errors compound, do not assume even double precision is enough!
• Mixing precision is very dangerous

Examples

Simple arithmetic $\delta < \epsilon_{\text{machine}}$

$$(1+\delta) - 1 = 1 - 1 = 0$$

$$1 - 1 + \delta = \delta$$

Catastrophic cancelation (assume $x + y \neq 0$):

This implies that the error made from doing the actual sum can be arbitrarily large! This is not too bad unless $x$ and $y$ have opposite sign, then if $x+y ~ 0$ then $\frac{1}{x+y}$ can be huge. This is refered to as catastrophic cancelation.

Matrix vector multiplication

$$ \begin{bmatrix} \epsilon & 1 \\ 1 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} \epsilon + 1 \\ 1 + 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}$$

Solve $A x = b$ where

$$A = \begin{bmatrix} 10 & 1 \\ 3 & 0.3 \end{bmatrix}~~~~~~\text{and}~~~~~~b=\begin{bmatrix} 11 \\ 3.3 \end{bmatrix}$$

and we know that

$$x = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$$

When used in a solver we get

$$x = \begin{bmatrix}-0.5 \\ 16 \end{bmatrix}$$

Evaluation of a polynomial

$$f(x) = x^7 - 7x^6 + 21 x^5 - 35 x^4 + 35x^3-21x^2 + 7x - 1$$

For more examples see the "100-digit Challenge" book.

Combination of Error

In general we need to concern ourselves with the combination of both truncation error and floating point error.

Example

Consider the finite difference approximation where $f(x) = e^x$ and we are evaluating at $x=1$?

Compare the error between decreasing $\Delta x$ and the true solution $f'(1) = e$.

Example 2

Evaluate $e^x$ with its taylor series.

Can we pick $N < \infty$ that can approximate $e^x$ over a give range $x \in [a,b]$ such that the relative error $E$ satisfies $E < 8 \cdot \varepsilon_{\text{machine}}$?

What might be a better way than simply evaluating the Taylor polynomial directly for various $N$?

Some other links that might be helpful regarding IEEE Floating Point:

• What Every Computer Scientist Should Know About Floating-Point Arithmetic
• IEEE 754 Floating Point Calculator
• Numerical Computing with IEEE Floating Point Arithmetic