Mathematical preliminaries and error analysis

(1)

師大

analysis

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

September 21, 2014

1 / 116

(2)

師大

Outline

1 Round-off errors and computer arithmetic IEEE standard floating-point format

Absolute and Relative Errors Machine Epsilon

Loss of Significance

2 Algorithms and Convergence Algorithm

Stability

Rate of convergence

2 / 116

(3)

師大

Outline

1 Round-off errors and computer arithmetic IEEE standard floating-point format

Absolute and Relative Errors Machine Epsilon

Loss of Significance

2 Algorithms and Convergence Algorithm

Stability

Rate of convergence

3 / 116

(4)

師大

What is the difference for the arithmetic in algebra and computer?

1 For arithmetic in algebra, 256 + 1 = 257, √

256 + 12

= 257

2 For arithmetic in computer (MATLAB), int8(256) +int8(1) =127 ???????

int16(256) +int16(1) = 257

sqrt(256+1)ˆ2 = ? The solution is equal to 257 or not.

(single(sqrt(5))+single(sqrt(3)))ˆ2 - (sqrt(3)+sqrt(5))ˆ2

4 / 116

(5)

師大

Example 1

Consider the following recurrence algorithm

x₀= 1, x₁ = ¹₃ xn+1= ¹³₃ xn−⁴₃xn−1

for computing the sequence of {x_n= (¹₃)ⁿ}.

Matlab program

n = 30; x = zeros(n,1); x(1) = 1; x(2) = 1/3;

for ii = 3:n

x(ii) = 13 / 3 * x(ii-1) - 4 / 3 * x(ii-2);

xn = (1/3)ˆ(ii-1); RelErr = abs(xn-x(ii)) / xn;

fprintf(’x(%2.0f) = %15.8e, x ast(%2.0f) = %14.8e,’, ...

’RelErr(%2.0f) = %11.4e \n’, ii,x(ii),ii,xn,ii,RelErr);

end

5 / 116

(6)

師大

Example 2

What is the binary representation of ²₃?

Solution: To determine the binary representation for ²₃, we write 2

3 = (0.a₁a₂a₃. . .)₂. Multiply by 2to obtain

4

3 = (a₁.a₂a₃. . .)₂.

Therefore, we get a₁ = 1by taking the integer part of both sides.

6 / 116

(7)

師大

Example 2

4

3 = (a₁.a₂a₃. . .)₂.

7 / 116

(8)

師大

Example 2

4

3 = (a₁.a₂a₃. . .)₂.

8 / 116

(9)

師大

Subtracting 1, we have 1

3 = (0.a₂a₃a₄. . .)₂. Repeating the previous step, we arrive at

2

3 = (0.101010 . . .)₂.

9 / 116

(10)

師大

Subtracting 1, we have 1

3 = (0.a₂a₃a₄. . .)₂. Repeating the previous step, we arrive at

2

3 = (0.101010 . . .)₂.

10 / 116

(11)

師大

In the computational world, each representable number has only afixedandfinitenumber of digits.

For any real number x, let

x = ±1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m,

denote the normalized scientific binary representation of x.

In 1985, the IEEE (Institute for Electrical and Electronic Engineers) published a report called Binary Floating Point Arithmetic Standard 754-1985. In this report, formats were specified for single, double, and extended precisions, and these standards are generally followed by microcomputer manufactures using floating-point hardware.

11 / 116

(12)

師大

x = ±1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m,

12 / 116

(13)

師大

x = ±1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m,

13 / 116

(14)

師大

Single precision

The single precision IEEE standard floating-point format allocates 32 bits for the normalized floating-point number

±q × 2^m as shown in the following figure.

23 bits sign of mantissa

normalized mantissa exponent

8 bits

0 1 8 9 31

Thefirst bitis asignindicator, denoted s. This is followed by an8-bit exponent cand a23-bit mantissa f.

The base for the exponent and mantissa is 2, and the actualexponent isc − 127. The value of c is restricted by the inequality 0 ≤ c ≤ 255.

14 / 116

(15)

師大

Single precision

8 bits

