Numerical Methods Week #1

Monday, June 14, 2004

Goals:

· Learn how to use the computer to solve and approximate linear/non-linear equations, numerical integration and differentiation, interpolation and differential equations

· The limitations as well as the potential of the computer to solve mathematical problems

Webpage: http://www.cs.purdue.edu/homes/cs314

Office Hours: Walk in for short questions or make an appointment. 
E-mail recommended.

Textbook: Numerical Methods using Matlab (3rd. ed.ition) By John H. Mathews and Kurtis D. Fink (Prentice Hall)

Grade Distribution:

Midterm, Final 50%

Projects and Homework 50%

Mailing List:

· Login to a CS machine (Lore)

· If you don’t have a CS account, send e-mail to grr@cs.purdue.edu

· Type “mailer add me to cs314”

· Send questions to cs314-ta@cs.purdue.edu. Don’t post it to the list.

· Answers will be posted depending on the importance of it to everyone

Newsgroups:

· Add yourself to Purdue Newsgroups CS314 Class

(Professor won’t be reading it. But, TA’s will be)

Syllabus:

· Representation of Floating Point Numbers and Error Analysis

· Solution of Nonlinear Equations

· The Solution of Linear Systems AX=B

· Interpolation and Polynomial Approximation

· Curve Fitting

· Numerical Differentiation

· Numerical Integration

· Numerical Optimization

· Solution of Differential Equations

Lecture Notes:

· Volunteers can take notes for one week and it will count up to 2% extra over the final grade depending on the quality of the class notes

· Notes in a HTML file and JPG/GIF format for diagrams

· E-mail it to zouh@cs.purdue.edu

· Do not e-mail them to the professor

· Sign up will be tomorrow, Tuesday June 15, 2004.

Class Time: MTWTF 1.00 PM – 2.00 PM at UNIV 303

Binary Numbers:

· Computers use binary numbers

· Easier to implement circuits with 2 states – “On” and “Off”

· Integer Representation:

o b0 * 2⁰ + b1 * 2¹ + b2 * 2² ……

Tuesday, June 15, 2004

NOTE: Homework 1 Available on the webpage due Tuesday, June 22, 2004 at 11.59 PM electronically

Decimal-Binary Conversion:

Separate integer part from the fraction: N.F

N = b₀ + (b₁x2¹) + (b₂x2²) + ... + (b_kx2^k)

We divide by 2 in order to find b_i (0 or 1).

N/2 = b₀/2 + (b₁x2⁰) + ... + (b_kx2^k-1)

The remainder is b₀/2. The integer part is the rest, which is (b₁x2⁰) + ... + (b_kx2^k-1)

After this operation, the new N is (b₁x2⁰) + ... + (b_kx2^k-1)

The remainder determines whether b₀ is 0 or 1.

The above process is repeated until the integer part is 0.

Example: 57₁₀

57/2 = 28, remainder = 1

28/2 = 14, remainder = 0

14/2 = 7, remainder = 0

7/2 = 3, remainder = 1

3/2 = 1, remainder = 1

1/2 = 0, remainder = 1

The answer is 111001

0.1₂ = 1x2^-1 = ½ = 0.5₁₀

0.01₂ = 0x2^-1 + 1x2^-2 = ¼ = 0.25₁₀

1/3 ₁₀ = 0x2^-1 + 1x2^-2 + 0x2^-3 + 1x2^-4 + …. = 0.0101010101…₂

Representation of Numbers in Machine

• Fixed Point and Floating Point
• Fixed Point:

The number of bits for the integer part and the fraction part are fixed.

Example: NNNN.FF₂ has 4 bits for the integer part and 2 bits for the fraction. In this case, the largest number is 1111.11₂ = 15.75_{10 and the} smallest Number – 0000.00₂ = 0₁₀

Advantages:

Fast arithmetic - Uses integer arithmetic

Used in signal processors and signal processing programs that requires real time operations. Example: MP3 Player, MP3 Encoder/Decoder, MPEG, etc

Disadvantages:

Small range of numbers

Bad precision

Have to be careful about potential overflow

• Floating Point:

The decimal point moves around to keep desired precision.

Uses scientific notation:
where Q is called the mantissa and N is called the exponent

Precision of the number depends on the number of the bits used to represent Q.

x = ±Q * 2^N [Numbers are represented in the computer]

Q, N: represented with a fixed number of bits.

To prevent multiple representations, Q and N are normalized so that there is only one representation.

Example: Making 0.5 ≤ q < 1, other than q = 0

Representation of Real Numbers

• Float and Double
• Float (Single Precision)

4 bytes (32 bits)

24 bits mantissa (1 bit sign)

8 bits exponent (1 bit sign)

Range is 2^-128 to 2¹²⁷

Precision is approximately 6 digits

• Double

8 bytes (64 bits)

53 bits mantissa (1 bit sign)

11 bits exponent (1 bit sign)

Range is 2 ¹⁰²⁴

Precision is approximately 15 digits

Errors

• Absolute Error = | p – p* | where p* is the approximation of p, p is quality
• Relative Error = | p – p* | / | p | for P ≠ 0

Truncation vs. Rounding

Use three digits for fraction.

Absolute Error

          Chop-off = 0.59 x 10^-3

          Round-off = 0.41 x 10^-3

