Summer 20-21

Numerical Methods for Science and Engineering

Lecture Note 7
Numerical Differentiation

7.1 Introduction
Numerical differentiation is the process of finding derivatives numerically for a function
whose values are given in data form generated from an experiment. For evenly distributed
data points and if we need the derivative at data points we may use the derivative formulas
called finite differences. When the data points are not even or required derivatives are at
points other than data points, we may use interpolating polynomials.
Numerical differentiation formulas can be derived by using the Taylor series expansion or by
differentiating the interpolating polynomials. Here we shall consider both way of deriving the
derivative formulas.

7.2 Method of Numerical Differentiation

Recall the definition of the derivative of a function
f ( x + Δx )−f ( x )
f ' ( x )= lim
Δx →0 Δx
For finite Δx ,
' f ( x + Δx )−f (x )
f ( x )≈
Δx
To find the derivative at
x=x
i , we choose another point x i+1 =xi +h ahead of x . This
i
gives two point forward difference formula
' f ( x i +h )−f ( x i ) '
f ( x i )≈ =f ( x i , h )
h (1)
If Δx is chosen as a negative number, say Δx=−h (h> 0), we have
f ( x i−h)−f ( x i ) f ( x i )−f ( x i−h) '
f ' ( x i )≈ ≈ =f (x i , h )
−h h (2)
This is backward difference formula for first derivative.
Adding Eq.(1) and Eq.(2), we have
' f ( x i +h )−f ( x i−h ) '
f ( x i )≈ =f ( x i , h)
2h (3)
which is a-point central difference formula for first derivative.

7.3 Derivative Formula from Taylor Series

For clear idea about the different formulas and their order of errors we may use the Taylor
series expansion of f (x ) .
From Taylor series expansion for h > 0, we have

Lec 6 - 1
 Summer 20-21

' h2 '' h3 '' '

f (x 0 +h )=f (x 0 )+h f (x 0 )+ f ( x 0 )+ f (x 0 )+ ¿ ⋅⋅¿ ¿
2! 3! (1)
2 3
' h '' h '' '
f (x 0 −h )=f ( x 0 )−h f ( x 0 )+ f ( x 0 )− f ( x 0 )+ ¿ ⋅⋅¿ ¿
2! 3! (2)
From the expansion of f (x 0 +h ) , we have
' f ( x 0 +h )−f (x 0 ) h '' h2 ' ''
f ( x 0 )= − f ( x 0 )− f ( x 0 )− ¿ ⋅⋅¿ ¿
h 2! 3!
'
which leads to the two-pont forward difference formula for f ( x 0 ) as
' f ( x 0 +h )−f ( x 0 )
f ( x 0 )= +E
h
where the error series is