0 1 8 9 31

15 / 116

(16)

師大

Single precision

8 bits

0 1 8 9 31

16 / 116

(17)

師大

The actual exponent of the number is restricted by the inequality−127 ≤ c − 127 ≤ 128.

A normalization is imposed that requires that the leading digit in fraction be 1, and this digit is not stored as part of the 23-bit mantissa.

Using this system gives a floating-point number of the form (−1)^s2^c−127(1 + f ).

17 / 116

(18)

師大

18 / 116

(19)

師大

19 / 116

(20)

師大

Example 3

What is the decimal number of the machine number 01000000101000000000000000000000?

1 The leftmost bit is zero, which indicates that the number is positive.

2 The next 8 bits, 10000001, are equivalent to c = 1 · 2⁷+ 0 · 2⁶+ · · · + 0 · 2¹+ 1 · 2⁰= 129.

The exponential part of the number is 2^129−127= 2².

3 The final 23 bits specify that the mantissa is

f = 0 · (2)⁻¹+ 1 · (2)⁻²+ 0 · (2)⁻³+ · · · + 0 · (2)⁻²³= 0.25.

4 Consequently, this machine number precisely represents the decimal number

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.25) = 5.

20 / 116

(21)

師大

Example 3

f = 0 · (2)⁻¹+ 1 · (2)⁻²+ 0 · (2)⁻³+ · · · + 0 · (2)⁻²³= 0.25.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.25) = 5.

21 / 116

(22)

師大

Example 3

f = 0 · (2)⁻¹+ 1 · (2)⁻²+ 0 · (2)⁻³+ · · · + 0 · (2)⁻²³= 0.25.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.25) = 5.

22 / 116

(23)

師大

Example 3

f = 0 · (2)⁻¹+ 1 · (2)⁻²+ 0 · (2)⁻³+ · · · + 0 · (2)⁻²³= 0.25.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.25) = 5.

23 / 116

(24)

師大

Example 3

f = 0 · (2)⁻¹+ 1 · (2)⁻²+ 0 · (2)⁻³+ · · · + 0 · (2)⁻²³= 0.25.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.25) = 5.

24 / 116

(25)

師大

Example 4

f = 0 · (2)⁻¹+ 0 · (2)⁻²+ 1 · (2)⁻³+ · · · + 1 · (2)⁻²³

= 0.2499998807907105.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2499998807907105)

= 4.999999523162842.

25 / 116

(26)

師大

Example 4

f = 0 · (2)⁻¹+ 0 · (2)⁻²+ 1 · (2)⁻³+ · · · + 1 · (2)⁻²³

= 0.2499998807907105.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2499998807907105)

= 4.999999523162842.

26 / 116

(27)

師大

Example 4

f = 0 · (2)⁻¹+ 0 · (2)⁻²+ 1 · (2)⁻³+ · · · + 1 · (2)⁻²³

= 0.2499998807907105.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2499998807907105)

= 4.999999523162842.

27 / 116

(28)

師大

Example 5

f = 0 · 2⁻¹+ 1 · 2⁻²+ 0 · 2⁻³+ · · · + 0 · 2⁻²²+ 1 · 2⁻²³

= 0.2500001192092896.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2500001192092896)

= 5.000000476837158.

28 / 116

(29)

師大

Example 5

f = 0 · 2⁻¹+ 1 · 2⁻²+ 0 · 2⁻³+ · · · + 0 · 2⁻²²+ 1 · 2⁻²³

= 0.2500001192092896.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2500001192092896)

= 5.000000476837158.

29 / 116

(30)

師大

Example 5

f = 0 · 2⁻¹+ 1 · 2⁻²+ 0 · 2⁻³+ · · · + 0 · 2⁻²²+ 1 · 2⁻²³

= 0.2500001192092896.

(−1)^s2^c−127(1 + f ) = 2²· (1 + 0.2500001192092896)

= 5.000000476837158.

30 / 116

(31)

師大

Summary

Above three examples

01000000100111111111111111111111 ⇒ 4.999999523162842 01000000101000000000000000000000 ⇒ 5

