Numerical Differentiation

Numerical Differentiation

On homework #2, we use the limit of difference quotient

f’(x) = lim f(x+h) – f(x)

h->0 h

We can approximate the derivative by using small h.

How small h has to be?

It depends on f(x). We can start with a large h and make it smaller and smaller until we sees that D_k (the approximation) is loosing precision.

Example: |D_k – D_k-1| < E then stop

Dk = f(x + h_k) – f(x) for k = 1, 2, …, N

H_k

Example: f(x) = e^x. We want to approximate f’(x) at x =1

Dk = e^(x+hk)– e^x

h_k

h_k D_k

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), In this case let E = .00001, the we stop at k = 5

This iterative approach to obtain the derivative is slow.

There are alternatives. One is called Central Differences Formulas

Central Differences Formulas of order O(h²)

- start with Taylor expansion: f(x+h) = f(x) + f’(x)h

We want x= x₀+h, therefore, x-x₀ = h

(I) f(x+h) = f(x) + f’(x) h + f’’(x) h² + f’’’(c₁) h³ which f’’’(c₁) h³ is the error

2! 3! 3!

f(x-h) = f(x) + f’(x) (-h) + f’’(x) (-h)² + f’’’(c2) (-h)³

2! 3!

= f(x) - f’(x) (h) + f’’(x) (h)² - f’’’(c2) (h)³ (II)

2! 3!

Subs (II) from (I)

f(x+h) – f(x-h) = 2f’(x) (h) + f’’’(c₁) h³ + f’’’(c2) (h)³

Our objective is to find f’(x)

f’(x) = 1 [ f(x+h) – f(x-h) – [f’’’(c₁) + f’’’(c2) ] h³ ]

2h 3!

f’(x) = f(x+h) – f(x-h) - [f’’’(c₁) + f’’’(c2) ] h²

2h 2.3!

Which [f’’’(c₁) + f’’’(c2) ] h²is the error

2.3!

So, f’(x) ≈ f(x+h) – f(x-h)

Example: f(x) = e^x and h=.0001

f’(x) = e^1+.0001– e^1-.0001 = 2.7183

2*(.0001)

Centered Formula of order O(h⁴)

Do Taylor Expansion:

- f(x+h) = f(x) + f’(x) h + f’’(x) h² + f’’’(x) h³+ f’’’’(x) h⁴+ f’’’’’(c₁) h⁵

2! 3! 4! 5!

Which f’’’’’(c₁) h⁵ is the error.

- f(x+h) = f(x) + f’(x) (-h) + f’’(x) (-h)² + f’’’(x) (-h)³+ f’’’’(x) (-h)⁴ + f’’’’’(c₂) (-h)⁵

2! 3! 4! 5!

Which f’’’’’(c₂) (-h)⁵ is the error.

So, (I) f(x+h) – f(x-h) = 2 f’(x) h + 2 f’’’(x) h³ + [ f’’’’’(c₁) + f’’’’’(c₂) ] (h)⁵

3! 5!

Which [ f’’’’’(c₁) + f’’’’’(c₂) ] (h)⁵ is the error

Using step size 2h instead of f. in (I), we get

f(x+2h) – f(x-2h) = 2 f’(x) 2h + 2 f’’’(x) (2h)³ + [ f’’’’’(c₁) + f’’’’’(c₂) ] (2h)⁵

3! 5!

(II) = 4 f’(x) h + 16 f’’’(x) (h)³ + 32 [ f’’’’’(c₁) + f’’’’’(c₂) ] (h)⁵

3! 5!

Which 32 [ f’’’’’(c₁) + f’’’’’(c₂) ] (h)⁵ is the error

Then multiple (I) by 8 and substract (II)

- 8(f(x+h) – f(x-h)) – (f(x+2h) – f(x-2h)) = 12 f’(x) h – 24 [ f’’’’’(c₁) + f’’’’’(c₂) ] h⁵

So we have f’(x) = 8 f(x+h) – 8 f(x-h) + f(x-2h) – f(x+2h)

12h

+ 24 [ f’’’’’(c₁) + f’’’’’(c₂) ] h⁵ 1

5! 12h

which 24 [ f’’’’’(c₁) + f’’’’’(c₂) ] h⁵ 1 is the error O(h4)

5! 12h

Let f(x) = sin(x); h=.001 centered formula O(h⁴); f’(x) = cos(x); f’(π/3) = cos(π/3) = .5

f’(π/3) ≈ -sin(π/3+.002) +8 sin (π/3+.001) – sin(π/3-.001) + sin(π/3-.002)

12*(.001)

= .5

Using centered formula O(h²) previously

f’(x) ≈ sin(π/3+.001) - sin (π/3-.001) = 4.99999917

