Numerical Differentiation

5 Numerical Differentiation
Given the function f (x), compute d n f/dx n at given x
5.1 Introduction
Numerical differentiation deals with the following problem: We are given the func-
tion y = f (x) and wish to obtain one of its derivatives at the point x = xk . The term
“given” means that we either have an algorithm for computing the function, or we
possess a set of discrete data points (xi , yi ), i = 0, 1, . . . , n. In either case, we have ac-
cess to a finite number of (x, y) data pairs from which to compute the derivative. If
you suspect by now that numerical differentiation is related to interpolation, you are
correct—one means of finding the derivative is to approximate the function locally
by a polynomial and then to differentiate it. An equally effective tool is the Taylor se-
ries expansion of f (x) about the point xk , which has the advantage of providing us
with information about the error involved in the approximation.
Numerical differentiation is not a particularly accurate process. It suffers from
a conflict between roundoff errors (caused by limited machine precision) and errors
inherent in interpolation. For this reason, a derivative of a function can never be com-
puted with the same precision as the function itself.
5.2 Finite Difference Approximations
The derivation of the finite difference approximations for the derivatives of f (x) is
based on forward and backward Taylor series expansions of f (x) about x, such as
h2 h3 h4 (4)
f (x + h) = f (x) + hf (x) + f (x) + f (x) + f (x) + · · · (a)
2! 3! 4!
h2 h3 h4 (4)
f (x − h) = f (x) − hf (x) + f (x) − f (x) + f (x) − · · · (b)
2! 3! 4!
183
Downloaded from https:/www.cambridge.org/core. University of Florida, on 20 Jan 2017 at 22:01:31, subject to the Cambridge Core terms of use, available at
https:/www.cambridge.org/core/terms. https://doi.org/10.1017/CBO9781139523899.006
(2h)2 (2h)3 (2h)4 (4)

f (x + 2h) = f (x) + 2hf (x) + f (x) + f (x) + f (x) + · · · (c)
2! 3! 4!
(2h)2 (2h)3 (2h)4 (4)

f (x − 2h) = f (x) − 2hf (x) + f (x) − f (x) + f (x) − · · · (d)
2! 3! 4!
We also record the sums and differences of the series:
h4 (4)
f (x + h) + f (x − h) = 2f (x) + h2 f (x) + f (x) + · · · (e)
12
h3
f (x + h) − f (x − h) = 2hf (x) + f (x) + . . . (f)
3
4h4 (4)
f (x + 2h) + f (x − 2h) = 2f (x) + 4h2 f (x) + f (x) + · · · (g)
3
8h3
f (x + 2h) − f (x − 2h) = 4hf (x) + f (x) + · · · (h)
3
Note that the sums contain only even derivatives, whereas the differences retain just
the odd derivatives. Equations (a)–(h) can be viewed as simultaneous equations that
can be solved for various derivatives of f (x). The number of equations involved and
the number of terms kept in each equation depend on the order of the derivative and
the desired degree of accuracy.
First Central Difference Approximations

The solution of Eq. (f) for f (x) is
f (x + h) − f (x − h) h2
f (x) = − f (x) − · · ·
2h 6
or
f (x + h) − f (x − h)
f (x) = + O(h2 ) (5.1)
2h
which is called the first central difference approximation for f (x). The term O(h2 )
reminds us that the truncation error behaves as h2 .
Similarly, Eq. (e) yields the first central difference approximation for f (x):
f (x + h) − 2f (x) + f (x − h) h2 (4)
f (x) = 2
+ f (x) + . . .
h 12
or
f (x + h) − 2f (x) + f (x − h)
f (x) = + O(h2 ) (5.2)
h2
185 5.2 Finite Difference Approximations
Central difference approximations for other derivatives can be obtained from

Eqs. (a)–(h) in the same manner. For example, eliminating f (x) from Eqs. (f) and
(h) and solving for f (x) yields
f (x + 2h) − 2f (x + h) + 2f (x − h) − f (x − 2h)
f (x) = + O(h2 ) (5.3)
2h3
The approximation
f (x + 2h) − 4f (x + h) + 6f (x) − 4f (x − h) + f (x − 2h)
f (4) (x) = + O(h2 ) (5.4)
h4
is available from Eq. (e) and (g) after eliminating f (x). Table 5.1 summarizes the
results.
f (x − 2h) f (x − h) f (x) f (x + h) f (x + 2h)

