Unavoidable Errors in Computing

Gerald W. Recktenwald
Department of Mechanical Engineering
Portland State University
gerry@[Link]
These slides are a supplement to the book Numerical Methods with Matlab:
Implementations and Applications, by Gerald W. Recktenwald, c 2000–2006, Prentice-Hall,
c 2000–2006 Gerald W. Recktenwald.
Upper Saddle River, NJ. These slides are copyright
The PDF version of these slides may be downloaded or stored or printed for noncommercial,
educational use. The repackaging or sale of these slides in any form, without written
consent of the author, is prohibited.

The latest version of this PDF file, along with other supplemental material for the book,
can be found at [Link]/recktenwald or [Link]/~gerry/nmm/.

Version 1.0 August 19, 2006

 Overview (1)

• Digital representation of numbers

⊲ Size limits
⊲ Resolution limits
⊲ The floating point number line

• Floating point arithmetic

⊲ roundoff
⊲ machine precision

NMM: Unavoidable Errors in Computing page 2

 Overview (2)

• Implications for routine computation

⊲ Use “close enough” instead of “equals”
⊲ loss of significance for addition
⊲ catastrophic cancellation for subtraction

• Truncation error
⊲ Demonstrate with Taylor series
⊲ Order Notation

NMM: Unavoidable Errors in Computing page 3

 What’s going on here?

Spontaneous generation of an insignificant digit:

>> format long e % display lots of digits
>> 2.6 + 0.2
ans =
2.800000000000000e+00
>> ans + 0.2
ans =
3.000000000000000e+00
>> ans + 0.2
ans =
3.200000000000001e+00 Why does the least significant digit appear?
>> 2.6 + 0.6
ans =
3.200000000000000e+00 Why does the small error not show up here?

NMM: Unavoidable Errors in Computing page 4

 Bits, Bytes, and Words

base 10 conversion base 2

1 1 = 20 0000 0001
2 2 = 21 0000 0010
4 4 = 22 0000 0100
8 8 = 23 0000 1000
9 8 + 1 = 23 + 20 0000 1001
10 8 + 2 = 23 + 21 0000 1010
27 16 + 8 + 2 + 1 = 24 + 23 + 21 + 20 0001
| {z1011}
one byte

NMM: Unavoidable Errors in Computing page 5

 Digital Storage of Integers (1)

As a prelude to discussing the binary storage of floating point values, first

consider the binary storage of integers.

• Integers can be exactly represented by base 2

• Typical size is 16 bits

• 216 = 65536 is largest 16 bit integer

• [−32768, 32767] is range of 16 bit integers in twos complement notation

• 32 bit and larger integers are available

NMM: Unavoidable Errors in Computing page 6

 Digital Storage of Integers (2)

Note: Unless explicitly specified otherwise, all mathematical calculations

in Matlab use double precision floating point numbers.

Expert’s Note: The built-in int8, int16, int32, uint8, uint16, and
uint32 classes are used to reduce memory usage for
very large data sets.

NMM: Unavoidable Errors in Computing page 7

 Digital Storage of Integers (3)

Let b be a binary digit, i.e. 1 or 0

(bbbb)2 ⇐⇒ |23|22|21|20|

The rightmost bit is the least significant bit (LSB)

The leftmost bit is the most significant bit (MSB)

Example:

(1001)2 = 1 × 23 + 0 × 22 + 0 × 21 + 1 × 20

=8+0+0+1=9

NMM: Unavoidable Errors in Computing page 8

 Digital Storage of Integers (4)

Limitations:

• A finite number of bits are used to store each value in computer

memory.

• Limiting the number of bits limits the size of integer that can be
represented

largest 3 bit integer: (111)2 =4+2+1=7 = 23 − 1

largest 4 bit integer: (1111)2 = 8 + 4 + 2 + 1 = 15 = 24 − 1
largest 5 bit integer: (11111)2 = 16 + 8 + 4 + 2 + 1 = 31 = 25 − 1
largest n bit integer: = 2n − 1

NMM: Unavoidable Errors in Computing page 9

 Digital Storage of Floating Point Numbers (1)