Relative Error

          Chop-off = 1.9 x 10^-4

          Round-off = 1.3 x 10^-4

Conclusion: Due to a smaller error, Round-off is better.

Propagation Of Error

• Assume p and q are approximated by p* and q*
• p* = p + ep
• q* = q + eq
• Sum: p* + q* = p + q + (ep + eq). Sum of errors = new error
• Product: p*q* = pq + p(eq) + q(ep) + (ep)(eq) New error is amplified by the magnitude of p and q.

An algorithm is stable if a small initial error stays small and unstable if a small initial error gets large.

Solution of Nonlinear Equations (f(x) = 0)

• Iteration to solve these equations
• We start with a value P₀ and a function P_k+1 = g(P_k)
• Stopping criteria: | P_k - P_k+1 | < E

· Fixed Point: a point p such that p = G(p)

· Fixed Point iteration: the sequence of values produced by p(n-1) = G(p(n)) for n = 0, 1, 2,

· Example:

x² -2x +1 = 0

x² +1 = 2x => x = (x² +1)/2

x_k+1 = (x_k² +1)/2

• Fixed Point Theorem

If | g’(x) | ≤ k < 1, for x in [a,b] and g(x) in [a,b], then the iteration P_n = g(p_n-1) will converge to a unique fixed point. In this case P is said to be an attractive fixed point

If | g’(x) | > 1, then the iteration P_n = g(p_n-1) will not converge

Convergence

• If g'(x) < 1 => convergence
• If g'(x) ≥ 1 => divergence
• Fixed points are the intersection of y = x and y = g(x), achieved by repeated iteration until some minimum error bounds around the point are achieved.
• Oscillating convergence occurs when -1 < g'(p) < 0
  Monotone convergence occurs when 0 < g'(p) < 1
  Monotone divergence occurs when 1 < g'(p)
  Oscillating divergence occurs when -1 > g'(p)
• Example: x² -2x +1 = 0

g(x) = (x² +1)/2

g(0) = 1/2
g(1) = 1
g'(x) = x
g'(x) < 1. Therefore, it is convergence.

Bisection Method:

Used to find a solution for f(x) = 0

Start with two values a, b (an interval) where a < b and such that

• f(a) <0 and f(b) > 0 or f(a) > 0 and f(b) < 0

Similar to binary search.

Method:

• Determine middle point, c = (a+b) /2
• If f(a) and f(c) have opposite signs then ‘0’ lies in [a,c]. Then substitute c for b => b ← c
• If f(b) and f(c) have opposite signs then ‘0’ lies in [c,b]. Then substitute c for a => a ← c
• If f(c) = 0 then the solution has been found as c is zero

Example:

sin(x) = ½ (x is in degrees)

f(x) = sin(x) – ½ = 0

We start with a = 0, b = 50

a	b	c	f(c)
0	50	25	-0.077
25	50	37.5	0.1087
25	37.25	31.25	0.018
25	31.25	28.125	-0.02

C approaches 30.

At each iteration, the range is reduced by half. If we want an approximation error < E, we need n iterations such that log[2] |b - a|/e ≤ n

Example:
sin(x) = ½ = 0, e = .01 n > 12
The solution will be at most 0.01 units from the approximation obtained if at least 13 iterations are performed.

False Position (Regula-Falsi)

• Similar to bisection method.
• Needs an interval [a, b] where a zero is fast
• Converges faster than bisection
• A line is subtended between (a, f(a)), (b, f(b))
• ‘c’ is the intersection of this line and x axis  
• If f(a) and f(c) have different signs, then assign b←c
• If f(b) and f(c) have different signs, then assign a←c

• If f(c) = 0 then the solution has been found as c is zero

• To compute c:

m = (f(a) – f(b)) / (a-b) ---> 1

m = (f(b) – 0) / (b-c)   ----> 2

1 = 2

f(b) / (b-c) = (f(a) – f(b)) / (a-b)

b-c = f(b) * (a-b) / (f(a) – f(b))

c = b - f(b) * (a-b) / (f(a) – f(b) )

Example:

sin(x) = ½  => f(x) = sin(x) – ½ = 0

Let a = 0 , b = 50

When a = 0, b = 50

f(a) = -0.5, f(b) = .766, c = 37.546, f(c) = .03
When a = 0, b = 37.546

f(a) = -0.5, f(b) = .03, c = 30.70, f(c) = .01

When a = 0, b = 30.70

f(a) = -0.5, f(b) = .01, c = 30.09, f(c) = .001

C seems to be approaching zero.

• Disadvantages

The approximation to a line using two points may be bad

Need an interval a, b where the solution is found

Newthon-Rapson

• Uses only one value as guess instead of on interval
• Uses the derivative

Convergence

• Horizontal Convergence

we expect f(x) to be between some small error E such that f(x) < E

• Vertical Convergence

We stop when | x – p | < S.

We already know p. Thus, this is not possible.

This is approximated by determining when | P_n – P_n-1 | < S

• Both Horizontal and Vertical Convergence

Stop when | P_n – P_n-1 | < E & f(x) < E

Troublesome Functions

If the graph y=f(x) is steep near the root (P, 0) then the root finding is well conditioned. Easier to find the solution. Very close to vertical. The intersection with the x-axis will be well-marked. If the graph y=f(x) is shallow near (P, 0) then the root finding may only have a few significant digits. Intersections are not well-marked.