Math'375 Fall 2005 I Numbers and Formats

1
Math.375 Fall 2005I - Numbers and Formats

• Vageli Coutsias

2
Introduction

• 1 1 0 or machine epsilon?
• eps 2.220446049250313e-016
• How does matlab produce its numbers?
• Where we learn about number formats, truncation
  errors and roundoff

3
Matlab real number formats
format long (default for ? ) pi
3.14159265358979 format short pi
3.1416 format short e pi
3.1416e000 format long e pi
3.141592653589793e000
4
Floating-point numbers

• Base or radix
• t Precision
• L,U Exponent range

5
The set F of f.p. numbers

• Basis b
• Significant digits t
• Range (U, L)
• F(b, t, U, L) F(2, 53, -1021, 1024)
• is the IEEE standard

6
sign bit
mantissa (52 bits) (it is assumed d1 1, so only
need d2,d53)
exponent (11 bits)
7
IEEE double precision standard If E2047 and F
is nonzero, then VNaN

("Not a number") If E2047 and F0 and S0,(1)
then VInf, (-Inf) If E0 and F0 and
S0,(1), then V0,(-0) If 0ltElt2047 then
V(-1)S 2 (E-1023) (1.F)
where "1.F" denotes the binary number
created by prefixing F with an
implicit leading 1 and a binary point. If E0
and F is nonzero, then V(-1)S
2 (-1022) (0.F) These are
"unnormalized" values.
8
contains unnormalized values, as allowed, e.g.,
in IEEE standard
9
Floating point numbers
The machine precision is the smallest nuber ?
such that
10
Machine epsilon

• The distance from 1 to the next larger float
• Gives the relative error in representing a real
  number in the system F

(RELATIVE) ROUNDOFF ERROR
11
Machine epsilon computed
a 1 b 1 while ab a b b/2 end
b b 1.110223024625157e-016 shows that ab
a is satisfied by numbers b not equal to 0
here b eps/2 is the largest such number for a
1
12

• Overflow does not only cause programs to crash!
• Arianne Vs short maiden flight on 7/4/96 was due
  to a floating exception.

13
FLOAT ? INTEGER

• During the conversion of a 64-bit floating-point
  number to a 16-bit signed integer
• Caused by the float being outside the range
  representable by such integers
• The programming philosophy employed did not guard
  against software errors-a fatal assumption!

14
COMPLEX NUMBERS

• z x iy
• x Re(z) is the real part
• y Im(z) is the imaginary part
• i2 -1 is the imaginary unit
• polar form
• complex conjugate

15

• Matlab commands
• gtgt z 3i4
• gtgt Cartesian form
• x real(z) y imag(z)
• gtgt Polar form
• theta angle(z) rho abs(z)
• So z abs(z)(cos(angle(z))
  isin(angle(z)))
• x-iy conj(z)

16
The complex plane
z12i, 3,-1i, -i compass(z,r) defines
an array of complex numbers which are plotted as
vectors in the x-y plane
17
Roots of complex numbers
x -1 x(1/3) ans 0.5000
0.8660i Matlab assumes complex arithmetic, and
returns automatically the root with the smallest
phase angle ( the other two roots are
-1 and 0.5000 - 0.8660i
18
Types of errors in numerical computation

• Roundoff errors
• Pi 3.14159
• Pi 3.1415926535897932384626
• Truncation errors
• Cos x 1 x2/2
• Cos x 1 x2/2 x4/4!
• Errors usually accumulate randomly
• (random walk)
• But they can also be systematic, and the reasons
• may be subtle!

19
x linspace(1-210-8,1210-8,401) f
x.7-7x.621x.5-35x.435x.3-21x.27x-1
plot(x,f)
CANCELLATION ERRORS
20
g -1x.(7x.(-21x.(35x.(-35x.(21x.(-7
x)))))) plot(x,g)
21
h (x-1).7 plot(x,h)
22
plot(x,h,x,g,x,f)
23
z(1) 0 z(2) 2 for k 229 z(k1)
2(k-1/2)(1-(1-4(1-k)z(k)2)(1/2))(1/2) end
semilogy(130,abs(z-pi)/pi)
24
Recurrence relations
25
Instability is inescapable But one can learn to
ride it!
26
(No Transcript)
27
(No Transcript)
28
Due to limitations in computer arithmetic, need
to practice defensive coding.

• Mathematics?numerics implementation
• The sequence of computations carried out for the
  solution of a mathematical problem can be so
  complex as to require a mathematical analysis of
  its own correctness.
• How can we determine efficiency?

Text
29
Vocabulary

• Consistency correctness -
• Convergence (rate)
• Efficiency (operation counts)
• Numerical Stability (error amplification)
• Accuracy (relative vs. absolute error)
• Roundoff vs. Truncation

30
The first monster
for k150, term xterm/k s stermerr(k)
abs(f(k) - s)end relerr err/exp(x)
semilogy(1nTerms,relerr) semilogy(1nTerms,err)
ylabel('Relative Error in Partial Sum.')
xlabel('Order of Partial Sum.')
title(sprintf('x 5.2f',x)) figure
semilogy(1nTerms,err)end term 1 s 1 f
exp(x)ones(nTerms,1) for k150, term
xterm/k s sterm err(k)
abs(f(k) - s) end relerr err/exp(x)
semilogy(1nTerms,relerr) semilogy(1nTerms,err)
ylabel('Relative Error in Partial Sum.')
xlabel('Order of Partial Sum.')
title(sprintf('x 5.2f',x)) figure
semilogy(1nTerms,err)end

• T(x) ex log(1e-1)
• For xgt30 a dramatic deviation from the
  theoretical value 1 is found.
• For xgt36.7 the plotted value drops (incorrectly)
  to zero
• Here
• x linspace(0,40,1000)
• y exp(x).log(1exp(-x))

31
(No Transcript)
32
(No Transcript)
33
(No Transcript)
34
Pet Monster I
35
(No Transcript)
36
(No Transcript)
37
(No Transcript)
38
Summary

• Roundoff and other errors
• Multiple plots
• Formats and floating point numbers

39
References

• C. Essex, M.Davidson, C. Schulzky,
• Numerical Monsters, preprint.
• Higham Higham, Matlab Guide, SIAM
• SIAM News, 29(8), 10/98 (Arianne V failure)
• B5 Trailer http//www.scifi.com/b5rangers/
• As usual, find m-files at
• aix coutsias/375/IV