2hf (x) −1 0 1
h2 f (x) 1 −2 1
2h3 f (x) −1 2 0 −2 1
h4 f (4) (x) 1 −4 6 −4 1
Table 5.1. Coefficients of Central Finite Difference

Approximations of O(h2 )
First Noncentral Finite Difference Approximations

Central finite difference approximations are not always usable. For example, con-
sider the situation where the function is given at the n discrete points x0 , x1 , . . . , xn .
Because central differences use values of the function on each side of x, we would
be unable to compute the derivatives at x0 and xn . Clearly, there is a need for finite
difference expressions that require evaluations of the function only on one side of x.
These expressions are called forward and backward finite difference approximations.
Noncentral finite differences can also be obtained from Eqs. (a)–(h). Solving
Eq. (a) for f (x) we get
f (x + h) − f (x) h h2 h3 (4)
f (x) = − f (x) − f (x) − f (x) − · · ·
h 2 6 4!
Keeping only the first term on the right-hand side leads to the first forward difference
approximation:
f (x + h) − f (x)
f (x) = + O(h) (5.5)
h
Similarly, Eq. (b) yields the first backward difference approximation:
f (x) − f (x − h)
f (x) = + O(h) (5.6)
h
Note that the truncation error is now O(h), which is not as good as O(h2 ) in central
difference approximations.
We can derive the approximations for higher derivatives in the same manner. For
example, Eqs. (a) and (c) yield
f (x + 2h) − 2f (x + h) + f (x)
f (x) = + O(h) (5.7)
h2
The third and fourth derivatives can be derived in a similar fashion. The results
are shown in Tables 5.2a and 5.2b.
f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h)

hf (x) −1 1
h2 f (x) 1 −2 1
h3 f (x) −1 3 −3 1
h4 f (4) (x) 1 −4 6 −4 1
Table 5.2a. Coefficients of Forward Finite Difference

Approximations of O(h)
f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

hf (x) −1 1
h2 f (x) 1 −2 1
h3 f (x) −1 3 −3 1
h4 f (4) (x) 1 −4 6 −4 1
Table 5.2b. Coefficients of Backward Finite Difference

Approximations of O(h)
Second Noncentral Finite Difference Approximations

Finite difference approximations of O(h) are not popular because of reasons that
are explained shortly. The common practice is to use expressions of O(h2 ). To ob-
tain noncentral difference formulas of this order, we have to retain more term in the
Taylor series. As an illustration, we derive the expression for f (x). We start with
Eqs. (a) and (c), which are
h2 h3 h4 (4)
f (x + h) = f (x) + hf (x) + f (x) + f (x) + f (x) + · · ·
2 6 24
4h3 2h4 (4)
f (x + 2h) = f (x) + 2hf (x) + 2h2 f (x) + f (x) + f (x) + · · ·
3 3
We eliminate f (x) by multiplying the first equation by 4 and subtracting it from the
second equation. The result is
h4 (4)
f (x + 2h) − 4f (x + h) = −3f (x) − 2hf (x) + f (x) + · · ·
2
187 5.2 Finite Difference Approximations
Therefore,
−f (x + 2h) + 4f (x + h) − 3f (x) h2 (4)
f (x) = + f (x) + · · ·
2h 4
or
−f (x + 2h) + 4f (x + h) − 3f (x)
f (x) = + O(h2 ) (5.8)
2h
Equation (5.8) is called the second forward finite difference approximation.
Derivation of finite difference approximations for higher derivatives involves ad-
ditional Taylor series. Thus the forward difference approximation for f (x) uses series
for f (x + h), f (x + 2h), and f (x + 3h); the approximation for f (x) involves Taylor ex-
pansions for f (x + h), f (x + 2h), f (x + 3h), f (x + 4h), and so on. As you can see, the
computations for high-order derivatives can become rather tedious. The results for
both the forward and backward finite differences are summarized in Tables 5.3a and
5.3b.
f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h) f (x + 5h)