01000000101000000000000000000001 ⇒ 5.000000476837158 Only a relativelysmall subsetof the real number system is used for the representation of all the real numbers.

This subset, which are called thefloating-point numbers, contains only rational numbers, both positive and negative.

When a number can not be represented exactly with the fixed finite number of digits in a computer, anear-by floating-point number is chosen for approximate representation.

31 / 116

(32)

師大

Summary

01000000100111111111111111111111 ⇒ 4.999999523162842 01000000101000000000000000000000 ⇒ 5

32 / 116

(33)

師大

Summary

01000000100111111111111111111111 ⇒ 4.999999523162842 01000000101000000000000000000000 ⇒ 5

33 / 116

(34)

師大

Summary

01000000100111111111111111111111 ⇒ 4.999999523162842 01000000101000000000000000000000 ⇒ 5

34 / 116

(35)

師大

The smallest positive number

Let s = 0, c = 1 and f = 0 which is equivalent to 2⁻¹²⁶· (1 + 0) ≈ 1.175 × 10⁻³⁸ The largest number

Let s = 0, c = 254 and f = 1 − 2⁻²³which is equivalent to 2¹²⁷· (2 − 2⁻²³) ≈ 3.403 × 10³⁸

Definition 6

If a number x with |x| < 2⁻¹²⁶· (1 + 0), then we say that an underflowhas occurred and is generally set to zero.

If |x| > 2¹²⁷· (2 − 2⁻²³), then we say that anoverflowhas occurred.

35 / 116

(36)

師大

Definition 6

36 / 116

(37)

師大

Definition 6

37 / 116

(38)

師大

Definition 6

38 / 116

(39)

師大

Double precision

A floating point number in double precision IEEE standard format uses two words (64 bits) to store the number as shown in the following figure.

1 sign of mantissa

52-bit mantissa 0 1

11-bit

11 12

63

Thefirstbit is a sign indicator, denoted s. This is followed by an11-bitexponent c and a52-bitmantissa f .

The actual exponent isc − 1023.

39 / 116

(40)

師大

Double precision

1 sign of mantissa

52-bit mantissa 0 1

11-bit

11 12

63

40 / 116

(41)

師大

Double precision

1 sign of mantissa

52-bit mantissa 0 1

11-bit

11 12

63

41 / 116

(42)

師大

Format of floating-point number

(−1)^s× (1 + f ) × 2^c−1023 The smallest positive number

Let s = 0, c = 1 and f = 0 which is equivalent to 2⁻¹⁰²²· (1 + 0) ≈ 2.225 × 10⁻³⁰⁸. The largest number

Let s = 0, c = 2046 and f = 1 − 2⁻⁵²which is equivalent to 2¹⁰²³· (2 − 2⁻⁵²) ≈ 1.798 × 10³⁰⁸.

42 / 116

(43)

師大

43 / 116

(44)

師大

44 / 116

(45)

師大

Chopping and rounding

x = ±1.a1a2· · · a_tat+1at+2· · · × 2^m,

1 chopping: simply discard the excess bits a_t+1, at+2, . . .to obtain

f l(x) = ±1.a1a2· · · a_t× 2^m.

2 rounding: add 2^−(t+1)× 2^m to x and then chop the excess bits to obtain a number of the form

f l(x) = ±1.δ₁δ₂· · · δ_t× 2^m.

In this method, if at+1= 1, we add 1 to atto obtain f l(x), and if a_t+1 = 0, we merely chop off all but the first t digits.

45 / 116

(46)

師大

Chopping and rounding

x = ±1.a1a2· · · a_tat+1at+2· · · × 2^m,

f l(x) = ±1.a1a2· · · a_t× 2^m.

f l(x) = ±1.δ₁δ₂· · · δ_t× 2^m.

46 / 116

(47)

師大

Chopping and rounding

x = ±1.a1a2· · · a_tat+1at+2· · · × 2^m,

f l(x) = ±1.a1a2· · · a_t× 2^m.

f l(x) = ±1.δ₁δ₂· · · δ_t× 2^m.

47 / 116