2*(.001)

Using The Diff Quotient f’(x) ≈ sin(π/3+.001) - sin(π/3) = .499566904

(.001)

Numerical Integration

Trapezoidal Rule

It approximates the area under the curve with trapezoid

M = 3

A1 = (f₀+f₁).h ; A2 = (f₂+f₁).h ; A3 = (f₂+f₃).h ;

2 2 2

A1 + A2 + A3 = h. [f₀+f₁+f₂+f₃ ]

2 2

In general,

_{b
m-1}

∫ f(x) dx = h [ f(x_o) + f(x_m) + ∑ f(x_k) ] for M+1 pts, x_o, .., x_m

^a 2^K=1

Example:

₁

1) ∫ e^{-x^2} dx M = 4

^-1

h= b-a = 1-(-1) = .5

m 4

₁

∫ e^{-x^2} dx ≈ .5 [e^{-(-1)^2} + e^{-(-1)^2} + e^{-(-.5)^2} + e^{-(0)^2} + e^{-(.5)^2} ] = 1.4627

^-12

_π/2

2) ∫ sin(x) dx M = 3

⁰

h= b-a = π/2-(0) = π/6

m 3

_π/2

∫ sin(x) dx ≈ π/6. [ sin(0) + sin(π/2) + sin(π/6) + sin(π/3) ] = .9770486

⁰2

Exact sol: 1

The Simpson Rule: it approximates the integral using theorem under a parabola every points.

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

We use the Lagrange polynomial, we can obtain.

P₂(x) = f_o (x-x₁) (x-x₂) + f₁ (x-x₀) (x-x₂) + f₂ (x-x₀) (x-x₁)

(x₀-x₁)(x₀-x₂) (x₁-x₀)(x₁-x₂) (x₂-x₀)(x₂-x₁)

x₁-x₀ = h ; x₂-x₀ = 2h ; x₂-x₁ = h ….. (1)

_{x2 x2}

∫ P₂(x) dx = ∫ [ f_o (x-x₁)(x-x₂) + f₁ (x-x₀)(x-x₂) + f₂ (x-x₀)(x-x₁) ] dx

^{x0 x0} 2h² (-h²) 2h²

We use substitution x = x₀+ ht_, therefore t = x - x_0, dx = h.dt

t = x_{0 -} x₀= 0 ; t = x_{2 -} x₀= 2

h h

_{x2 2}

∫ P₂(x) dx = ∫ [ f_o (x₀+ ht - x₁) (x₀+ ht - x₂) + f₁ (x₀+ ht - x₀) (x₀+ ht - x₂)

^{x0 0} 2h² h²

+ f₁ (x₀+ ht - x₀) (x₀+ ht - x₁) ] h dt

2h²

₂

= ∫ [ f_o (x₀+ ht - x₁) (x₀+ ht - x₂) + f₁ (ht) (x₀+ ht - x₂)

⁰2h h

+ f₁ (ht) (x₀+ ht - x₁) ] dt

from (1):

₂

= ∫ [ f_o (ht - h) (ht – 2h) + f₁ (ht) (ht – 2h) + f₁ (ht) (ht – h) ] dt

⁰2h h 2h

₂

= f_o(ht³ – 3ht² + 2ht) - f₁(ht³ – 2ht²) + f₂ (ht³-ht²) |

2 3 2 3 2 2 3 2 ⁰

= h . (f_o – 4 f₁ + f₂ )

This equation is only for 3 pts. If interval [a,b] is subdivided into 2M intervals [x_k, x_k-1] of equal width, we can use the Simpson Rule every 3 pts.

Example: 2M = 6 ; we use 2M so that it is even

M = 3

M is the # of parabolas

So 2M intervals & 2M+1 pts and M parabolas

So A = A1 + A2 + A3 = h (f_o + 4f₁ + f₂) + h (f₂+ 4f₃ + f₄) + h (f₄ + 4f₅ + f₆)

3 3 3

= h (f_o + 4f₁ + 2f₂+ 4f₃ + 2f₄+ 4f₅ + f₆)

In general:

A= h (f_o + 4f₁ + 2f₂+ 4f₃ + 2f₄+ 4f₅ + 2f₆+ …. + 4f_2M-1 + f_M)

Example:

_π/2

∫ sin(x) dx 2M = 2, M = 1

⁰

h= b - a = π/2-(0) = π/4

2M 2(1)

A = h (f₀ + 4f₁ + f₂) = (π/4) (1/3) [0 + 4 (1/2) √2 + 1] = 1.002279878

Now try 2M = 4 ; M =2