2hf (x) −3 4 −1
h2 f (x) 2 −5 4 −1
2h3 f (x) −5 18 −24 14 −3
h4 f (4) (x) 3 −14 26 −24 11 −2
Table 5.3a. Coefficients of Forward Finite Difference Approximations of O(h2 )
f (x − 5h) f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

2hf (x) 1 −4 3
h2 f (x) −1 4 −5 2
2h3 f (x) 3 −14 24 −18 5
h4 f (4) (x) −2 11 −24 26 −14 3
Table 5.3b. Coefficients of Backward Finite Difference Approximations of O(h2 )
Errors in Finite Difference Approximations

Observe that in all finite difference expressions the sum of the coefficients is zero.
The effect on the roundoff error can be profound. If h is very small, the values of f (x),
f (x ± h), f (x ± 2h), and so on, will be approximately equal. When they are multiplied
by the coefficients and added, several significant figures can be lost. Yet we cannot
make h too large, because then the truncation error would become excessive. This
unfortunate situation has no remedy, but we can obtain some relief by taking the
following precautions:
• Use double-precision arithmetic.

• Employ finite difference formulas that are accurate to at least O(h2 ).
To illustrate the errors, let us compute the second derivative of f (x) = e−x at x =
1 from the central difference formula, Eq. (5.2). We carry out the calculations with
six- and eight-digit precision, using different values of h. The results, shown in Table
5.4, should be compared with f (1) = e−1 = 0.367 879 44.
h Six-Digit Precision Eight-Digit Precision

0.64 0.380 610 0.380 609 11
0.32 0.371 035 0.371 029 39
0.16 0.368 711 0.368 664 84
0.08 0.368 281 0.368 076 56
0.04 0.368 75 0.367 831 25
0.02 0.37 0.3679
0.01 0.38 0.3679
0.005 0.40 0.3676
0.0025 0.48 0.3680
0.00125 1.28 0.3712
Table 5.4. (e−x ) at x = 1 from Central Finite Difference

Approximation
In the six-digit computations, the optimal value of h is 0.08, yielding a result ac-
curate to three significant figures. Hence three significant figures are lost due to a
combination of truncation and roundoff errors. Above optimal h, the dominant er-
ror is caused by truncation; below it, the roundoff error becomes pronounced. The
best result obtained with the eight-digit computation is accurate to four significant
figures. Because the extra precision decreases the roundoff error, the optimal h is
smaller (about 0.02) than in the six-figure calculations.
5.3 Richardson Extrapolation
Richardson extrapolation is a simple method for boosting the accuracy of certain nu-
merical procedures, including finite difference approximations (we also use it later in
other applications).
Suppose that we have an approximate means of computing some quantity G.
Moreover, assume that the result depends on a parameter h. Denoting the approxi-
mation by g(h), we have G = g(h) + E (h), where E (h) represents the error. Richard-
son extrapolation can remove the error, provided that it has the form E (h) = chp , c
and p being constants. We start by computing g(h) with some value of h, say h = h1 .
In that case we have
p
G = g(h1 ) + ch1 (i)
Then we repeat the calculation with h = h2 , so that

p
G = g(h2 ) + ch2 (j)
189 5.3 Richardson Extrapolation
Eliminating c and solving for G, Eqs. (i) and (j) yield
(h1 / h2 ) p g(h2 ) − g(h1 )

G= (5.8)
(h1 / h2 ) p − 1
which is the Richardson extrapolation formula. It is common practice to use h2 =

h1 /2, in which case Eq. (5.8) becomes
2 p g(h1 /2) − g(h1 )

G= (5.9)
2p − 1
Let us illustrate Richardson extrapolation by applying it to the finite difference

approximation of (e−x ) at x = 1. We work with six-digit precision and use the results
in Table 5.4. Because the extrapolation works only on the truncation error, we must
confine h to values that produce negligible roundoff. In Table 5.4 we have
g(0.64) = 0.380 610 g(0.32) = 0.371 035
The truncation error in central difference approximation is E (h) = O(h2 ) = c1 h2 +

