CS 314

July 19th  – July 23rd

 

Grades are in mentor (so far)

 

Login to mentor

Cd ~cs314/Summer2004/MAILBOX/$user

Ls

Cat grades

Hw1   hw2-part1 midterm

 

Homework 4 is posted

“Curve fitting polynomials and splines”

 

Midterm solution is posted as well

 

Numerical Differentiation

 

Limit of Difference Quotient

 

f’(x) = lim h-> 0 (f(x+h) – f(x)) / h

 

we can approximate the derivative by choosing a small h

 

How small does h have to be?

 

It depends on f(x).

 

We can start with a large h and make it smaller until we see that Dk (the approximation) is losing precision.

 

Example: | D­k – Dk-1 | < Epsilon then stop

 

Dk = f(x+hk) – f(x) / hk for k = 1, 2, 3…

 

Example: f(x) = ex   - we want to approximate f(x)

 

             hk                           Dk

0          .1                     2.8588

1          .01                   2.7319

2          .001                 2.7196

3          .0001               2.7184

4          .00001             2.7183

5          .000001           2.7183

6          .0000001         2.7183

 

hk = 10-(k+1)

 

Dk = ex + hk – ex / hk

At x= 1 -ΰ Dk = (e1 + h – e) / hk

 

Let Epsilon = 0.0001 then we stop at k = 5

 

The iteration approach to obtain the derivate is slow. There are other alternatives.

 

Central Differences Formulas

 

Centered Formula of order O(h2)

  Start with the Taylor expansion

 

1) F(x+h) = f(x) + f’(x)h + f’’(x)h2 / 2! + f’’’(c1)h3 / 3!

 

F(x-h) = f(x) + f’(x)(-h) + f’’(x)(-h)2 /  2! + f’’’(x)(-h)3 / 3!

 

2) F(x-h) = f(x) - f’(x)h  + f’’(x)h2 / 2! - f’’’(c2)h3 / 3!

 

Subtracting 2 from 1

 

F(x+h) – f(x-h) = f(x) + f’(x)h + f’’(x)h2 / 2! + f’’’(x)h3 / 3! – (f(x) - f’(x)h  + f’’(x)h2 / 2! - f’’’(x)h3 / 3!)

 

F(x+h) – f(x-h) = 2f’(x)h + (f ‘’’(c1) + f’’’(cs) h3) / 3!  Truncation error

 

f’(x) = f(x+h) – f(x-h) / 2h

 

 

Example: f(x) = ex h = 0.0001

 

F’(1) = e1+0.0001 – e1-0.0001 / 2(0.0001) = 2.7183

----------------------------------------------------------------------------------------------

 

Centered formula of order O(h4)

 

Do Taylor expansion:

 

f(x+h) = f(x) + f’(x) = f’’(x)h2 / 2! = f’’’(x)h3 / 3! + f(4)(x)h4 / 4! + f(5)(c1) h5 / 5!

 

F(x-h) = f(x) – f’(x)h + f’’(x)h2 / 2! – f’’’(x)h3 / 3! + f(4)(x)h4 / 4! – f(5)(c2) / 5!

 

1) F(x+h) – f(x-h) = 2f’(x) 2f’’’x) h3 / 3! + f(5)(c1) h5 / 5! + f(5)(c2) / 5!

 

Using step size 2h instead of h in (1) we get:

f(x+2h) – f(x-2h) = 2f’(x)(2h) + 2f’’’x) (2h)3 / 3! + f(5)(c1) (2h)5 / 5! + f(5)(c2) / 5!

 

f(x+2h) – f(x-2h) = 4f’(x)h + 16f’’’(x) h3 / 3! + 32f(5)(c1) h5 / 5! + f(5)(c2) / 5!

 

Now multiply (1) by 8 and subtract (2):

 

F(x+h) – 8f(x-h) – (f(x+2h) - f(x-2h) = 1f’(x)h + 16f’’’(x)h3 / 3! + 8f5(c1)f5(c2) / 5!

 

8f(x+h) – 8f(x-h) – f(x+2h) – f(x-2h) = 12f’(x)h + kf(5) h5

 

So we have:

 

F’(x) = -f(x+2h) + 8f(x+h) – 8f(x-h) + f(x-2h) / 12h + kf9v)(c3)h5 / 12h

 

 

Numerical Integration

 

Trapezoidal Rule – approximates area under the curve with trapezoids

 