A = h (f_o + 4f₁ + 2f₂+ 4f₃ + f₄) = (π/8) [ sin (0) + 4 sin (π/8) + 2 sin (π/4) + 2 sin (3π/8) +

sin(π/2) = 1.000134585

Numerical Optimization

- Minimization of a function f(x). Let I = [a,b] an interval

- Local maximum value at x = p in the interval I = [a,b], if f(x) ≤ f(p) for all x in I

- Local Minimum value at x = p in the interval I = [a,b], if f(x) ≥ f(p) for all x in I

- A function is increasing in interval I = [a,b], if x₁ < x₂ and f(x₁) < f(x₂) for all x₁,x₂ in I

- If f’(x) > 0 for all x in I, then f(x) is increasing in the interval I

A function is decreasing in interval I = [a,b], if x₁ < x₂and f(x₁) > f(x₂) for all x₁,x₂ in I

If f’(x) <0 for all x in I, then f(x) is decreasing in the interval I

If there is a max or min value at x = p, then f’(p) = 0

If f’’(p) < 0, f(p) is local min

If f’’(x) > 0, f(p) is local max

Minimization using derivative

Assume we want to minimize f(x) and it has a unique min at pm a<p<b. we can solve the function f’(x) = 0 using any methods for non-linear equation – Newton, Secant, Bisection, Regula-False. f’(x) can be computed numerically, using the centered

The Golden Ratio

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

Then we evaluate f(x) at c and d

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

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

c & d 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.

We choose c, d such that c = a +(1-γ) (b-a)

d = a + γ (b-a)

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

If the new interval becomes [a,d] then the odd d will be the new b, and the old c becomes new d. The new c will have to be recomputed as well as f(c).

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:

f(x) = x² +1 [-2,1]

i =1; a = -2; b =1

c = a +(1-γ) (b-a) d = a + γ (b-a)

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

= -.8541 = -.14380

f(a₁) = 5 f(b₁) = f(1) = 2

f(c₁) = f(-.8541) = 1.72948

f(d₁) = f(-.14589) = 1.02188

f(c) > f(d)

Now interval is [-.8541, 1]

a = -.8541; b = 1; c = a + (1-γ) (b-a) = -.8341 + (1-.61803) (1-(-.8541)) = -.14164

d = a + γ (b-a) = -.8541 + .61803 (1-(-.8541)) = .291789

f(a₂) = 1.72949; f(b₂) = f(1) = 2; f(c₂) = f(-.14589) = 1.02188

f(d₂) = f(.291789) = 1.08514

Since f(c) < f(d); the minimum is the left side of d.

Now, a₃ = a₂ and b₃ = d₂

Look in the interval [a₂,d₂]

c = 3

a₃ = a₂ = -.8541; b₃ = d₂ = .291789; c₂ = a + (1-γ) (b-a) = -.4164

d₃ = c₂ = -.14584

f(a₃) = f(a₂) = f(-.8541) = 1.72949

f(b₃) = f(d₂) = f(.291789) 1.08514

f(c₃) = f(-.4164) = 1.17339

Since f(c₃) < f(d₃); a4 = c₃and b₄ = b₃. Now interval [-.1464, .291]

The Gradient Method or Steepset Desant

- Assume we want to minimize f(x) of N variables where x = (x₁, x₂, x₃, …, x_n)

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

∆f(x) = ( df(x), df(x), …, df(x) )

dx₁ dx₂ dx_n

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

Then -∆f(x) will point to the direction of the greatest decrease of f(x)

To obtain the minimum using the gradient method, starting point p₀ we move along the line in the direction of greatest decrease, that is -∆f(x).

In the simplest form, p₁ = p₀ – G_o .h

Where h = a small constant < 1

G is the gradient at P₀

P_k+1 = P_k – G_k. h; G_k= -∆f(x_k)

Example: y = x² +1 and y’ = 2x, gradient = G = 2x

Let h = .1

Start at

x₀ = -1, G₀ = 2(-1) = -2

x₁ = x₀ - G₀.h; G₁ = 2(-.8) = -1.6

= -1 - (-2)(.1) = -.8

x₂ = x₁ - G₁. h = (-.8) - (-1.6)(.1) = -.64, G₂ = 2(-.64) = -1.281

x₃ = x₂ - G₂. h = (-.64) - (-1.281)(.1) = -.512, G₃ = 2(-.512) = -

x₄ = x₃ - G₃. h = -.4096

x₅ = x₄ - G₄. h = -.32768

Numerical Method for Diff Equations

Consider the diff equation

dy = ky

dy = k.dt

∫ dy = ∫ k.dt

ln(y) = k₁t + k₂

Some diff equation does not have an analytical solution, so they have to be approximated using numerical methods.

1^st method: Euler’s Method

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