c2 h4 + c3 h6 + . . . . Therefore, we can eliminate the first (dominant) error term if we
substitute p = 2 and h1 = 0.64 in Eq. (5.9). The result is
22 g(0.32) − g(0.64) 4(0.371 035) − 0.380 610

G= = = 0. 367 84 3
22 − 1 3
which is an approximation of (e−x ) with the error O(h4 ). Note that it is as accurate as
the best result obtained with eight-digit computations in Table 5.4.
EXAMPLE 5.1
Given the evenly spaced data points
x 0 0.1 0.2 0.3 0.4

f (x) 0.0000 0.0819 0.1341 0.1646 0.1797
compute f (x) and f (x) at x = 0 and 0.2 using finite difference approximations of
O(h2 ).
Solution. We will use finite difference approximations of O(h2 ). From the forward
difference tables Table 5.3a we get
−3f (0) + 4f (0.1) − f (0.2) −3(0) + 4(0.0819) − 0.1341

f (0) = = = 0.967
2(0.1) 0.2
2f (0) − 5f (0.1) + 4f (0.2) − f (0.3)

f (0) =
(0.1)2
2(0) − 5(0.0819) + 4(0.1341) − 0.1646

= = −3.77
(0.1)2
The central difference approximations in Table 5.1 yield
−f (0.1) + f (0.3) −0.0819 + 0.1646

f (0.2) = = = 0.4135
2(0.1) 0.2
f (0.1) − 2f (0.2) + f (0.3) 0.0819 − 2(0.1341) + 0.1646

f (0.2) = = = −2.17
(0.1)2 (0.1)2
EXAMPLE 5.2
Use the data in Example 5.1 to compute f (0) as accurately as you can.
Solution. One solution is to apply Richardson extrapolation to finite difference ap-

proximations. We start with two forward difference approximations of O(h2 ) for f (0):
one using h = 0.2 and the other one using h = 0.1. Referring to the formulas of O(h2 )
in Table 5.3a, we get
−3f (0) + 4f (0.2) − f (0.4) 3(0) + 4(0.1341) − 0.1797

g(0.2) = = = 0.8918
2(0.2) 0.4
−3f (0) + 4f (0.1) − f (0.2) −3(0) + 4(0.0819) − 0.1341

g(0.1) = = = 0.9675
2(0.1) 0.2
Recall that the error in both approximations is of the form E (h) = c1 h2 + c2 h4 +

c3 h6 + . . .. We can now use Richardson extrapolation, Eq. (5.9), to eliminate the dom-
inant error term. With p = 2 we obtain
22 g(0.1) − g(0.2) 4(0.9675) − 0.8918

f (0) ≈ G = = = 0.9927
22 − 1 3
which is a finite difference approximation of O(h4˙).
EXAMPLE 5.3
b C
B β
c
a
α d
A D
The linkage shown has the dimensions a = 100 mm, b = 120 mm, c = 150 mm, and
d = 180 mm. It can be shown by geometry that the relationship between the angles
α and β is
α (deg) 0 5 10 15 20 25 30
β (rad) 1.6595 1.5434 1.4186 1.2925 1.1712 1.0585 0.9561
191 5.4 Derivatives by Interpolation
If link A B rotates with the constant angular velocity of 25 rad/s, use finite difference
approximations of O(h2 ) to tabulate the angular velocity dβ/dt of link BC against α.
Solution. The angular speed of BC is