E=−¿ ¿
From the expansion of f (x 0 −h ) , we have 2-point backward difference formula
' f ( x 0 )−f ( x 0 −h )
f ( x 0 )= +E
h
2
h '' h '''
E= f (x 0 )− f ( x 0 )+ ¿ ⋅⋅¿ ¿
with error term 2! 3!
In the two point formula the error series is of the form
E=a1 h+a 2 h 2 +a 3 h3 +⋅⋅⋅¿ ¿
where a’s does not depend on h.

By subtraction, we obtain
2 3 ' '' 2
f (x 0 +h )−f (x 0 −h )=2 hf ( x 0 )+ h f ( x 0 )+ h5 f ( v)( x 0 )+⋅⋅⋅¿ ¿
3! 5!
'
This leads to the 3-point central formula for approximating
f ( x0 )
' f ( x 0 +h )−f ( x 0 −h )
f ( x 0 )= +E
2h

With
E=−¿ ¿
Adding the Taylor series for f (x 0 +h ) and f (x 0 −h ) , we get
2 4 (4)
f (x 0 +h )+f ( x 0 −h )=2 f ( x 0 )+h2 f ( x 0 )+
''
h f ( x 0 )+⋅⋅⋅¿ ¿
4!
''
When this is rearranged, we get 3-point central difference formula for f ( x 0 )
f ( x 0 + h)−2 f ( x 0 )+ f ( x 0 −h )
f ' ' ( x 0 )= 2
+E
h
1 2 (4 ) 1
where the error series is
E=−2
[ 4!
h f ( x 0 )+ h 4 f (6) ( x 0 )⋅⋅¿
6! ]
In the three point central difference formula the error series is of the form

Lec 6 - 2
 Summer 20-21

E=a2 h2 +a 4 h 4 +a 6 h6+⋅⋅⋅¿ ¿
7.4 Formulas for Computing Derivatives
First Derivatives
f 1−f 0
f ' ( x 0 )≈ Ο(h)
h , 2-points forward difference
f 0 −f −1
f ' ( x 0 )≈ Ο(h)
h , 2-points backward difference
f 1−f −1
f ' ( x 0 )≈
2h , Ο(h 2 ) 3-points central difference
' 1
f ( x 0 )≈ [−3 f 0 + 4 f 1−f 2 ]
2h , Ο(h 2 ) 3-points forward difference
1
f ' ( x 0 )≈ [ 3 f 0 −4 f −1 +f −2 ] Ο(h 2 )
2h , 3-points backward difference
1
f ' ( x 0 )≈ [−f 2+8 f 1 −8 f −1 +f −2]
12 h , Ο(h 4 ) 5-points central difference

Second Derivatives
1
f ' ' ( x 0 )≈ [ f −1 −2 f 0 + f 1 ] 2
h2 , Ο(h ) 3-point central difference
'' 1
f ( x 0 )≈ 2 [ f 0 −2 f 1 + f 2 ]
h , Ο(h) 3-point forward difference
1
f ' ' ( x 0 )≈ 2 [ f 0 −2 f −1 + f −2 ]
h , Ο(h) 3-point backward difference
'' 1
f ( x 0 )≈ [−f 2+16 f 1−30 f 0 +16 f −1−f −2 ]
12 h2 , Ο(h 2 )
5-point central difference

7.5 Richardson Extrapolation

n
If the two approximations of order Ο(h ) for M are M (h1 ) and M (h2 ) , then the
Richardson’s extrapolated estimate MR of M can be written as
n
M R=M (h 1 )+ A (h1 ) (1)
n
M R =M (h 2 )+ A (h2 ) (2)
where it is assumed that the constant multiplicative factor A is same for both cases.
Subtracting (1) from (2),
0=M( h2 )−M (h 1 )+ A(h n2−hn1 )
M (h1 )−M (h2 )
A= n n
or h2 −h 1
Substituting in (1), we have

Lec 6 - 3
 Summer 20-21

n
h1 [ M (h1 )−M (h2 ) ]
M R =M (h 1 )+
h n2 −hn1 .
which can be written as
M (h1 )−M (h2 )
M R =M (h 1 )+
(h2 /h1 )n−1
or
M ( h1 )−M ( h2 )
M R =M ( h 1 )+ n
r −1
h2
r=
where h1 .
This is known as the Richardson extrapolation formula.

Lower order formula and Richardson extrapolation can be used to deduce the higher order
'
formula. For convenience we have used the notation f ( x 0 , h ) to indicate clearly the
'
approximation of f ( x 0 ) with step size h and f (x 0 +rh)=f r .

Thus the 3-point central difference formula for first derivative will ne written as
f −f
f ' ( x 0 , h )= 1 −1
2h

Example 7.1 The values of distance at various times are given below
Time (t) 4 6 8 10 12
Distance(s) 7.38 12.07 18.37 26.42 36.40
2

v = ds
a= d s2
The speed and acceleration can be calculated by dt and acceleration
. dt
(a) Using three point central difference formula estimate the speeds at (i) t=8, (ii) t=7,
and (iii) t=9.
(b) Using two point formulas and extrapolation estimate the speeds at (i) t=4, and (ii)
t=12.
(c) Use three points central difference formula and extrapolation to estimate speed at t=8
.
(d) Use three points central or forward or backward formula to estimate the accelerations
at (i) t=8, (ii) t=4, and (iii) t=12.
(e) Write down MATLAB code to estimate the speed and acceleration at time t=8 using
three point central difference formulas.
(f) Use MATAB functions “sp=spline(x,y)”, “fnder(sp, dorder)” and “fnval(sp, xo)”
to estimate the speed and acceleration at time t=6.5∧10.4.

Solution
(a) Three point central derivative formula for first derivative is
1
f ' ( x 0 , h )= f ( x 0+ h )−f ( x 0−h ) ]
2h [
(i) Speed at t=8

Lec 6 - 4
 Summer 20-21

1 1
v ( 8 , 2 )= [ s ( 10 )−s ( 6 ) ] = [26.42−12.07]=3.5875 .
2× 2 2
8−6
(ii) Speed at t=7. Here h= =1.
2
1
v ( 7 , 1 )= [ s ( 8 )−s ( 6 ) ]= 1 [18.37−12.07 ]=3.15 .
2× 1 4
10−8
(iii) Speed at t=9. Here h= =1.
2
1
v ( 9 , 1 )= [ s ( 10 )−s ( 8 ) ] = 1 [26.42−18.37]=4.025 .
2× 1 2
(b)
(i) For t=4, we have to use forward difference formula
1
f ' (x 0)≈ f ' ( x 0 , h )= [f ( x 0+ h )−f ( x 0 ) ].
h
Speed at t=4 we need to use h=2∧h=4.
1 1
v ( 4 , 2 )= [ s (6 )−s ( 4 ) ]= [ 12.07−7.38 ] =2.345 .
2 2
1
v ( 4 , 4 )= [ s ( 8 )−s ( 4 ) ] = 14 [ 18.37−7.38 ] =2.7475 .
4
Extrapolated value is
v ( 8,2 )−v (8,4)
v R ( 4 )=v ( 8,2 ) + 1
=2.345+ ( 2.345−2.7475 )=2.9425 .
2 −1
(ii) For t=12, we have to use backward difference formula
1
f ' (x 0)≈ f ' ( x 0 , h )= f ( x 0 )−f ( x0 −h ) ] .
h[
Speed at t=12 we need to use h=2∧h=4.
1 1
v ( 12 ,2 )= [ s (12)−s (10 ) ]= [ 36.40−26.42 ] =4.99.
2 2
1 1
v ( 12 , 4 )= [ s ( 12 ) −s ( 8 ) ]= [ 36.42−18.37 ] =4.5075 .
4 4
Extrapolated value is
v ( 12,2 )−v (12,4)
v R ( 12 )=v ( 12,2 ) + =4.99+ ( 4.99−4.5075 ) =5.4725 .
21−1

(c) Three points central difference formula is

1
f (x 0)≈ f ' ( x 0 , h ) = f ( x +h )−f ( x0 −h ) ] .
2h [ 0
Speed at t=8 we need to use h=2∧h=4.

Lec 6 - 5
 Summer 20-21

1 1
v ( 8 , 2 )= [ s ( 10 )−s ( 6 ) ] = [ 26.42−12.07 ] =3.5875 .
2× 2 4
1
v ( 8 , 4 )= [ s ( 12 )−s ( 4 ) ]= 1 [ 36.40−7.38 ] =3.6275 .
2 ×4 8
Extrapolated value is
v ( 8 , 2 )−v (8 , 4) 3.5875−3.6275
v R ( 8 )=v ( 8 ,2 ) + =3.5875+ =3.5742 .
2
2 −1 3

(d)
(i) For t=8, we have to use central difference formula (it gives better approximation).
Three point central derivative formula for second derivative is
1
f ' ' (x 0 ,h)= ( f ( x 0 +h ) −2 f ( x 0 ) + f ( x 0−h ) ) .
h2
Acceleration at t=8 is
1
a ( 8 , 2 )= 2
[ s ( 10 )−2 s ( 8 ) +s (6)]= 14 [26.42−2 ( 18.37 )+12.07 ]=0.4375.
2
(ii) For ¿ 4, three point forward difference formula for second derivative is
1
f ' ' (x 0 ,h)= ( f ( x 0 ) −2 f ( x 0+ h ) + f ( x 0 +2 h ) ) .
h2
Acceleration at t=4 is
1
a ( 4 , 2 )= 2
[ s ( 4 ) −2 s ( 6 )+ s (8)]= 14 [7.38−2 ( 12.07 ) +18.37 ]=0.4025 .
2
(iii) For ¿ 12, three point backward difference formula for second derivative is
1
f ' ' (x 0 ,h)= ( f ( x 0 ) −2 f ( x 0−h ) + f ( x 0−2 h ) ) .
h2
Acceleration at t=12 is
1
a ( 12 , 2 )= 2
[ s ( 12 )−2 s ( 10 ) +s (8)]= 14 [36.40−2 ( 26.42 )+18.37 ]=0.4825 .
2
(d)
>> clear
>> x=[6 8 10];
>> y=[12.07 18.37 26.42];
>> h=x(2)-x(1);
>> D1=(y(3)-y(1))/(2*h);
>> D2=(y(3)-2*y(2)+y(1))/h^2;

Lec 6 - 6
 Summer 20-21

(e) >> clear

>> x=[4 6 8 10 12];
>> y=[7.38 12.07 18.37 26.42 36.40];
>> % syntax for derivative is “fnder(f, dorder)”
>> sp=spline(x,y); % generates spline function sp
>> D1sp=fnder(sp,1); % generate first derivative of spline function sp
>> ValD1=fnval(D1sp,[8.4, 11]) % gives values from D1sp

ValD1 =

3.7501 4.9860

>> D2sp=fnder(sp, 2); % gererates second derivative

>> ValD2=fnval(D2sp, [8.4, 11]) % gives values from D2sp

ValD2 =

0.4445 0.5063

Lec 6 - 7
 Summer 20-21

Exercise 7
Numerical Differentiation
1. The distance s of a runner from a fixed point is measured (in metres) at intervals of half a
second. The data obtained is
Time t (s) 0 0.5 1 1.5 2
Distance s (m) 0 3.65 6.8 9.9 12.15
(a) Use two points difference formulas to approximate the runner’s speed at times t = 0s
and t = 2s.
(b) Use three points central difference formula to approximate the runner’s speed at times
t = 0.5s and t = 1.25s.
(c) Use three points central difference formula to approximate the runner’s acceleration at
times t = 1 s.
(d) Write down MATLAB code to estimate the speed and acceleration at time at time
t=1.5 s using three point central difference formulas.

3. The speed v (in m/s) of a rocket measured at half second intervals is

Time t (s) 0 0.5 1 1.5 2
speed v (in m/s) 0 11.860 26.33 41.075 59.05
5
(a) Use central difference formula to approximate the acceleration of the rocket at times t
= 1 s and t = 1.75s.
(b) Use two-point backward difference formula and Richardson extrapolation to estimate
the acceleration of the rocket at time t = 2 s.
(c) Use three-point central difference formula and extrapolation to estimate the
acceleration of the rocket at time t = 1 s.
(d) Use MATAB to estimate the acceleration of the rocket at time t=0.5 , 1.25∧2 using
spline interpolation.

4. The voltage E=E (t ) in an electric circuit obeys the differential equations

di
E(t )=L +Ri
dt , where R is the resistance and L is the inductance. Use L =0.05,
R = 2 and the value of I(t ) in the table’

t 1.0 1.1 1.2 1.3 1.4

i(t) 8.2277 7.2428 5.9908 4.5260 2.9122

di
(1.2)
(a) Find dt using three point central difference formula and extrapolation to
compute E(1.2)
(b) Compare your result with the exact solution I (t )=10 exp(−t/10 )sin 2 t .
(c) Use MATAB to estimate I ' (t) at each value of t using spline interpolation.

Lec 6 - 8
 Summer 20-21

Write down the MATLAB commands to find the corresponding voltage E(t ).

5. The distance traveled by an object is given in the table below:

t (s) 8 9 10 11 12
s(t) (m) 17.453 21.460 25.752 30.302 35.084

ds d2 s
The speed and acceleration can be calculated by v= and acceleration a= 2 .
dt dt
(a) Using three point central difference formula estimate the speeds at (i) t=9, and (ii) 
t=10.5.
(b) Using two point formulas estimate the speeds at (i) t=8, and (ii) t=12.
(c) Use three points central difference formula and extrapolation to estimate speed at
t=10.
(d) Use three points central or forward or backward difference formula to estimate the
accelerations at (i) t=10, (ii) t=8 and (iii) t=12.
(e) Use MATAB to estimate the speed and acceleration at time t¿ 8.5 , 10.5 ,∧11.2 using
spline interpolation.

6. The table below shows the values of f (x) at different values of x:

x 0.8 1.0 1.2 1.4 1.6
f(x) 0.954 1.648 2.623 3.947 5.697
(a) Using three point central difference formula estimate f ' ( 1 )∧f '(1.3)
(b) Using two point forward difference formula and extrapolation estimate f ' (0.8).
(c) Use three points central difference formula and extrapolation to estimate f ' ( 1.2 ).
(d) Use three points backward formula to estimate f ' ' (1.6).
(e) Write down MATAB codes using “sp=spline(x,y)”, “fnder(sp, dorder)” and
“fnval(sp, xo)” to estimate the values of f ' ( x )∧f ' ' (x) at
x=0.9 , 1.1∧1.42 .

7. A rod is rotating in a plane. The follow table gives the angle θ( in radians) through which
the rod has turned for various values of time t.
Time t 0 0.2 0.4 0.6 0.8
(s)
Angle θ 0 0.12 0.49 1.12 2.02
(a) Use two points difference formula to approximate the angular velocity of the rod at t
= 0.6s.
(b) Use three points central difference formula to approximate the angular velocity of the
rod at t = 0.6s.
(c) Use three points central difference formula to approximate the angular acceleration of
the rod at t = 0.6s.
(d) Write down MATLAB code to estimate the speed and acceleration at time t=0.5 s
using three point central difference formulas.

Lec 6 - 9
 Summer 20-21

8. The table below gives the results of an observation, θ is the observed temperature in
degree celsius of a vessel of cooling water, t is the time in minutes from the beginning of
observation.
Time t (min) 1 3 5 7 9
temperature θ 85.3 74.5 67. 60.5 54.3
0
(a) Use two points difference formula to approximate the rate of cooling at t=3 and t =5
(b) Use three points central difference formula to approximate the rate of cooling at t=3
and t=4
(c) Write down MATLAB code to estimate the approximate rate of cooling at t=3 and
t=4 three point central difference formulas.

Lec 6 - 10