Numerical Differentiation CIS 541 - Differentiation The - - PowerPoint PPT Presentation

numerical differentiation cis 541 differentiation
SMART_READER_LITE
LIVE PREVIEW

Numerical Differentiation CIS 541 - Differentiation The - - PowerPoint PPT Presentation

Numerical Differentiation CIS 541 - Differentiation The mathematical definition: Roger Crawfis + f x ( h ) f x ( ) '( ) lim = f x h h 0 Can also be thought of as the tangent line. x x+h August 12, 2005 OSU/CIS


slide-1
SLIDE 1

CIS 541 - Differentiation

Roger Crawfis

August 12, 2005 OSU/CIS 541 2

Numerical Differentiation

  • The mathematical definition:
  • Can also be thought of as the tangent line.

( ) ( ) '( ) lim

h

f x h f x f x h

+ − =

x x+h

August 12, 2005 OSU/CIS 541 3

Numerical Differentiation

  • We can not calculate the limit as h goes to zero, so

we need to approximate it.

  • Apply directly for a non-zero h leads to the slope
  • f the secant curve.

x x+h

August 12, 2005 OSU/CIS 541 4

Numerical Differentiation

  • This is called Forward Differences and can

be derived using Taylor’s Series:

2 2

( ) ( ) ( ) ( ) 2! ( ) ( ) ( ) ( ) 2! ( ) ( ) ( ) ( ) 2! ( ) ( ) ( ) h f x h f x f x h f h f x h f x f x h f f x h f x h f x f h f x h f x f x as h h ξ ξ ξ ′ ′′ + = + + ′ ′′ ∴ + − = + + − ′ ′′ ∴ = − + − ′ ∴ → → Theoretically speaking

slide-2
SLIDE 2

August 12, 2005 OSU/CIS 541 5

Truncation Errors

  • Let f(x) = a+e, and f(x+h) = a+f.
  • Then, as h approaches zero, e<<a and f<<a.
  • With limited precision on our computer, our

representation of f(x) ≈ a ≈ f(x+h).

  • We can easily get a random round-off bit as

the most significant digit in the subtraction.

  • Dividing by h, leads to a very wrong answer

for f’(x).

August 12, 2005 OSU/CIS 541 6

Error Tradeoff

  • Using a smaller step size reduces truncation error.
  • However, it increases the round-off error.
  • Trade off/diminishing returns occurs: Always think and test!

Log error Log step size Truncation error Round off error Total error Point of diminishing returns

August 12, 2005 OSU/CIS 541 7

Numerical Differentiation

  • This formula favors (or biases towards) the

right-hand side of the curve.

  • Why not use the left?

x x+h x-h

August 12, 2005 OSU/CIS 541 8

Numerical Differentiation

  • This leads to the Backward Differences

formula.

2

( ) ( ) ( ) ( ) 2! ( ) ( ) ( ) ( ) 2! ( ) ( ) ( ) h f x h f x f x h f f x f x h h f x f h f x f x h f x as h h ξ ξ ′ ′′ − = − + − − ′ ′′ ∴ = + − − ′ ∴ → →

slide-3
SLIDE 3

August 12, 2005 OSU/CIS 541 9

Numerical Differentiation

  • Can we do better?
  • Let’s average the two:
  • This is called the Central Difference

formula.

1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 f x h f x f x f x h f x h f x h f x h h h + − − − + − − ⎛ ⎞ ′ ≈ + = ⎜ ⎟ ⎝ ⎠

Forward difference Backward difference

August 12, 2005 OSU/CIS 541 10

Central Differences

  • This formula does not seem very good.

– It does not follow the calculus formula. – It takes the slope of the secant with width 2h. – The actual point we are interested in is not even evaluated.

x x+h x-h

August 12, 2005 OSU/CIS 541 11

Numerical Differentiation

  • Is this any better?
  • Let’s use Taylor’s Series to examine the error:

2 3 2 3 3 3