dβ dβ dα dβ
= = 25 rad/s
dt dα dt dα
where dβ/dα can be computed from finite difference approximations using the data
in the table. Forward and backward differences of O(h2 ) are used at the endpoints,
and central differences elsewhere. Note that the increment of α is
1 π 2
h = 5 deg rad / deg = 0.087 266 rad
180
The computations yield
−3β(0◦ ) + 4β(5◦ ) − β(10◦ ) −3(1.6595) + 4(1.5434) − 1.4186
β̇(0◦ ) = 25 = 25
2h 2 (0.087 266)
= −32.01 rad/s.
β(10◦ ) − β(0◦ ) 1.4186 − 1.6595
β̇(5◦ ) = 25 = 25 = −34.51 rad/s
2h 2(0.087 266)
etc.
The complete set of results is
α (deg) 0 5 10 15 20 25 30
β̇ (rad/s) −32.01 −34.51 −35.94 −35.44 −33.52 −30.81 −27.86
5.4 Derivatives by Interpolation
If f (x) is given as a set of discrete data points, interpolation can be a very effective
means of computing its derivatives. The idea is to approximate the derivative of f (x)
by the derivative of the interpolant. This method is particularly useful if the data
points are located at uneven intervals of x, when the finite difference approximations
listed in Section 5.2 are not applicable1 .
Polynomial Interpolant
The idea here is simple: Fit the polynomial of degree n
Pn−1 (x) = a 0 + a 1 x + a 2 x 2 + · · · + a n x n
through n + 1 data points and then evaluate its derivatives at the given x. As pointed
out in Section 3.2, it is generally advisable to limit the degree of the polynomial to
less than six in order to avoid spurious oscillations of the interpolant. Because these
oscillations are magnified with each differentiation, their effect can devastating.
1 It is possible to derive finite difference approximations for unevenly spaced data, but they would
not be as accurate as the formulas derived in Section 5.2.
In view of this limitation, the interpolation is usually a local one, involving no more
than a few nearest-neighbor data points.
For evenly spaced data points, polynomial interpolation and finite difference
approximations produce identical results. In fact, the finite difference formulas are
equivalent to polynomial interpolation.
Several methods of polynomial interpolation were introduced in Section 3.2. Un-
fortunately, none are suited for the computation of derivatives of the interpolant.
The method that we need is one that determines the coefficients a 0 , a 1 , . . . , a n of the
polynomial. There is only one such method discussed in Chapter 3: the least-squares
fit. Although this method is designed mainly for smoothing of data, it will carry out
interpolation if we use m = n in Eq. (3.22)—recall that m is the degree of the inter-
polating polynomial and n + 1 represents the number of data points to be fitted. If
the data contain noise, then the least-squares fit should be used in the smoothing
mode; that is, with m < n. After the coefficients of the polynomial have been found,
the polynomial and its first two derivatives can be evaluated efficiently by the func-
tion evalPoly listed in Section 4.7.
Cubic Spline Interpolant

Because of its stiffness, a cubic spline is a good global interpolant; moreover, it is easy
to differentiate. The first step is to determine the second derivatives ki of the spline at
the knots by solving Eqs. (3.11). This can be done with the function curvatures in
the module cubicSpline listed in Section 3.3. The first and second derivatives are
then be computed from
' (
ki 3(x − xi+1 )2
fi,i+1 (x) = − (xi − xi+1 )
6 xi − xi+1
' (
ki+1 3(x − xi )2 yi − yi+1
− − (xi − xi+1 ) + (5.10)
6 xi − xi+1 xi − xi+1
x − xi+1 x − xi
fi,i+1 (x) = ki − ki+1 (5.11)
xi − xi+1 xi − xi+1
which are obtained by differentiation of Eq. (3.10).
EXAMPLE 5.4
Given the data
x 1.5 1.9 2.1 2.4 2.6 3.1

f (x) 1.0628 1.3961 1.5432 1.7349 1.8423 2.0397
compute f (2) and f (2) using (1) polynomial interpolation over three nearest-
neighbor points, and (2) a natural cubic spline interpolant spanning all the data
points.
Solution of Part (1). The interpolant is P2 (x) = a 0 + a 1 x + a 2 x 2 passing through the