(48)

師大

Definition 7 (Roundoff error)

The error results from replacing a number with its floating-point form is calledroundoff errororrounding error.

Definition 8 (Absolute Error and Relative Error)

If x is an approximation to the exact value x^∗, theabsolute error is|x^∗− x|and therelative erroris ^|x_|x^∗^−x|∗| , provided that x^∗ 6= 0.

Example 9

(a) If x^∗ = 0.3000 × 10⁻³and x = 0.3100 × 10⁻³, then the absolute error is 0.1 × 10⁻⁴and the relative error is 0.3333 × 10⁻¹.

(b) If x^∗ = 0.3000 × 10⁴and x = 0.3100 × 10⁴, then the absolute error is 0.1 × 10³ and the relative error is 0.3333 × 10⁻¹.

48 / 116

(49)

師大

Example 9

49 / 116

(50)

師大

Example 9

50 / 116

(51)

師大

Example 9

51 / 116

(52)

師大

Remark 1

As a measure of accuracy, the absolute error may be misleading and the relative error more meaningful.

Definition 10

The number x is said to approximate x^∗to tsignificant digitsif t is the largest nonnegative integer for which

|x − x^∗|

|x^∗| ≤ 5 × 10^−t.

52 / 116

(53)

師大

Remark 1

As a measure of accuracy, the absolute error may be misleading and the relative error more meaningful.

Definition 10

The number x is said to approximate x^∗to tsignificant digitsif t is the largest nonnegative integer for which

|x − x^∗|

|x^∗| ≤ 5 × 10^−t.

53 / 116

(54)

師大

If the floating-point representation f l(x) for the number x is obtained by using t digits and chopping procedure, then the relative error is

|x − f l(x)|

|x| = |0.00 · · · 0a_t+1at+2· · · × 2^m|

|1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m|

= |0.a_t+1at+2· · · |

|1.a₁a2· · · a_tat+1at+2· · · | × 2^−t. The minimal value of the denominator is 1. The numerator is bounded above by 1. As a consequence

x − f l(x) x

≤ 2^−t.

54 / 116

(55)

師大

|x − f l(x)|

|x| = |0.00 · · · 0a_t+1at+2· · · × 2^m|

|1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m|

= |0.a_t+1at+2· · · |

x − f l(x) x

≤ 2^−t.

55 / 116

(56)

師大

|x − f l(x)|

|x| = |0.00 · · · 0a_t+1at+2· · · × 2^m|

|1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m|

= |0.a_t+1at+2· · · |

x − f l(x) x

≤ 2^−t.

56 / 116

(57)

師大

|x − f l(x)|

|x| = |0.00 · · · 0a_t+1at+2· · · × 2^m|

|1.a₁a₂· · · a_ta_t+1a_t+2· · · × 2^m|

= |0.a_t+1at+2· · · |

x − f l(x) x

≤ 2^−t.

57 / 116

(58)

師大

If t-digit rounding arithmetic is used and

a_t+1= 0, then f l(x) = ±1.a1a₂· · · a_t× 2^m.A bound for the relative error is

|x − f l(x)|

|x| = |0.at+1at+2· · · |

|1.a1a₂· · · ata_t+1a_t+2· · · | × 2^−t≤ 2^−(t+1), since the numerator is bounded above by ¹₂ due to at+1= 0.

at+1= 1, then f l(x) = ±(1.a1a2· · · at+ 2^−t) × 2^m.The upper bound for relative error becomes

|x − f l(x)|

|x| = |1 − 0.at+1a_t+2· · · |

|1.a1a2· · · atat+1at+2· · · | × 2^−t≤ 2^−(t+1), since the numerator is bounded by ¹₂due to at+1= 1.

Therefore the relative error for rounding arithmetic is

x − f l(x) x

≤ 2^−(t+1)= 1 2× 2^−t.

58 / 116

(59)

師大

a_t+1= 0, then f l(x) = ±1.a1a₂· · · a_t× 2^m. A bound for the relative error is

|x − f l(x)|

|x| = |0.at+1at+2· · · |