Numeric values with non-zero fractional parts are stored as floating point
numbers.

All floating point values are represented with a normalized scientific

notation1.

Example:

12.2792 = 0.123792 102

Mantissa Exponent

1
The IEEE Standard on Floating Point arithmetic defines a normalized binary format. Here we use a
simplified decimal (base ten) format that, while abusing the standard notation, expresses the essential ideas
behind the decimal to binary conversion.

NMM: Unavoidable Errors in Computing page 10

 Digital Storage of Floating Point Numbers (2)

Floating point values have a fixed number of bits allocated for storage of
the mantissa and a fixed number of bits allocated for storage of the
exponent.

Two common precisions are provided in numeric computing languages

Bits for Bits for
Precision mantissa exponent
Single 23 8
Double 53 11

NMM: Unavoidable Errors in Computing page 11

 Digital Storage of Floating Point Numbers (3)

A double precision (64 bit) floating point number can be schematically

represented as
64 bits
z }| {
b bb . . . . . . bbb bbbbbbbbbbb
|{z} | {z } | {z }
sign 52 bit value 11 bit exponent,
of mantissa including sign

The finite number of bits in the exponent limits the magnitude or range of
the floating point numbers.

The finite number of bits in the mantissa limits the number of significant
digits or the precision of the floating point numbers.

NMM: Unavoidable Errors in Computing page 12

 Digital Storage of Floating Point Numbers (4)

1
The floating point mantissa is expressed in powers of
2
 0
1
= 1 is not used
2
 1  2  3
1 1 1
= 0.5 = 0.25 = 0.125 ···
2 2 2

NMM: Unavoidable Errors in Computing page 13

 Digital Storage of Floating Point Numbers (5)

Algorithm 5.1 Convert Floating-Point to Binary

r0 = x
for k = 1, 2, . . . , m
if rk−1 ≥ 2−k
bk = 1
rk = rk−1 − 2−k
else
bk = 0
rk = rk−1
end if
end for

NMM: Unavoidable Errors in Computing page 14

 Digital Storage of Floating Point Numbers (6)

Example: Binary mantissa for x = 0.8125 — Apply Algorithm 5.1

k 2−k bk rk = rk−1 − bk 2−k

0 NA NA 0.8125
1 0.5 1 0.3125
2 0.25 1 0.0625
3 0.125 0 0.0625
4 0.0625 1 0.0000

Therefore, the binary mantissa for 0.8125 is (exactly) (1101)2

NMM: Unavoidable Errors in Computing page 15

 Digital Storage of Floating Point Numbers (7)

Example: Binary mantissa for x = 0.1 — Apply Algorithm 5.1

k 2−k bk rk = rk−1 − bk 2−k

0 NA NA 0.1
1 0.5 0 0.1
2 0.25 0 0.1
3 0.125 0 0.1
4 0.0625 1 0.1 - 0.0625 = 0.0375
5 0.03125 1 0.0375 - 0.03125 = 0.00625
6 0.015625 0 0.00625
7 0.0078125 0 0.00625
8 0.00390625 1 0.00625 - 0.00390625 = 0.00234375
9 0.001953125 1 0.0234375 - 0.001953125 = 0.000390625
10 0.0009765625 0 0.000390625
.. ..
. .

NMM: Unavoidable Errors in Computing page 16

 Digital Storage of Floating Point Numbers (8)

Calculations on the preceding slide show that

the binary mantissa for 0.1 is (00011 0011 . . .)2.

The decimal value of 0.1 cannot be represented by a finite

number of binary digits.

NMM: Unavoidable Errors in Computing page 17

 Consequences of Finite Storage (1)

Limiting the number of bits =⇒ Upper and lower limits on

allocated for storage of the the range (or magnitude) of
exponent floating point numbers

Limiting the number of bits =⇒ Limit on the precision (or

allocated for storage of the number of significant digits)
mantissa for any floating point number.

NMM: Unavoidable Errors in Computing page 18

 Consequences of Finite Storage (2)

Most real numbers cannot be stored exactly (they do not exist on the
floating point number line)

• Integers less than 252 can be stored exactly.

Try this:
>> x = 2^51
>> s = dec2bin(x)
>> x2 = bin2dec(s)
>> x2 - x

• Numbers with 15 (decimal) digit mantissas that are the exact sum of
powers of (1/2) can be stored exactly.

NMM: Unavoidable Errors in Computing page 19

 Floating Point Number Line
Compare floating point numbers to real numbers.

Real numbers Floating point numbers

Range Infinite: arbitrarily large and Finite: the number of bits
arbitrarily small real numbers allocated to the exponent limit
exist. the magnitude of floating point
values.

Precision Infinite: There is an infinite set Finite: there is a finite number

of real numbers between any (perhaps zero) of floating point
two real numbers. values between any two floating
point values.

In other words: The floating point number line is a subset of the real
number line.

NMM: Unavoidable Errors in Computing page 20

 Floating Point Number Line

denormal

under- under-
flow flow
overflow usable range usable range overflow

–10+308 –10-308 0 10-308 10+308

–realmax –realmin realmin realmax

zoom-in view

NMM: Unavoidable Errors in Computing page 21

 Symbolic versus Numeric Calculation (1)

Commercial software for symbolic computation

• DeriveTM

• MACSYMATM

• MapleTM

• MathematicaTM

Symbolic calculations are exact. No rounding occurs because symbols and

algebraic relationships are manipulated without storing numerical values.

NMM: Unavoidable Errors in Computing page 22

 Symbolic versus Numeric Calculation (2)

Example: Evaluate f (θ) = 1 − sin2 θ − cos2 θ

>> theta = 30*pi/180;

>> f = 1 - sin(theta)^2 - cos(theta)^2
f =
-1.1102e-16

f is close to, but not exactly equal to zero because of roundoff. Also note
that f is a single value, not a formula.

NMM: Unavoidable Errors in Computing page 23

 Symbolic versus Numeric Calculation (3)

Symbolic computation using the Symbolic Math Toolbox in Matlab

>> t = sym(’t’) % declare t as a symbolic variable
t =
t

>> f = 1 - sin(t)^2 - cos(t)^2 % create a symbolic expression

f =
1-sin(t)^2-cos(t)^2

>> simplify(f) % ask Maple engine to make algebraic simplifications

f =
0

In the symbolic computation, f is exactly zero for any value of t. There is

no roundoff error in symbolic computation.

NMM: Unavoidable Errors in Computing page 24

 Numerical Arithmetic

Numerical values have limited range and precision. Values created by

adding, subtracting, multiplying, or dividing floating point values will also
have limited range and precision.

Quite often, the result of an arithmetic operation between two floating

point values cannot be represented as another floating point value.

NMM: Unavoidable Errors in Computing page 25

 Integer Arithmetic

Operation Result
2+2=4 integer
9 × 7 = 63 integer
12
=4 integer
3
29
=2 exact result is not an integer
13
29
=0 exact result is not an integer
1300

NMM: Unavoidable Errors in Computing page 26

 Floating Point Arithmetic

Operation Floating Point Value is . . .

2.0 + 2.0 = 4 exact
9.0 × 7.0 = 63 exact
12.0
=4 exact
3.0
29
= 2.230769230769231 approximate
13
29
= 2.230769230769231 × 10−2 approximate
1300

NMM: Unavoidable Errors in Computing page 27

 Floating Point Arithmetic in Matlab (1)

>> format long e

>> u = 29/13
u =
2.230769230769231e+00

>> v = 13*u
v =
29

>> v-29
ans =
0

Two rounding errors are made: (1) during computation and storage of u,
and (2) during computation and storage of v. Fortuitously, the
combination of rounding errors produces the exact result.

NMM: Unavoidable Errors in Computing page 28

 Floating Point Arithmetic in Matlab (2)

>> x = 29/1300
x =
2.230769230769231e-02
>> y = 29 - 1300*x
y =
3.552713678800501e-015

In exact arithmetic, the value of y should be zero.

The roundoff error occurs when x is stored. Since 29/1300 cannot be

expressed with a finite sum of the powers of 1/2, the numerical value
stored in x is a truncated approximation to 29/1300.

When y is computed, the expression 1300*x evaluates to a number slightly

different than 29 because the bits lost in the computation and storage of x
are not recoverable.

NMM: Unavoidable Errors in Computing page 29

 Roundoff in Quadratic Equation (1)
(See Example 5.3 in the text)
The roots of
ax2 + bx + c = 0 (1)
are √
−b ± b2 − 4ac
x= (2)
2a
Consider
x2 + 54.32x + 0.1 = 0 (3)
which has the roots (to eleven digits)

x1 = 54.3218158995, x2 = 0.0018410049576.

Note that b2 ≫ 4ac

b2 = 2950.7 ≫ 4ac = 0.4

NMM: Unavoidable Errors in Computing page 30

 Roundoff in Quadratic Equation (2)

Compute roots with four digit arithmetic

p q
2
b2 − 4ac = (−54.32) − 0.4000
√
= 2951 − 0.4000
√
= 2951
= 54.32

The result of each intermediate mathematical operation is rounded to four

digits.

NMM: Unavoidable Errors in Computing page 31

 Roundoff in Quadratic Equation (3)

Use x1,4 to designate the first root computed with four-digit arithmetic:
√
−b + b2 − 4ac
x1,4 =
2a
+54.32 + 54.32
=
2.000
108.6
=
2.000
= 54.30

Correct root is x1 = 54.3218158995. Four digit arithmetic leads to 0.4

percent error in this example.

NMM: Unavoidable Errors in Computing page 32

 Roundoff in Quadratic Equation (4)

Using four-digit arithmetic the second root, x2,4, is

√
−b − b2 − 4ac
x2,4 =
2a
+54.32 − 54.32
= (i)
2.000
0.000
= (ii)
2.000
=0 (iii)

An error of 100 percent!

The poor approximation to x2,4 is caused by roundoff in the calculation of
√
b2 − 4ac. This leads to the subtraction of two equal numbers in line (i).

NMM: Unavoidable Errors in Computing page 33

 Roundoff in Quadratic Equation (5)

A solution: rationalize the numerators of the expressions for the two roots:
√ √ !
−b + b2 − 4ac −b − b2 − 4ac 2c
x1 = √ = √ , (4)
2a 2
−b − b − 4ac 2
−b − b − 4ac

√ √ !
−b − b2 − 4ac −b + b2 − 4ac 2c
x2 = √ = √ (5)
2a 2
−b + b − 4ac −b + b2 − 4ac

NMM: Unavoidable Errors in Computing page 34

 Roundoff in Quadratic Equation (6)

Now use Equation (5) to compute the troublesome second root with four
digit arithmetic

2c 0.2000 0.2000
x2,4 = √ = = = 0.001842.
2
−b + b − 4ac +54.32 + 54.32 108.6

The result is in error by only 0.05 percent.

NMM: Unavoidable Errors in Computing page 35

 Roundoff in Quadratic Equation (7)

Compare the formulas for x2

√
−b − b2 − 4ac
x2,std =
2a
2c
x2,new = √
−b + b2 − 4ac

The two formulations for x2 are algebraically equivalent. The difference in

the computed values of x2,4 is due to roundoff alone.

NMM: Unavoidable Errors in Computing page 36

 Roundoff in Quadratic Equation (8)

Repeat the calculation of x1,4 with the new formula

2c
x1,4 = √
−b − b2 − 4ac
0.2000
= (i)
+54.32 − 54.32
0.2000
= (ii)
0
= ∞.
√
Limited precision in the calculation of b2 + 4ac leads to a catastrophic
cancellation error in step (i)

NMM: Unavoidable Errors in Computing page 37

 Roundoff in Quadratic Equation (9)

A robust solution uses a formula that takes the sign of b into account in a
way that prevents catastrophic cancellation.
The ultimate quadratic formula:
1h p i
q ≡ − b + sign(b) b2 − 4ac
2
where (
1 if b ≥ 0,
sign(b) =
−1 otherwise
Then roots to quadratic equation are
q c
x1 = x2 =
a q

NMM: Unavoidable Errors in Computing page 38

 Roundoff in Quadratic Equation (10)

Summary

• Finite-precision causes roundoff in individual calculations

• Effects of roundoff usually accumulate slowly, but . . .

• Subtracting nearly equal numbers leads to severe loss of precision. A

similar loss of precision occurs when two numbers of very different
magnitude are added.

• Roundoff is inevitable: good algorithms minimize the effect of roundoff.

NMM: Unavoidable Errors in Computing page 39

 Catastrophic Cancellation Errors (1)

The errors in
c=a+b and c=a−b
will be large when a ≫ b or a ≪ b.

Consider c = a + b with

a = [Link] . . . × 100
b = [Link] . . . × 10−8

where x and y are decimal digits.

NMM: Unavoidable Errors in Computing page 40

 Catastrophic Cancellation Errors (1)

Evaluate c = a + b with a = [Link] . . . × 100 and b = [Link] . . . × 10−8

Assume for convenience of exposition that z = x + y < 10.

available precision
z }| {
[Link] xxxx xxxx xxxx
+ 0.000 0000 yyyy yyyy yyyy yyyy
= [Link] xxxx zzzz zzzz yyyy yyyy
| {z }
lost digits

The most significant digits of a are retained, but the least significant digits
of b are lost because of the mismatch in magnitude of a and b.

NMM: Unavoidable Errors in Computing page 41

 Catastrophic Cancellation Errors (2)

For subtraction: The error in

c=a−b

will be large when a ≈ b.

Consider c = a − b with

a = x.xxxxxxxxxxx1ssssss
b = x.xxxxxxxxxxx0tttttt

where x, y, s and t are decimal digits. The digits sss . . . and ttt . . . are
lost when a and b are stored in double-precision, floating point format.

NMM: Unavoidable Errors in Computing page 42

 Catastrophic Cancellation Errors (3)

Evaluate a − b in double precision floating point arithmetic when

a = [Link] xxxx xxxx 1 and b = [Link] xxxx xxxx 0
available precision
z }| {
[Link] xxxx xxxx 1
− [Link] xxxx xxxx 0
= 0.000 0000 0000 1 uuuu uuuu uuuu
| {z }
unassigned digits
−12
= [Link] uuuu uuuu × 10
The result has only one significant digit. Values for the uuuu digits are not
necessarily zero. The absolute error in the result is small compared to
either a or b. The relative error in the result is large because
ssssss − tttttt 6= uuuuuu (except by chance).

NMM: Unavoidable Errors in Computing page 43

 Catastrophic Cancellation Errors (4)

Summary

• Occurs in addition α + β or subtraction α − β when α ≫ β or α ≪ β

• Occurs in subtraction: α − β when α ≈ β

• Error caused by a single operation (hence the term “catastrophic”) not

a slow accumulation of errors.

• Can often be minimized by algebraic rearrangement of the troublesome

formula. (Cf. improved quadratic formula.)

NMM: Unavoidable Errors in Computing page 44

 Machine Precision (1)

The magnitude of roundoff errors is quantified by machine precision εm.

There is a number, εm > 0, such that

1+δ =1

whenever |δ| < εm.

In exact arithmetic 1 + δ = 1 only when δ = 0, so in exact arithmetic εm is
identically zero.
Matlab uses double precision (64 bit) arithmetic. The built-in variable
eps stores the value of εm.

eps = 2.2204 × 10−16

NMM: Unavoidable Errors in Computing page 45

 Machine Precision (2)

Matlab code for Computing Machine Precision:

epsilon = 1;
it = 0;
maxit = 100;
while it < maxit
epsilon = epsilon/2;
b = 1 + epsilon;
if b == 1
break;
end
it = it + 1;
end

NMM: Unavoidable Errors in Computing page 46

 Implications for Routine Calculations

• Floating point comparisons should test for “close enough” instead of

exact equality.

• Express “close” in terms of

absolute difference, |x − y|
or
|x − y|
relative difference,
|x|

NMM: Unavoidable Errors in Computing page 47

 Floating Point Comparison

Don’t ask “is x equal to y”.

if x==y % Don’t do this

...
end

Instead ask, “are x and y ‘close enough’ in value”

if abs(x-y) < tol

...
end

NMM: Unavoidable Errors in Computing page 48

 Absolute and Relative Error (1)

“Close enough” can be measured with either absolute difference or relative

difference, or both

Let

α = some exact or reference value

α
b = some computed value

Absolute error  
Eabs(b 
α) = αb − α

NMM: Unavoidable Errors in Computing page 49

 Absolute and Relative Error (1)

Relative error  
αb − α
α) = 
Erel(b 
αref 

Often we choose αref = α so that

 
αb − α
α) =  
Erel(b
α

NMM: Unavoidable Errors in Computing page 50

 Absolute and Relative Error (2)

Example: Approximating sin(x) for small x

Since
x3 x5
sin(x) = x − + − ...
3! 5!
we can approximate sin(x) with

sin(x) ≈ x

for small enough |x| < 1

NMM: Unavoidable Errors in Computing page 51

 Absolute and Relative Error (3)

The absolute error in approximating sin(x) ≈ x for small x is

x3 x5
Eabs = x − sin(x) = − + ...
3! 5!

And the relative error is

x − sin(x) x
Eabs = = −1
sin(x) sin(x)

NMM: Unavoidable Errors in Computing page 52

 Absolute and Relative Error (4)

Plot relative and absolute error in approximating sin(x) with x.

−3
x 10 Error in approximating sin(x) with x
Although the absolute error is 20
Absolute Error
Relative Error
relatively flat around x = 0,
the relative error grows more 15

quickly.
10
The relative error grows

Error
quickly because the absolute 5
value of sin(x) is small near
x = 0. 0

−5
−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3
x (radians)

NMM: Unavoidable Errors in Computing page 53

 Iteration termination (1)

An iteration generates a sequence of scalar values xk , k = 1, 2, 3, . . ..

The sequence converges to a limit ξ if

|xk − ξ| < δ, for all k > N,

where δ is a small.

In practice, the test is usually expressed as

|xk+1 − xk | < δ, when k > N.

NMM: Unavoidable Errors in Computing page 54

 Iteration termination (2)

Absolute convergence criterion

Iterate until |x − xold| < ∆a where ∆a is the absolute convergence

tolerance.

In Matlab:
x = ... % initialize
xold = ...

while abs(x-xold) > deltaa

xold = x;
update x
end

NMM: Unavoidable Errors in Computing page 55

 Iteration termination (3)

Relative convergence criterion

 
 x − xold 
Iterate until   < δr , where δr is the absolute convergence
xold 
tolerance.

In Matlab:
x = ... % initialize
xold = ...

while abs((x-xold)/xold) > deltar

xold = x;
update x
end

NMM: Unavoidable Errors in Computing page 56

 Example: Solve cos(x) = x (1)

1.6

Find the value of x that 1.4

satisfies cos(x) = x. 1.2

y = cos(x)
1

0.8

and
0.6

y=x
0.4

0.2

0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
x (radians)

NMM: Unavoidable Errors in Computing page 57

 Example: Solve cos(x) = x (2)

The fixed point iteration as a method for obtaining a numerical

approximation to the solution of a scalar equation. For now, trust that the
follow algorithm will eventually give the solution.

1. Guess x0

2. Set xold = x0

3. Update guess
xnew = cos(xold)

4. If xnew ≈ xold stop; otherwise set xold = xnew and return to step 3

NMM: Unavoidable Errors in Computing page 58

 Solve cos(x) = x (3)

MATLAB implementation

x0 = ... % initial guess

k = 0;
xnew = x0;
while NOT_CONVERGED & k < maxit
xold = xnew;
xnew = cos(xold);
it = it + 1;
end

Let’s examine someways of defining the logical value NOT_CONVERGED.

NMM: Unavoidable Errors in Computing page 59

 Solve cos(x) = x (4)

Bad test # 1

while xnew ~= xold

This test will be true unless xnew and xold are exactly equal. In other
words, xnew and xold are equal only when their bit patterns are identical.
This is bad because

• Test may never be met because of oscillatory bit patterns

• Even if test is eventually met, the iterations will probably do more work
than needed

NMM: Unavoidable Errors in Computing page 60

 Solve cos(x) = x (5)

Bad test # 2

while (xnew-xold) > delta

This test evaluates to false whenever (xnew-xold) is negative,

even if |(xnew − xold)| ≫ delta.
Example:

>> xold = 100; xnew = 1; delta = 5e-9;

>> (xnew-xold) > delta
ans =
0

These values of xnew and xold are not close, but the erroneous
convergence criterion would cause the iterations to stop.

NMM: Unavoidable Errors in Computing page 61

 Solve cos(x) = x (6)

Workable test # 1: Absolute tolerance

while abs(xnew-xold) > delta

An absolute tolerance is useful when the iterative sequence converges to a

value with magnitude much less than one.

What value of delta to use?

NMM: Unavoidable Errors in Computing page 62

 Solve cos(x) = x (7)

Workable test # 2: Relative tolerance

while abs( (xnew-xold)/xref ) > delta

The user supplies appropriate value of xref. For this particular iteration
we could use xref = xold.

while abs( (xnew-xold)/xold ) > delta

Note: For the problem of solving cos(x) = x, the solution is O(1) so

the absolute and relative convergence tolerance will terminate
the calculations at roughly the same iteration.

NMM: Unavoidable Errors in Computing page 63

 Solve cos(x) = x (8)

Using the relative convergence tolerance, the code becomes

x0 = ... % initial guess

k = 0;
xnew = x0;
while ( abs( (xnew-xold)/xold ) > delta ) & k < maxit
xold = xnew;
xnew = cos(xold);
it = it + 1;
end

Note: Parentheses around abs( (xnew-xold)/xold ) > delta are

not needed for correct Matlab implementation. The
parenthesis are added to make the meaning of the clear to
humans reading the code.

NMM: Unavoidable Errors in Computing page 64

 Truncation Error

Consider the series for sin(x)

x3 x5
sin(x) = x − + − ···
3! 5!

For small x, only a few terms are needed to get a good approximation to
sin(x). The · · · terms are “truncated”

ftrue = fsum + truncation error

The size of the truncation error depends on x and the number of terms
included in fsum.

NMM: Unavoidable Errors in Computing page 65

 Truncation of series for sin(x) (1)

function ssum = sinser(x,tol,n)

% sinser Evaluate the series representation of the sine function
%
% Input: x = argument of the sine function, i.e., compute sin(x)
% tol = tolerance on accumulated sum. Default: tol = 5e-9
% n = maximum number of terms. Default: n = 15
%
% Output: ssum = value of series sum after nterms or tolerance is met

term = x; ssum = term; % Initialize series

fprintf(’Series approximation to sin(%f)\n\n k term ssum\n’,x);
fprintf(’%3d %11.3e %12.8f\n’,1,term,ssum);
for k=3:2:(2*n-1)
term = -term * x*x/(k*(k-1)); % Next term in the series
ssum = ssum + term;
fprintf(’%3d %11.3e %12.8f\n’,k,term,ssum);
if abs(term/ssum)<tol, break; end % True at convergence
end
fprintf(’\nTruncation error after %d terms is %g\n\n’,(k+1)/2,abs(ssum-sin(x)));

NMM: Unavoidable Errors in Computing page 66

 Truncation of series for sin(x) (2)

For small x, the series for sin(x) converges in a few terms

>> s = sinser(pi/6,5e-9,15);
Series approximation to sin(0.523599)

k term ssum
1 5.236e-001 0.52359878
3 -2.392e-002 0.49967418
5 3.280e-004 0.50000213
7 -2.141e-006 0.49999999
9 8.151e-009 0.50000000
11 -2.032e-011 0.50000000

Truncation error after 6 terms is 3.56382e-014

NMM: Unavoidable Errors in Computing page 67

 Truncation of series for sin(x) (3)

The truncation error in the series is small relative to the true value of
sin(π/6)

>> s = sinser(pi/6,5e-9,15);
.
.
>> err = (s-sin(pi/6))/sin(pi/6)
err =
-7.1276e-014

NMM: Unavoidable Errors in Computing page 68

 Truncation of series for sin(x) (4)

For larger x, the series for sin(x) converges more slowly

>> s = sinser(15*pi/6,5e-9,15);
Series approximation to sin(7.853982)

k term ssum
1 7.854e+000 7.85398163
3 -8.075e+001 -72.89153055
5 2.490e+002 176.14792646
.
. .
. .
.
. . .
25 1.537e-003 1.00012542
27 -1.350e-004 0.99999038
29 1.026e-005 1.00000064

Truncation error after 15 terms is 6.42624e-007

Increasing the number of terms will allow the series to converge with the
tolerance of 5 × 10−9. A better solution to the slow convergence of the
series are explored in Exercise 23.

NMM: Unavoidable Errors in Computing page 69

 Taylor Series (1)

For a sufficiently continuous function f (x) defined on the interval

x ∈ [a, b] we define the nth order Taylor Series approximation Pn(x)
 2 2 
 

df  (x − x0) d f  n n 
(x − x0) d f 
Pn(x) = f (x0)+(x−x0)  + +· · ·+
dx x=x0 2 dx2 x=x0 n! dxn x=x0

Then there exists ξ with x0 ≤ ξ ≤ x such that

f (x) = Pn(x) + Rn(x)

where 
(n+1) (n+1)
(x − x0) d f 
Rn(x) =
(n + 1)! dx(n+1)  x=ξ

NMM: Unavoidable Errors in Computing page 70

 Taylor Series (2)

Big “O” notation

 (n+1)

(x − x0)
f (x) = Pn(x) + O
(n + 1)!

or, for x − x0 = h we say

 
(n+1)
f (x) = Pn(x) + O h

NMM: Unavoidable Errors in Computing page 71

 Taylor Series Example

Consider the function

1
f (x) =
1−x
The Taylor Series approximations to f (x) of order 1, 2 and 3 are

1 x − x0
P1(x) = +
1 − x0 (1 − x0)2

1 x − x0 (x − x0)2
P2(x) = + +
1 − x0 (1 − x0)2 (1 − x0)3

1 x − x0 (x − x0)2 (x − x0)3
P3(x) = + + +
1 − x0 (1 − x0)2 (1 − x0)3 (1 − x0)4

NMM: Unavoidable Errors in Computing page 72

 Taylor Series: Pi(x) near x = 1.6 (3)

−0.5

−1

−1.5
Approximations to f(x) = 1/(1−x)

−2

−2.5

−3

−3.5

−4 exact
P1(x)
P (x)
−4.5 2
P (x)
3
−5
1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
x

NMM: Unavoidable Errors in Computing page 73

 Roundoff and Truncation Errors (1)

Roundoff and truncation errors occur in numerical computation.

Example:
Finite difference approximation to f ′(x) = df /dx

′ f (x + h) − f (x) h ′′
f (x) = − f (x) + . . .
h 2
This approximation is said to be first order because the leading term in the
truncation error is linear in h. Dropping the truncation error terms we
obtain

′ f (x + h) − f (x) ′
ffd (x) = or ffd (x) = f ′(x) + O(h)
h

NMM: Unavoidable Errors in Computing page 74

 Roundoff and Truncation Errors (2)

To study the roles of roundoff and truncation errors, compute the finite
difference2 approximation to f ′(x) when f (x) = ex.

′ d x
The relative error in the ffd (x) approximation to e is
dx
′
ffd (x) − f ′(x) ffd
′
(x) − ex
Erel = =
f ′(x) ex

2
The finite difference approximation is used to obtain numerical solutions to ordinary and partial
differentials equations where f (x) and hence f ′ (x) is unknown.

NMM: Unavoidable Errors in Computing page 75

 Roundoff and Truncation Errors (3)

Evaluate Erel at x = 1 for a range of h.

0
Truncation error 10

dominates at large h. 10
−1

Roundoff error Truncation error

−2
10 dominates dominates
Roundoff error in
−3
f (x + h) − f (h) 10

dominates as h → 0. Relative error 10

−4

−5
10

−6
10

−7
10

−8
10 −12 −10 −8 −6 −4 −2 0
10 10 10 10 10 10 10
Stepsize, h

NMM: Unavoidable Errors in Computing page 76