points at x = 1.9, 2.1, and 2.4. The normal equations, Eqs. (3.22), of the least-squares
fit are
⎡ ⎤⎡ ⎤ ⎡ ⎤
n xi xi2 a0 yi
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ xi xi3 ⎥ ⎢ ⎥ ⎢ yi xi ⎥
⎣ xi2 ⎦ ⎣ a1 ⎦ = ⎣ ⎦
2 4
xi xi3 xi a2 yi xi2
After substituting the data, we get
⎡ ⎤⎡ ⎤ ⎡ ⎤
3 6.4 13.78 a0 4.6742
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 6.4 29.944 ⎥ ⎢ ⎥ ⎢ ⎥
⎣ 13.78 ⎦ ⎣ a 1 ⎦ = ⎣ 10.0571 ⎦
13.78 29.944 65.6578 a2 21.8385
T
which yields a = −0.7714 1.5075 −0.1930 .
The derivatives of the interpolant are P2 (x) = a 1 + 2a 2 x and P2 (x) = 2a 2 . There-
fore,
f (2) ≈ P2 (2) = 1.5075 + 2(−0.1930)(2) = 0.7355
f (2) ≈ P2 (2) = 2(−0.1930) = −0.3860
Solution of Part (2). We must first determine the second derivatives ki of the spline
at its knots, after which the derivatives of f (x) can be computed from Eqs. (5.10) and
(5.11). The first part can be carried out by the following small program:
#!/usr/bin/python
## example5_4
from cubicSpline import curvatures
from LUdecomp3 import *
import numpy as np
xData = np.array([1.5,1.9,2.1,2.4,2.6,3.1])
yData = np.array([1.0628,1.3961,1.5432,1.7349,1.8423, 2.0397])
print(curvatures(xData,yData))
input("Press return to exit")
The output of the program, consisting of k 0 to k 5 , is
[ 0. -0.4258431 -0.37744139 -0.38796663 -0.55400477 0. ]
Because x = 2 lies between knots 1 and 2, we must use Eqs. (5.10) and (5.11) with
i = 1. This yields
' (
k 1 3(x − x2 )2
f (2) ≈ f1,2

(2) = − (x1 − x2 )
6 x1 − x2
' (
k 2 3(x − x1 )2 y1 − y2
− − (x1 − x2 ) +
6 x1 − x2 x1 − x2
' (
(−0.4258) 3(2 − 2.1)2
= − (−0.2)
6 (−0.2)
' (
(−0.3774) 3(2 − 1.9)2 1.3961 − 1.5432
− − (−0.2) +
6 (−0.2) (−0.2)
= 0.7351
x − x2 x − x1
f (2) ≈ f1,2

(2) = k 1 − k2
x1 − x2 x1 − x2
2 − 2.1 2 − 1.9
= (−0.4258) − (−0.3774) = −0. 4016
(−0.2) (−0.2)
Note that the solutions for f (2) in parts (1) and (2) differ only in the fourth signif-
icant figure, but the values of f (2) are much farther apart. This is not unexpected,
considering the general rule: The higher the order of the derivative, the lower the pre-
cision with which it can be computed. It is impossible to tell which of the two results
is better without knowing the expression for f (x). In this particular problem, the data
points fall on the curve f (x) = x 2 e−x/2 , so that the “true” values of the derivatives are
f (2) = 0.7358 and f (2) = −0.3679.
EXAMPLE 5.5
Determine f (0) and f (1) from the following noisy data:
x 0 0.2 0.4 0.6

f (x) 1.9934 2.1465 2.2129 2.1790
x 0.8 1.0 1.2 1.4
f (x) 2.0683 1.9448 1.7655 1.5891
Solution. We used the program listed in Example 3.10 to find the best polynomial fit
(in the least-squares sense) to the data. The program was run three times with the
following results:
Degree of polynomial ==> 2

Coefficients are:
[ 2.0261875 0.64703869 -0.70239583]
Std. deviation = 0.0360968935809

Coefficients are:
[ 1.99215 1.09276786 -1.55333333 0.40520833]
Std. deviation = 0.0082604082973

Coefficients are:
[ 1.99185568 1.10282373 -1.59056108 0.44812973 -0.01532907]
Std. deviation = 0.00951925073521
Based on standard deviation, the cubic seems to be the best candidate for the
interpolant. Before accepting the result, we compare the plots of the data points and
the interpolant—see the following figure. The fit does appear to be satisfactory.
2.3
2.2
2.1
2.0
f(x)
1.9
1.8
1.7
1.6
1.5
0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40
x
Approximating f (x) by the interpolant, we have
f (x) ≈ a 0 + a 1 x + a 2 x 2 + a 3 x 3
so that
f (x) ≈ a 1 + 2a 2 x + 3a 3 x 2
Therefore,
f (0) ≈ a 1 = 1.093
f (1) = a 1 + 2a 2 + 3a 3 = 1.093 + 2(−1.553) + 3(0.405) = −0. 798
In general, derivatives obtained from noisy data are at best rough approxima-
tions. In this problem, the data represent f (x) = (x + 2)/ cosh x with added random