at+1= 1, then f l(x) = ±(1.a1a2· · · at+ 2^−t) × 2^m. The upper bound for relative error becomes

|x − f l(x)|

|x| = |1 − 0.at+1a_t+2· · · |

x − f l(x) x

≤ 2^−(t+1)= 1 2× 2^−t.

59 / 116

(60)

師大

|x − f l(x)|

|x| = |0.at+1at+2· · · |

at+1= 1, then f l(x) = ±(1.a1a2· · · at+ 2^−t) × 2^m.The upper bound for relative error becomes

|x − f l(x)|

|x| = |1 − 0.at+1a_t+2· · · |

x − f l(x) x

≤ 2^−(t+1)= 1 2× 2^−t.

60 / 116

(61)

師大

|x − f l(x)|

|x| = |0.at+1at+2· · · |

|x − f l(x)|

|x| = |1 − 0.at+1a_t+2· · · |

x − f l(x) x

≤ 2^−(t+1)= 1 2× 2^−t.

61 / 116

(62)

師大

|x − f l(x)|

|x| = |0.at+1at+2· · · |

|x − f l(x)|

|x| = |1 − 0.at+1a_t+2· · · |

x − f l(x) x

≤ 2^−(t+1)= 1 2× 2^−t.

62 / 116

(63)

師大

Definition 11 (Machine epsilon)

The floating-point representation, f l(x), of x can be expressed as

f l(x) = x(1 + δ), |δ| ≤ ε_M, (1) whereεM ≡ 2^−t is referred to as theunit roundoff erroror machine epsilon.

Single precision IEEE standard floating-point format The mantissa f corresponds to 23 binary digits (i.e., t = 23), the machine epsilon is

εM = 2⁻²³≈ 1.192 × 10⁻⁷.

This approximately corresponds to7accurate decimal digits

63 / 116

(64)

師大

εM = 2⁻²³≈ 1.192 × 10⁻⁷.

64 / 116

(65)

師大

εM = 2⁻²³≈ 1.192 × 10⁻⁷.

65 / 116

(66)

師大

Double precision IEEE standard floating-point format The mantissa f corresponds to 52 binary digits (i.e., t = 52), the machine epsilon is

εM = 2⁻⁵²≈ 2.220 × 10⁻¹⁶.

which provides between15and16decimal digits of accuracy.

Summary of IEEE standard floating-point format

single precision double precision

ε_M 1.192 × 10⁻⁷ 2.220 × 10⁻¹⁶

smallest positive number 1.175 × 10⁻³⁸ 2.225 × 10⁻³⁰⁸ largest number 3.403 × 10³⁸ 1.798 × 10³⁰⁸

decimal precision 7 16

66 / 116

(67)

師大

εM = 2⁻⁵²≈ 2.220 × 10⁻¹⁶.

ε_M 1.192 × 10⁻⁷ 2.220 × 10⁻¹⁶

67 / 116

(68)

師大

εM = 2⁻⁵²≈ 2.220 × 10⁻¹⁶.

ε_M 1.192 × 10⁻⁷ 2.220 × 10⁻¹⁶

68 / 116

(69)

師大

Let stand for any one of the four basic arithmetic operators +, −, ?, ÷.

Whenever twomachine numbersxand y are to be combined arithmetically, the computer will produce f l(x y)instead of x y.

Under (1), the relative error of f l(x y) satisfies

f l(x y) = (x y)(1 + δ), δ ≤ ε_M, (2)

where ε_M is the unit roundoff.

But if x, y arenotmachine numbers, then they must first rounded to floating-point format before the arithmetic operation and the resulting relative error becomes

f l(f l(x) f l(y)) = (x(1 + δ₁) y(1 + δ₂))(1 + δ₃),

where δ_i ≤ ε_M, i = 1, 2, 3.

69 / 116

(70)

師大

f l(x y) = (x y)(1 + δ), δ ≤ ε_M, (2)

f l(f l(x) f l(y)) = (x(1 + δ₁) y(1 + δ₂))(1 + δ₃),

where δ_i ≤ ε_M, i = 1, 2, 3.