A = (f0 + f1 / 2) * h for each section of the curve

 

A1 + A2 … = (f0 + f1 / 2) * h + (f1 + f2 / 2) * h   …..

 

= h * (f0 / 2 + f1 / 2 + f1 / 2 + f2 / 2 + f2 / 2 + f3 / 2 + f3 / 2)

 

= h * (f0 / 2 + f1 + f2 + f3 / 2)

 

 

For example int( e-x^2)    m = 4

 

H = b – a / M = (1) – (-1) / 4 = 2/4 = ½ = .5

 

int( e-x^2) = .5 [(e-(-1)^2 + e-(1)^2]) / 2 + e-(-.5)^2] + e-(0)^2] + e-(0.5)^2]

 

= 1.4627

 

 

int(sin(x) from 0 to pi/2

 

= pi / 6 * ((sin(o) + sin(pi/2)) / 2 + sin(pi/6) + sin(pi/3)

 

= .9770486j

 

exact solution:

 

int(sin(x) from 0 to pi/2 = - cos(x) from pi/2 to 0 = -cos(pi/2) – cos(o) = 1

 

 

 

Simpson Rule – It approximates the integral using the area under a parapola every three points.

 

The quadratic polynomial gives a better approximation than the line used in the trapezoidal rule.

 

Using Lagrange Polynomials we can obtain P2(x):

 

P2(x) = (f0 * (x – x1) (x – x2) / (x0 – x1) (x0 – x2)) + (f1 * (x – x0) (x – x2) / (x1 – x0) (x1 – x2)) + (f2 * (x – x0) (x – x1) / (x2 – x0) (x2 – x1))

 

Int(P2(x)) from x2 to x1 = int( ((f0 * (x – x1) (x – x2)) / (-h)(-2h))  + ((f1 * (x – x0) (x – x2)) / (h)(-h)) + ((f2 * (x – x0) (x – x1)) / (2h)(h))

 

To do the integral, we do the substitution x = x0 + ht

 

T = x – x0 / h

 

Dx = h dt

 

Int(P2(x)) from 0 to 2 = Int( (f0 * (x0 + ht – x1)(x0 + ht – x2) / 2h2) + (f1 * (x0 + ht – x0)(x0 + ht – x2) / -h2) + (f3 * (x0 + ht – x0)(x0 + ht – x1) / 2h2)

 

Int(P2(x)) from 0 to 2 = (f0 * (ht –h )(ht – 2h) / 2h2) + f1 ((ht)(ht-2h)/(-h2) + f2 * (ht)(ht –h )(2h2) hdt

 

Int(P2(x)) from 0 to 2 = (f0 / 2)(t -1)(t-2) + (f1 / -1)(t)(t-2) + (f2 / 2)(t)(t-1)

 

= f0h / 2 [ t3 – 3t2 + 2t]  – f1h [t3/3 – 2t2/2] + f2h [t3 / 3 – t2/2] from 0 to 2

 

With all the math is comes out to:

 

Int(P2(x)) = h/3 * [f0 + 4f1 + f2]

 

This equation is only for three points. If the interval [a,b] is subdivided into 2M intervals [xk, xk+1] of equal width we can use the simpson rule every three points.

 

 

A = A1 + A2 + A3

 

A = h/3 * [f0 + 4f1 + f2] + h/3 * [f2 + 4f3 + f4] + h/3 * [f4 + 4f5 + f6]

 

A = h/3 [f0 + 4f1 + 2f2 + 4f3 + 2f4 + 4f5 + f6]

 

In general:

 

A = h/3 [f0 + 4f1 + 2f2 + 4f3 + ... 4f2m-1 + f2m]

 

Example:

 

Int(sin(x)) from 0 to pi/2 w/ 2m = 2

 

H = b-a / 2m

 

H= pi/2 – 0 / 2(1) = pi / 4

 

A = pi/4 / 3 [sin(0) + 4sin(pi/4) + sin(pi/2)]

 

A = 1.002279878

 

exact solution:

 

int(sin(x) from 0 to pi/2 = - cos(x) from pi/2 to 0 = -cos(pi/2) – cos(o) = 1

 

Now using 2m = 4

 

H = b – a / 2m = (pi/2 – 0) / 4 = pi/8

 

A = (pi/8) / 3 [sin(o) + 4sin(pi/8) + 2sin(pi/4) + 4sin(3pi/8) + sin(pi/2)]

 

A = 1.000134585

 

 

The Golden Ratio

 

Assume p(x) has one minimum in the interval [a,b], we divide the interval into 3 sub intervals a < c < d < b

 

We can evaluate f(x) at c and d

 

If f(c) < f(d), then the new interval will be reduced to [a,d]

 

If f(c) > f(d), then the new interval will be [c,b]

 

[c,b]  would be any 2 values such that a < c < d < b, we choose [c,d] to such a way that we reduce the number computations done.

 

By using golden ratio to choose c, d. we reduce the number of operations that we need to do.

 

If the new interval becomes [c,b] then old c will be new a, and the old d becomes new c. The new d will be have to be recomputed.

 

So, at every iteration, we need to recomputed only 1 value either f(c) or f(d) instead of 2 values.

 

Example:

 

i = 1     a1 = -2   b1 = 1    

 

c1 = a + (1-r)(b-a)

=2 + (1- .61803)(1-(-2)

=-.14 589

 

d1 = a + r(b-a)

=-2 + .61803  (1-(-2))

=-.14589

 

f(a) = f(-2) = (-2)2 +1 = 5

f(b) = f(1) = (1)2 + 1 = 2

f© = f(-.8541) = 1.72949

f(d) = f(-.14589) = 1.08128

f(d) < f(c) ΰ the new interval [ c, b] = [-.8541, 1]

 

old c becomes a

old d becomes c

 

i=2

 

a2 = -.8541     b2 = 1

 

c2 = a + (1-r)(b-a)

= -.8541 + (1-.61803)(1-(-.8541))

=-.14589 = d1 as e3xpected

d2 = a + r(b-a) = -/8541 + .61803(1-(-.8541)

= .291789

 

f(a2) = f(-.8541) = 1.72949

f(b2) = f(1) = 2

f(c2) = f(-.14589) = 1.02128

f(d2) = f(.291789) = (.291789)2 +1 = 1.08514

f(c2) < f(d2) – The minimum is in the left side of d2

 

a3 <- a2

b3 <- b2

 

i=3

 

a3 = -.8541

b3 = -.291789

c3 = a+(1-r)(b-a)

=-.8541 + (1- .61803)(.291789 – (-.8541)

=-.4164

d3 = c2 = -.14584

 

f(a3) = f(a2) = f(-.8541) = 1.72949

f(b3) = f(d2) = f(-.291789) = 1.08514

f(c3) = f(-4164) = 1.17339

f(d3) = f(-1.4584) = 1.02128

 

f(c3) > f(d3)

 

new interval

a4 <- c3

b4 <- b3                  [-.4164, .291]

 

 

 

 

Steepest Descent or Gradient Method

 

-Assume we want to minimize f(x) of N variables where x = (x1, x2, x3 …)

 

The gradient f(x) is a vector function defined as

 

Upsidedown trianfle f(x) = del(fx) / del(x1) + del(fx) / del(x2) + del(fx) / del(x3) ...

 

From the concept of gradient we know that the gradient vector points in the direction of greatest increase of f(x)

 

Then –gradf(x) will point to the direction of greatest decrease.

 

To obtain the minimum using the gradient method, starting at point P we move along the line in the direction of greatest decrease, that is –gradf(x)

In the simplest form

 

P1 = Po – G0h

..

..

Pk+1 = Pk – Gkh

 

Example:

 

Y = x2 +1

Y’ = 2x                G = 2x

G0 = 2(-1) = -2

 

X1= x0 – g0h

=-1 – (-2)(.1)

=-1+.2=-.8

 

G1 = 2(x1) = 2(-.8) = -1.6

X2  = x1 – G1h

=-.8-(-1.6)(.1)

=.63

 

G2 = 2(-.64) = -1.28

X3 = -.64-(-1.28)(.1)3

=-.512

 

G3 = 2(-.512) = -1.034

X4 = -.512 – (-1.024)(.1)

=-.4096

 

Ginf  = 0

Xinf = 0

 

 

 

Numerical Methods for Differential Equations

 

Consider the differential equation

Dy/dt = ky

Ln(y) = k1t + k2

 

Some differential equations do not have an analytical solution so they have to be approximated using numerical methods.

 

Euler’s Method

Let [a,b] be the interval over which we want to find the solution y’ = f(t,y)