noise. Thus f (x) = 1 − (x + 2) tanh x / cosh x, so that the “correct” derivatives are
f (0) = 1.000 and f (1) = −0.833.
PROBLEM SET 5.1

1. Given the values of f (x) at the points x, x − h1 , and x + h2 , where h1 = h2 , de-
termine the finite difference approximation for f (x). What is the order of the
truncation error?
2. Given the first backward finite difference approximations for f (x) and f (x), de-
rive the first backward finite difference approximation for f (x) using the oper-

ation f (x) = f (x) .
3. Derive the central difference approximation for f (x) accurate to O(h4 ) by apply-
ing Richardson extrapolation to the central difference approximation of O(h2 ).
4. Derive the second forward finite difference approximation for f (x) from the
Taylor series.
5. Derive the first central difference approximation for f (4) (x) from the Taylor series.
6. Use finite difference approximations of O(h2 ) to compute f (2.36) and f (2.36)
from the following data:
x 2.36 2.37 2.38 2.39

f (x) 0.85866 0.86289 0.86710 0.87129
7. Estimate f (1) and f (1) from the following data:
x 0.97 1.00 1.05

f (x) 0.85040 0.84147 0.82612
8. Given the data

x 0.84 0.92 1.00 1.08 1.16
f (x) 0.431711 0.398519 0.367879 0.339596 0.313486
calculate f (1) as accurately as you can.

9. Use the data in the table to compute f (0.2) as accurately: as possible:
x 0 0.1 0.2 0.3 0.4

f (x) 0.000 000 0.078 348 0.138 910 0.192 916 0.244 981
10. Using five significant figures in the computations, determine d(sin x)/dx at x =
0.8 from (a) the first forward difference approximation and (b) the first central
difference approximation. In each case, use h that gives the most accurate result
(this requires experimentation).
11. Use polynomial interpolation to compute f and f at x = 0, using the data
x −2.2 −0.3 0.8 1.9

f (x) 15.180 10.962 1.920 −2.040
Given that f (x) = x 3 − 0. 3x 2 − 8. 56x + 8. 448, gauge the accuracy of the result.
12.
B
2.5
R R
θ x
A C
The crank A B of length R = 90 mm is rotating at a constant angular speed of

dθ /dt = 5 000 rev/min. The position of the piston C can be shown to vary with
the angle θ as
% &
x = R cos θ + 2.52 − sin2 θ
Write a program that plots the acceleration of the piston at θ = 0◦ , 5◦ , 10◦ , . . . ,

180◦ . Use numerical differentiation to compute the acceleration.
13.
v γ
C
A αB β
a x
The radar stations A and B, separated by the distance a = 500 m, track the plane
C by recording the angles α and β at one-second intervals. If three successive
readings are
t (s) 9 10 11
α 54.80◦ 54.06◦ 53.34◦
β 65.59◦ 64.59◦ 63.62◦
calculate the speed v of the plane and the climb angle γ at t = 10 s. The coordi-
nates of the plane can be shown to be
tan β tan α tan β

x =a y =a
tan β − tan α tan β − tan α
14.
20
D
70 β
Dimensions C
in mm 190
0
19
A
α θ
B 60
Geometric analysis of the linkage shown resulted in the following table relating
the angles θ and β:
θ (deg) 0 30 60 90 120 150

β (deg) 59.96 56.42 44.10 25.72 −0.27 −34.29
Assuming that member A B of the linkage rotates with the constant angular ve-
locity dθ/dt = 1 rad/s, compute dβ/dt in rad/s at the tabulated values of θ . Use
cubic spline interpolation.
15. The relationship between stress σ and strain ε of some biological materials in
uniaxial tension is
dσ
= a + bσ
dε
where a and b are constants (dσ /dε is called tangent modulus). The following
table gives the results of a tension test on such a material:
Strain ε Stress σ (MPa)