70 / 116

(71)

師大

f l(x y) = (x y)(1 + δ), δ ≤ ε_M, (2)

f l(f l(x) f l(y)) = (x(1 + δ₁) y(1 + δ₂))(1 + δ₃),

where δ_i ≤ ε_M, i = 1, 2, 3.

71 / 116

(72)

師大

f l(x y) = (x y)(1 + δ), δ ≤ ε_M, (2)

f l(f l(x) f l(y)) = (x(1 + δ₁) y(1 + δ₂))(1 + δ₃),

where δ_i ≤ ε_M, i = 1, 2, 3.

72 / 116

(73)

師大

Example

Let x = 0.54617 and y = 0.54601. Using rounding and four-digit arithmetic, then

x^∗= f l(x) = 0.5462is accurate tofoursignificant digits since

|x − x^∗|

|x| = 0.00003

0.54617 = 5.5 × 10⁻⁵ ≤ 5 × 10⁻⁴.

y^∗ = f l(y) = 0.5460is accurate tofivesignificant digits since

|y − y^∗|

|y| = 0.00001

0.54601 = 1.8 × 10⁻⁵≤ 5 × 10⁻⁵.

73 / 116

(74)

師大

Example

|x − x^∗|

|x| = 0.00003

0.54617 = 5.5 × 10⁻⁵ ≤ 5 × 10⁻⁴.

|y − y^∗|

|y| = 0.00001

0.54601 = 1.8 × 10⁻⁵≤ 5 × 10⁻⁵.

74 / 116

(75)

師大

Example

|x − x^∗|

|x| = 0.00003

0.54617 = 5.5 × 10⁻⁵ ≤ 5 × 10⁻⁴.

|y − y^∗|

|y| = 0.00001

0.54601 = 1.8 × 10⁻⁵≤ 5 × 10⁻⁵.

75 / 116

(76)

師大

The exact value of subtraction is r = x − y = 0.00016.

But

r^∗≡ x y = f l(f l(x) − f l(y)) = 0.0002.

Since

|r − r^∗|

|r| = 0.25 ≤ 5 × 10⁻¹ the result has onlyonesignificant digit.

Loss of accuracy

76 / 116

(77)

師大

The exact value of subtraction is r = x − y = 0.00016.

But

r^∗≡ x y = f l(f l(x) − f l(y)) = 0.0002.

Since

|r − r^∗|

|r| = 0.25 ≤ 5 × 10⁻¹ the result has onlyonesignificant digit.

Loss of accuracy

77 / 116

(78)

師大

Loss of Significance

One of the most common error-producing calculations involves the cancellation of significant digits due to the subtraction of nearly equal numbersor theaddition of one very large number and one very small number.

Sometimes, loss of significance can be avoided by rewriting the mathematical formula.

Example 12

The quadratic formulas for computing the roots of ax²+ bx + c = 0, when a 6= 0, are

x₁= −b +√

b²− 4ac

2a and x₂= −b −√

b²− 4ac

2a .

Consider the quadratic equationx²+ 62.10x + 1 = 0and discuss the numerical results.

78 / 116

(79)

師大

Example 12

x₁= −b +√

b²− 4ac

2a and x₂= −b −√

b²− 4ac

2a .

79 / 116

(80)

師大

Example 12

x₁= −b +√

b²− 4ac

2a and x₂= −b −√

b²− 4ac

2a .

80 / 116

(81)

師大

Solution

Using the quadratic formula and 8-digit rounding arithmetic, one can obtain

x₁ = −0.01610723 and x₂ = −62.08390.

Now we perform the calculations with 4-digit rounding arithmetic. First we have

pb²− 4ac =p

62.10²− 4.000 =√

3856 − 4.000 = 62.06, and

f l(x₁) = −62.10 + 62.06

2.000 = −0.04000

2.000 = −0.02000.

The relative error in computing x₁is

|f l(x₁) − x1|

|x₁| = | − 0.02000 + 0.01610723|

| − 0.01610723| ≈ 0.2417 ≤ 5×10⁻¹.

81 / 116