( ) ( ) ( ) ( ) ( ) 2 3! ( ) ( ) ( ) ( ) ( ) 2 3! ( ) ( ) 2 ( ) ( ) ( ) 3! 3! h h f x h f x f x h f x f h h f x h f x f x h f x f subtracting h h f x h f x h f x h f f ξ ς ξ ς ′ ′′ ′′′ + = + + + ′ ′′ ′′′ − = − + − ⎛ ⎞ ′ ′′′ ′′′ + − − = + + ⎜ ⎟ ⎝ ⎠

August 12, 2005 OSU/CIS 541 12

Central Differences

  • The central differences formula has much

better convergence.

  • Approaches the derivative as h2 goes to

zero!!

[ ]

2

( ) ( ) 1 ( ) ( ) , , 2 6 f x h f x h f x f h x h x h h ζ ζ + − − ′ ′′′ = − ∈ − +

( )

2

( ) ( ) ( ) 2 f x h f x h f x O h h + − − ′ = +

slide-4
SLIDE 4

August 12, 2005 OSU/CIS 541 13

Warning

  • Still have truncation error problem.
  • Consider the case of:
  • Build a table with

smaller values of h.

  • What about large

values of h for this function?

( ) 100 100 100 ( ) 2 1, 0.000333, 6 0.0100033 0.0099966 ( ) 0.010050 0.000666666 Relative error: 0.01-0.010050 0.5% 0.01 x f x x h x h f x h at x h with significant digits f x = + − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ′ = = − ′ = = ฀ ฀

August 12, 2005 OSU/CIS 541 14

Richardson Extrapolation

  • Can we do better?
  • Is my choice of h a good one?
  • Let’s subtract the two Taylor Series expansions again:

( ) ( ) ( ) ( ) ( )

2 3 4 5 4 5 2 3 4 5 4 5 5 5 3 3

( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3! 4! 5! ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3! 4! 5! ( ) ( ) ( ) ( ) 2 ( ) 2 2 2 ( ) 3! 3! 5! h h h h f x h f x f x h f x f x f x f x h h h h f x h f x f x h f x f f x f x subtracting f x f x h f x h f x h f x h h h f x ς ′ ′′ ′′′ + = + + + + + + ′ ′′ ′′′ − = − + − + − + ′′′ ′′′ ′ + − − = + + + + L L L

August 12, 2005 OSU/CIS 541 15

Richardson Extrapolation

  • Assuming the higher derivatives exist, we can

hold x fixed (which also fixes the values of f(x)), to obtain the following formula.

  • Richardson Extrapolation examines the operator

below as a function of h.

[ ]

2 4 6 2 4 6

1 ( ) ( ) ( ) 2 f x f x h f x h a h a h a h h ′ = + − − + + + +L

[ ]

1 ( ) ( ) ( ) 2 h f x h f x h h ϕ = + − −

August 12, 2005 OSU/CIS 541 16

Richardson Extrapolation

  • This function approximates f’(x) to O(h2) as

we saw earlier.

  • Let’s look at the operator as h goes to zero.

2 4 6 2 4 6 2 4 6 2 4 6

( ) ( ) ( ) ( ) 2 2 2 2 h f x a h a h a h h h h h f x a a a ϕ ϕ ′ = − − − − ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ′ = − − − − ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ L L

Same leading constants

slide-5
SLIDE 5

August 12, 2005 OSU/CIS 541 17

Richardson Extrapolation

  • Using these two formula’s, we can come up

with another estimate for the derivative that cancels out the h2 terms.

( )

4 6 4 6 4

3 15 ( ) 4 ( ) 3 ( ) 2 4 16 1 ( ) ( ) ( ) ( ) 2 3 2 h h f x a h a h

  • r

h h f x h O h ϕ ϕ ϕ ϕ ϕ ′ − = − − − − ⎡ ⎤ ′ = + − + ⎢ ⎥ ⎣ ⎦ L

new estimate difference between old and new estimates

Extrapolates by assuming the new estimate undershot.

August 12, 2005 OSU/CIS 541 18

Richardson Extrapolation

  • If h is small (h<<1), then h4 goes to zero

much faster than h2.

  • Cool!!!
  • Can we cancel out the h6 term?

– Yes, by using h/4 to estimate the derivative.

August 12, 2005 OSU/CIS 541 19

Richardson Extrapolation

  • Consider the following property:
  • where L is unknown,
  • as are the coefficients, a2k.

2 2 1 2 2 1

( ) ( )

k k k k k k

h f x a h L a h ϕ

∞ = ∞ =

′ = − = −

∑ ∑

( )

lim ( )

h

L h f x ϕ

′ = =

August 12, 2005 OSU/CIS 541 20

Richardson Extrapolation

  • Do not forget the formal definition is simply

the central-differences formula:

  • New symbology (is this a word?):

[ ]

1 ( ) ( ) ( ) 2 h f x h f x h h ϕ = + − −

( ) ( )

2 1

,0 2 ,0 2

n k n k

h D n h L A k ϕ

∞ =

⎛ ⎞ ≡ ⎜ ⎟ ⎝ ⎠ ⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠

From previous slide

slide-6
SLIDE 6

August 12, 2005 OSU/CIS 541 21

Richardson Extrapolation

  • D(n,0) is just the central differences
  • perator for different values of h.
  • Okay, so we proceed by computing D(n,0)

for several values of n.

  • Recalling our cancellation of the h2 term.

( )

[ ]

( )

4 4

1 ( ) ( ) ( ) ( ) 2 3 2 1 (1,0) (1,0) (0,0) 4 1 h h f x h O h D D D O h ϕ ϕ ϕ ⎡ ⎤ ′ = + − + ⎢ ⎥ ⎣ ⎦ = + − + −

August 12, 2005 OSU/CIS 541 22

Richardson Extrapolation

  • If we let h→h/2, then in general, we can

write:

  • Let’s denote this operator as:

[ ]

4

1 ( ) ( ,0) ( ,0) ( 1,0) 4 1 2n h f x D n D n D n O ⎛ ⎞ ⎛ ⎞ ′ = + − − + ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ − ⎝ ⎠ ⎝ ⎠

( ) ( ) ( )

1

1 ,1 ( ,0) ,0 1,0 4 1 D n D n D n D n = + − − ⎡ ⎤ ⎣ ⎦ −

August 12, 2005 OSU/CIS 541 23

Richardson Extrapolation

  • Now, we can formally define Richardson’s

extrapolation operator as:

  • or

( ) ( ) ( )

4 1 , ( , 1) 1, 1 , 1 4 1 4 1

m m m

D n m D n m D n m m n = − − − − ≤ ≤ − −

new estimate

  • ld estimate

( ) ( ) ( )

1 , ( , 1) , 1 1, 1 4 1

m

D n m D n m D n m D n m = − + − − − − ⎡ ⎤ ⎣ ⎦ −

August 12, 2005 OSU/CIS 541 24

Richardson Extrapolation

  • Now, we can formally define Richardson’s

extrapolation operator as:

Memorize me!!!!

( ) ( ) ( )

1 , ( , 1) , 1 1, 1 4 1

m

D n m D n m D n m D n m = − + − − − − ⎡ ⎤ ⎣ ⎦ −

slide-7
SLIDE 7

August 12, 2005 OSU/CIS 541 25

Richardson Extrapolation Theorem

  • These terms approach f’(x) very quickly.

( ) ( )

2 1

, , 2

k n k m

h D n m L A k m

∞ = +

⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠

Order starts much higher!!!!

August 12, 2005 OSU/CIS 541 26

Richardson Extrapolation

  • Since m≤ n, this leads to a two-dimensional

triangular array of values as follows:

  • We must pick an initial value of h and a

max iteration value N.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

0,0 1,0 1,1 2,0 2,1 2,2 ,0 ,1 ,2 , D D D D D D D N D N D N D N N M M M O L

August 12, 2005 OSU/CIS 541 27

Example

( )

( ) ( ) ( ) ( ) ( ) ( )

5 2 3

cos(100 ( ) 1 1.3, 128 0,0 16.696386 1,0 40.583393 2,0 109.322528 3,0 135.031747 4,0 142.068615 5,0 143.866937 x f x x x h D D D D D D = = = = = = = = =

1 0.00244 4096 h = =

August 12, 2005 OSU/CIS 541 28

Example

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

0,0 16.696386 1,0 40.583393 1,1 48.583393 2,0 109.322528 2,1 132.235574 3,0 135.031747 3,1 143.601487 4,0 142.068615 4,1 144.414238 5,0 143.866937 5,1 144.466377 D D D D D D D D D D D = = = = = = = = = = =

1 3 extrapolate

slide-8
SLIDE 8

August 12, 2005 OSU/CIS 541 29

Example

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

0,0 16.696386 1,0 40.583393 1,1 48.583393 2,0 109.322528 2,1 132.235574 2,2 137.814897 3,0 135.031747 3,1 143.601487 3,2 144.359214 4,0 142.068615 4,1 144.414238 4,2 144.468421 5,0 143.866937 5,1 144 D D D D D D D D D D D D D D = = = = = = = = = = = = = =

( )

.466377 5,2 144.469853 D =

1 15 extrapolate

August 12, 2005 OSU/CIS 541 30

Example

  • Which converges up to eight decimal

places.

  • Is it accurate?

16.696386 40.583393 48.583393 109.322528 132.235574 137.814897 135.031747 143.601487 144.359214 144.463092 142.068615 144.414238 144.468421 144.470154 144.470182 143.866937 144.466377 144.469853 144.469876 144.469875 5, D(

)

5 144.469875 =

1 1023 extrapolate 1 255 extrapolate 1 63 extrapolate 1 15 extrapolate 1 3 extrapolate August 12, 2005 OSU/CIS 541 31

Example

  • We can look at the (theoretical) error term
  • n this example.
  • Taking the derivative:

( ) ( ) ( ) ( )

2 5 5 1 12 2 5 7

5,5 ,5 2 1 1.3 (6,5) ,5 4096 2

k k k k

h D L A k h f A A k

∞ = + ∞ =

⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ′ = + + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠

∑ ∑

2-144 (1.3) 144.469874253 f ′ = K

Round-off error

August 12, 2005 OSU/CIS 541 32

Second Derivatives

  • What if we need the second derivative?
  • Any guesses?

( ) ( ) ( ) ( )

2 3 4 5 4 5 2 3 4 5 4 5

( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3! 4! 5! ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3! 4! 5! h h h h f x h f x f x h f x f x f x f x h h h h f x h f x f x h f x f f x f x ς ′ ′′ ′′′ + = + + + + + + ′ ′′ ′′′ − = − + − + − + L L

slide-9
SLIDE 9

August 12, 2005 OSU/CIS 541 33

Second Derivatives

  • Let’s cancel out the odd derivatives and

double up the even ones:

– Implies adding the terms together.

( )

2 4 4

( ) ( ) 2 ( ) 2 ( ) 2 ( ) 2 4! h h f x h f x h f x f x f x ′′ + + − = + + +L

August 12, 2005 OSU/CIS 541 34

Second Derivatives

  • Isolating the second derivative term yields:
  • With an error term of:

2

( ) 2 ( ) ( ) ( ) f x h f x f x h f x h + − + − ′′ ≈

( ) ( )

4 2

1 12 E h f ξ = −

August 12, 2005 OSU/CIS 541 35

Partial Derivatives

  • Remember: Nothing special about partial

derivatives:

( ) ( ) ( ) ( ) ( ) ( )

, , , 2 , , , 2 f x h y f x h y f x y x h f x y h f x y h f x y y h + − − ∂ ≈ ∂ + − − ∂ ≈ ∂

August 12, 2005 OSU/CIS 541 36

Calculating the Gradient

  • For lab 2, you need to calculate the

gradient.

  • Just use central differences for each partial

derivative.

  • Remember to normalize it (divide by its

length).