0 0
0.05 0.252
0.10 0.531
0.15 0.840
0.20 1.184
0.25 1.558
0.30 1.975
0.35 2.444
0.40 2.943
0.45 3.500
0.50 4.115
Determine the parameters a and b by linear regression.

Numerical Differentiation

Uploaded by

Copyright:

Available Formats

Numerical Differentiation

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Differentiation

Uploaded by

Copyright:

Available Formats

5 Numerical Differentiation

Given the function f (x), compute d n f/dx n at given x

5.2 Finite Difference Approximations

(2h)2 (2h)3 (2h)4 (4)

(2h)2 (2h)3 (2h)4 (4)

We also record the sums and differences of the series:

First Central Difference Approximations

Central difference approximations for other derivatives can be obtained from

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

Table 5.1. Coefficients of Central Finite Difference

First Noncentral Finite Difference Approximations

f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h)

Table 5.2a. Coefficients of Forward Finite Difference

f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

Table 5.2b. Coefficients of Backward Finite Difference

Second Noncentral Finite Difference Approximations

f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h) f (x + 5h)

Table 5.3a. Coefficients of Forward Finite Difference Approximations of O(h2 )

f (x − 5h) f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

Table 5.3b. Coefficients of Backward Finite Difference Approximations of O(h2 )

Errors in Finite Difference Approximations

• Use double-precision arithmetic.

h Six-Digit Precision Eight-Digit Precision

Table 5.4. (e−x ) at x = 1 from Central Finite Difference

5.3 Richardson Extrapolation

Then we repeat the calculation with h = h2 , so that

Eliminating c and solving for G, Eqs. (i) and (j) yield

(h1 / h2 ) p g(h2 ) − g(h1 )

which is the Richardson extrapolation formula. It is common practice to use h2 =

2 p g(h1 /2) − g(h1 )

Let us illustrate Richardson extrapolation by applying it to the finite difference

g(0.64) = 0.380 610 g(0.32) = 0.371 035

The truncation error in central difference approximation is E (h) = O(h2 ) = c1 h2 +

22 g(0.32) − g(0.64) 4(0.371 035) − 0.380 610

x 0 0.1 0.2 0.3 0.4

−3f (0) + 4f (0.1) − f (0.2) −3(0) + 4(0.0819) − 0.1341

2f (0) − 5f (0.1) + 4f (0.2) − f (0.3)

2(0) − 5(0.0819) + 4(0.1341) − 0.1646

The central difference approximations in Table 5.1 yield

−f (0.1) + f (0.3) −0.0819 + 0.1646

f (0.1) − 2f (0.2) + f (0.3) 0.0819 − 2(0.1341) + 0.1646

Solution. One solution is to apply Richardson extrapolation to finite difference ap-

−3f (0) + 4f (0.2) − f (0.4) 3(0) + 4(0.1341) − 0.1797

−3f (0) + 4f (0.1) − f (0.2) −3(0) + 4(0.0819) − 0.1341

Recall that the error in both approximations is of the form E (h) = c1 h2 + c2 h4 +

22 g(0.1) − g(0.2) 4(0.9675) − 0.8918

which is a finite difference approximation of O(h4˙).

Solution. The angular speed of BC is

The complete set of results is

5.4 Derivatives by Interpolation

Cubic Spline Interpolant

which are obtained by differentiation of Eq. (3.10).

x 1.5 1.9 2.1 2.4 2.6 3.1

Solution of Part (1). The interpolant is P2 (x) = a 0 + a 1 x + a 2 x 2 passing through the

After substituting the data, we get

f (2) ≈ P2 (2) = 1.5075 + 2(−0.1930)(2) = 0.7355

f (2) ≈ P2 (2) = 2(−0.1930) = −0.3860

The output of the program, consisting of k 0 to k 5 , is

[ 0. -0.4258431 -0.37744139 -0.38796663 -0.55400477 0. ]