Spring 2024-2025
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 ways 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 = xi , we choose another point xi +1 = xi + h ahead of xi . This
gives two-point forward difference formula
f ( xi + h) − f ( xi )
f ( xi ) = f ( xi , h) (1)
h
If x is chosen as a negative number, say x = −h (h> 0), we have
f ( xi − h) − f ( xi ) f ( xi ) − f ( xi − h)
f ( xi ) = f ( xi , h) (2)
−h h
This is the backward difference formula for first derivative.
Adding Eq. (1) and Eq. (2), we have
f ( xi + h) − f ( xi − h)
f ( xi ) = f ( xi , h) (3)
2h
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
h2 h3
f ( x0 + h) = f ( x0 ) + hf ( x0 ) +
f ( x0 ) + f ( x0 ) + (1)
2! 3!
1
Spring 2024-2025
h2 h3
f ( x0 − h) = f ( x0 ) − hf ( x0 ) + f ( x0 ) − f ( x0 ) + (2)
2! 3!
From the expansion of f ( x0 + h) , we have
f ( x0 + h) − f ( x0 ) h h2
f ( x0 ) = − f ( x0 ) − f ( x0 ) −
h 2! 3!
which leads to the two-point forward difference formula for f ( x0 ) as
f ( x0 + h) − f ( x0 )
f ( x0 ) = +E
h
where the error series is
h h2
E = − f ( x0 ) + f ( x0 ) +
2! 3!
From the expansion of f ( x0 − h) , we have 2-point backward difference formula,
f ( x0 ) − f ( x0 − h)
f ( x0 ) = +E
h
h h2
with error term E = f ( x0 ) − f ( x0 ) +
2! 3!
In the two-point formula the error series is of the form
E = a1h + a 2 h 2 + a3 h3 +
where a’s does not depend on h.
By subtraction, we obtain,
2 3 2
f ( x0 + h) − f ( x0 − h) = 2hf ( x0 ) + h f ( x0 ) + h5 f (v ) ( x0 ) +
3! 5!
This leads to the 3-point central formula for approximating f ( x0 )
f ( x0 + h) − f ( x0 − h)
f ( x0 ) = +E
2h
1 1
With E = − h2 f ( x0 ) + h4 f (v) ( x0 ) +
3! 5!
Adding the Taylor series for f ( x0 + h) and f ( x0 − h) , we get
2
f ( x0 + h) + f ( x0 − h) = 2 f ( x0 ) + h 2 f ( x0 ) + h 4 f ( 4) ( x0 ) +
4!
When this is rearranged, we get 3-point central difference formula for f ( x0 )
f ( x0 + h) − 2 f ( x0 ) + f ( x0 − h)
f ( x0 ) = +E
h2
1 1
where the error series is E = −2 h2 f (4) ( x0 ) + h4 f (6) ( x0 )
4! 6!
In the three-point central difference formula the error series is of the form
E = a2h 2 + a 4 h 4 + a 6 h6 +
2
Spring 2024-2025
7.4 Formulas for Computing Derivatives
First Derivatives
f1 − f 0
f ( x0 ) , (h) 2-points forward difference
h
f − f −1
f ( x0 ) 0 , (h) 2-points backward difference
h
f − f −1
f ( x0 ) 1 , (h 2 ) 3-points central difference
2h
f ( x0 )
1
− 3 f 0 + 4 f1 − f 2 , (h 2 ) 3-points forward difference
2h
f ( x0 )
1
3 f 0 − 4 f −1 + f −2 , (h 2 ) 3-points backward difference
2h
f ( x0 )
1
− f 2 + 8 f1 − 8 f −1 + f −2 , (h4 ) 5-points central difference
12h
Second Derivatives
f ( x0 )
1
f −1 − 2 f 0 + f1 , (h 2 ) 3-point central difference
h2
f ( x 0 )
1
f 0 − 2 f1 + f 2 , (h) 3-point forward difference
h2
f ( x 0 )
1
f 0 − 2 f −1 + f −2 , (h) 3-point backward difference
h2
f ( x0 )
1
− f 2 + 16 f1 − 30 f 0 + 16 f −1 − f −2 , (h 2 )
12h 2
5-point central difference
7.5 Richardson Extrapolation
If the two approximations of order (h n ) for M are M (h1) and M (h2 ) , then the
Richardson’s extrapolated estimate M R of M can be written as
M R = M (h1) + A(h1)n (1)
M R = M (h2 ) + A(h2 ) n (2)
where it is assumed that the constant multiplicative factor A is same for both cases.
Subtracting (1) from (2),
0 = M (h2 ) − M (h1) + A(h2n − h1n )
M (h1 ) − M (h2 )
or A=
h2n − h1n
Substituting in (1), we have
h n M (h1 ) − M (h2 )
M R = M (h1 ) + 1 .
h2n − h1n
which can be written as
3
Spring 2024-2025
M (h1) − M (h2 )
M R = M (h1 ) +
(h2 / h1 ) n − 1
or
M (h1) − M (h2 )
M R = M (h1 ) +
r n −1
h2
where r= .
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 ( x0 , h) to clearly indicate the
approximation of f ( x0 ) with step size h and f ( x0 + rh) = f r .
Thus the 3-point central difference formula for first derivative will be written as
f −f
f ( x0 , h) = 1 −1
2h
Example 7.1 The values of distance at various times are given below
Time (𝑡) 4 6 8 10 12
Distance(𝑠) 7.38 12.07 18.37 26.42 36.40
𝑑𝑠 𝑑2 𝑠
The speed and acceleration can be calculated by 𝑣 = 𝑑𝑡 and acceleration 𝑎 = 𝑑𝑡 2 .
(a) Using three-point central difference formula estimate the speeds at (i) 𝑡 = 8,
(ii) 𝑡 = 7, and (iii) 𝑡 = 9.
(b) Using two-point formulas and extrapolation estimate the speeds at (i) 𝑡 = 4, and
(ii) 𝑡 = 12.
(c) Use three points central difference formula and extrapolation to estimate speed at
𝑡 = 8.
(d) Use three points central or forward or backward formula to estimate the
accelerations at (i) 𝑡 = 8, (ii) 𝑡 = 4, and (iii) 𝑡 = 12.
(e) Write down MATLAB code to estimate the speed and acceleration at time 𝑡 = 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 𝑡 = 6.5 and 10.4.
Solution
(a) Three-point central derivative formula for first derivative is
1
𝑓 ′ (𝑥0 , ℎ) = [𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 − ℎ)]
2ℎ
(i) Speed at 𝑡 = 8
1 1
𝑣(8, 2) = [𝑠(10) − 𝑠(6)] = [26.42 − 12.07] = 3.5875.
2×2 2
8−6
(ii) Speed at 𝑡 = 7. Here ℎ = = 1.
2
4
Spring 2024-2025
1 1
𝑣(7, 1) = [𝑠(8) − 𝑠(6)] = [18.37 − 12.07] = 3.15.
2×1 4
10−8
(iii) Speed at 𝑡 = 9. Here ℎ = = 1.
2
1 1
𝑣(9, 1) = [𝑠(10) − 𝑠(8)] = [26.42 − 18.37] = 4.025.
2×1 2
(b)
(i) For 𝑡 = 4, we must use forward difference formula
1
𝑓′(𝑥0 ) ≈ 𝑓 ′ (𝑥0 , ℎ) = [𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 )].
ℎ
Speed at 𝑡 = 4 we need to use ℎ = 2 and ℎ = 4.
1 1
𝑣(4, 2) = [𝑠(6) − 𝑠(4)] = [12.07 − 7.38] = 2.345.
2 2
1 1
𝑣(4, 4) = [𝑠(8) − 𝑠(4)] = [18.37 − 7.38] = 2.7475.
4 4
Extrapolated value is
𝑣(4,2) − 𝑣(4,4)
𝑣𝑅 (4) = 𝑣(4,2) + = 2.345 + (2.345 − 2.7475) = 2.9425.
21 − 1
(ii) For 𝑡 = 12, we have to use backward difference formula
1
𝑓′(𝑥0 ) ≈ 𝑓 ′ (𝑥0 , ℎ) = [𝑓(𝑥0 ) − 𝑓(𝑥0 − ℎ)].
ℎ
Speed at 𝑡 = 12 we need to use ℎ = 2 and ℎ = 4.
1 1
𝑣(12, 2) = [𝑠(12) − 𝑠(10)] = [36.40 − 26.42] = 4.99.
2 2
1 1
𝑣(12, 4) = [𝑠(12) − 𝑠(8)] = [36.42 − 18.37] = 4.5075.
4 4
Extrapolated value is
𝑣(12,2) − 𝑣(12,4)
𝑣𝑅 (12) = 𝑣(12,2) + = 4.99 + (4.99 − 4.5075) = 5.4725.
21 − 1
(c) Three points central difference formula is
1
𝑓′(𝑥0 ) ≈ 𝑓 ′ (𝑥0 , ℎ) = [𝑓(𝑥0 + ℎ) − 𝑓(𝑥0 − ℎ)].
2ℎ
Speed at 𝑡 = 8 we need to use ℎ = 2 and ℎ = 4.
1 1
𝑣(8, 2) = [𝑠(10) − 𝑠(6)] = [26.42 − 12.07] = 3.5875.
2×2 4
1 1
𝑣(8, 4) = [𝑠(12) − 𝑠(4)] = [36.40 − 7.38] = 3.6275.
2×4 8
Extrapolated value is
5
Spring 2024-2025
𝑣(8, 2) − 𝑣(8, 4) 3.5875 − 3.6275
𝑣𝑅 (8) = 𝑣(8, 2) + = 3.5875 + = 3.5742.
22 − 1 3
(d)
(i) For 𝑡 = 8, we have to use central difference formula (it gives better approximation).
Three-point central derivative formula for second derivative is
1
𝑓′′(𝑥0 , ℎ) = (𝑓(𝑥0 + ℎ) − 2𝑓(𝑥0 ) + 𝑓(𝑥0 − ℎ)).
ℎ2
Acceleration at 𝑡 = 8 is
1 1
𝑎(8, 2) = [𝑠(10) − 2 𝑠(8) + 𝑠(6)] = [26.42 − 2(18.37) + 12.07] = 0.4375.
22 4
(ii) For = 4, three-point forward difference formula for second derivative is
1
𝑓′′(𝑥0 , ℎ) = (𝑓(𝑥0 ) − 2𝑓(𝑥0 + ℎ) + 𝑓(𝑥0 + 2ℎ)).
ℎ2
Acceleration at 𝑡 = 4 is
1 1
𝑎(4, 2) = 2
[𝑠(4) − 2 𝑠(6) + 𝑠(8)] = [7.38 − 2(12.07) + 18.37] = 0.4025.
2 4
(iii) For = 12, three-point backward difference formula for second derivative is
1
𝑓′′(𝑥0 , ℎ) = (𝑓(𝑥0 ) − 2𝑓(𝑥0 − ℎ) + 𝑓(𝑥0 − 2ℎ)).
ℎ2
Acceleration at 𝑡 = 12 is
1 1
𝑎(12, 2) = [𝑠(12) − 2 𝑠(10) + 𝑠(8)] = [36.40 − 2(26.42) + 18.37] = 0.4825.
22 4
(e)
>> 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;
(f) >> 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,[6.5, 10.4]) % gives values from D1sp
6
Spring 2024-2025
ValD1 =
3.7501 4.9860
>> D2sp=fnder(sp, 2); % gererates second derivative
>> ValD2=fnval(D2sp, [6.5, 10.4]) % gives values from D2sp
ValD2 =
0.4445 0.5063
7
Spring 2024-2025
Exercise 7
Numerical Differentiation
1. The distance s of a runner from a fixed point is measured (in meters) 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
𝑡 = 1.5 𝑠 using three-point central difference formulas.
2. 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.335 41.075 59.05
(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 MATLAB to estimate the acceleration of the rocket at time 𝑡 =
0.5, 1.25 and 2 using spline interpolation.
3. The voltage E = E (t ) in an electric circuit obeys the differential equations
di
E (t ) = L + Ri , where R is the resistance and L is the inductance. Use L =0.05,
dt
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
(a) Find (1.2) using three point central difference formula and extrapolation to
dt
compute E(1.2)
(b) Compare your result with the exact solution I (t ) = 10exp( −t / 10) sin 2t .
(c) Use MATLAB to estimate 𝐼′(𝑡) at each value of t using spline interpolation.
Write down the MATLAB commands to find the corresponding voltage 𝐸(𝑡).
8
Spring 2024-2025
4. 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
𝑑𝑠 𝑑2 𝑠
The speed and acceleration can be calculated by 𝑣 = 𝑑𝑡 and acceleration 𝑎 = 𝑑𝑡 2 .
(a) Using three-point central difference formula estimate the speeds at (i) 𝑡 = 9, and
(ii) 𝑡 = 10.5.
(b) Using two-point formulas estimate the speeds at (i) 𝑡 = 8, and (ii) 𝑡 = 12.
(c) Use three points central difference formula and extrapolation to estimate speed at
𝑡 = 10.
(d) Use three points central or forward or backward difference formula to estimate the
accelerations at (i) 𝑡 = 10, (ii) 𝑡 = 8 and (iii) 𝑡 = 12.
(e) Use MATLAB to estimate the speed and acceleration at time t= 8.5,
10.5, and 11.2 using spline interpolation.
5. The table below shows the values of 𝑓(𝑥) 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 𝑓 ′ (1) and 𝑓′(1.3)
(b) Using two-point forward difference formula and extrapolation estimate 𝑓′(0.8).
(c) Use three points central difference formula and extrapolation to estimate 𝑓 ′ (1.2).
(d) Use three points backward formula to estimate 𝑓′′(1.6).
(e) Write down MATAB codes using “sp=spline(x,y)”, “fnder(sp, dorder)” and
“fnval(sp, xo)” to estimate the values of 𝑓 ′ (𝑥) and 𝑓′′(𝑥) at
𝑥 = 0.9, 1.1 and 1.42 .
6. 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 𝑡.
Time t (s) 0 0.2 0.4 0.6 0.8
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
𝑡 = 0.4 𝑠 using three-point central difference formulas.
9
Spring 2024-2025
7. The table below gives the results of an observation, 𝜃 is the observed temperature in
degree celsius of a vessel of cooling water, 𝑡 is the time in minutes from the
beginning of observation.
Time t (min) 1 3 5 7 9
temperature 𝜃 85.3 74.5 67.0 60.5 54.3
(a) Use two points difference formula to approximate the rate of cooling at 𝑡 = 3 and
𝑡=5
(b) Use three points central difference formula to approximate the rate of cooling at
𝑡 = 3 and 𝑡 = 4
(c) Write down MATLAB code to estimate the approximate rate of cooling at 𝑡 = 3
and 𝑡 = 4 three-point central difference formulas.
8. Consider the viscous flow of air over a flat plate. At a given station in the flow
direction, the velocity profile in the direction perpendicular to the plate is given by:
𝑦
𝑢 = 1582 (1 − 𝑒 −𝐿 ), where 𝐿 = 1 𝑓𝑡, and 𝑢 is in feet per second. Given that shear
𝑑𝑢
stress is 𝜏𝑤 = 𝜇 (𝑑𝑦 ), the dynamic viscosity of air is 𝜇 = 3.7373 × 10−5 𝑠𝑙𝑢𝑔/
(𝑓𝑡. 𝑠). The velocity measurements at different distances from the plate are given
in the table below
𝑦(𝑓𝑡) 0 0.20 0.40 0.60
𝑢(𝑓𝑡/𝑠) 0 286.77 521.55 713.78
(a) Find shear stress 𝜏𝑤 from the above data at 𝑦 = 0.20, using 3-points central
difference formula.
(b) Find the exact value and compare it with your result.
9. The following data show estimates of the population of Liberia in selected years
between 1960 to 2010
Year 1960 1970 1980 1990 2000 2010
Population (millions) 1.1 1.4 1.9 2.1 2.8 4
Calculate the rate of growth of the population in millions per year for 2010, using
two-point backward difference formula.
[Ref. Numerical Methods for engineers and Scientists – Amos Gilat, Vish
Subramaniam, Problem # 8.3 (a), Page # 330]
10. The following data is given for the stopping distance of a car on a wet road verses
the speed at which it begins barking.
𝑣 (mil/h) 12.5 25 37.5 50 62.5 75
𝑑 (ft) 20 59 118 197 299 420
Calculate the rate of change of the stopping distance at a speed of 62.5 mph using
the two-point backward difference formula.
[Ref. Numerical Methods for engineers and Scientists – Amos Gilat, Vish
Subramaniam, Problem # 8.4 (a)(i), Page # 330]
10