Department of Civil Engineering

CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

Numerical Differentiation
7.0 Numerical Differentiation

The analytic form of the derivative for a given function, in principle, is always possible to
determine. However, in some cases where analytic form is very complicated, the derivative may be
found by numerical approximation. For functions which are described only in terms of tabulated data,
numerical differentiation is the only means of computing a derivative

7.1 Direct Approximation Approach (Curve Fitting)

One of the oldest and most applied problems in mathematics is the problem of constructing an
approximation to a given function from a discrete set of data pairs. The fundamental two-dimensional
curve-fitting problem from a discrete set of (x, y) data pairs (x1, y1) through (xn, yn) is to predict the value
of y when the value of x is specified. There are a number of ways to do this curve-fitting, depending on
the kind of data available and the specific conditions placed on the resulting curve.

In this chapter, we consider methods which do not account for random errors in the x and y
values as in the least square’s approximation. The data will be treated as accurate and determine a curve
that passes through the data points exactly. The most general and powerful tools for curve-fitting is
approximation theory and the simpler but less robust technique is interpolation. Classical methods for
fitting curves involve polynomial models.

7.1.1 POLYNOMIAL INTERPOLATION

If a set of n data values (x1, y1) through (xn, yn) are given to represent y as a single-valued
function of x, one can find a unique polynomial of degree (n-1) that passes through the data points. For
example, we can find a unique straight line through two points, and we can find a unique quadratic that
passes through 3 points.

The polynomial model for y may be given as

y = a0 + a1x + a2x2 +….+an-1x(n-1)

and use the data pairs (x1, y1) to (xn, yn) to write the n equations and form a system of linear algebraic
equations to determine the coefficients ai. These equations are

and the solution of the linear system may be determined by the methods discussed in the previous
modules. However, solving the system of linear equations is an inefficient way of obtaining a
representation for y. Other polynomial models provide a more efficient ways of predicting y for a given
value of x. These models appear quite different from the above model, nevertheless, they produce the
same unique curve through n data points.

[Link] Langrange Interpolation

The Langrange formula is perhaps the best known of the classical formulas for interpolating
among n data pairs (x1, y1) through (xn, yn). It is used to represent y as a function of x in the form

in which the multiplying polynomials Lj are given by

1
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS
The popularity of Langrange formula is well known because it is easy to code and the data are not
required to be specified with x in ascending or descending order. Although the computation of y is simple,
but for large values of n, the method is still not efficient. The Langrange’s form is useful for interpolating
tables and for deriving formulas for numerical differentiation and integration.

Example [Link].1. For the data pairs (x1, y1), (x2, y2) (x3, y3), fit a curve using the Langrange interpolation.
Solution:

When j =1, i = 2, 3; j = 2, i = 1, 3; j = 3, i = 1, 2

Example [Link].2. Compute ln 9.2 when ln 9.0 = 2.1972, ln 9.5 = 2.2513, ln 11.0 = 2.3979

y(x) = 2.1972(x2 − 20.5x + 104.5) − 3.0017(x2 – 20x + 99) + 0.7993(x2 − 18.5x + 85.5)
y(x) = −0.0052x2 + 0.2050x + 0.756

hence,

ln(9.2) = y(9.2) = 2.2192

[Link] Newton’s General Interpolating Formula

Newton’s general interpolating formula for data pairs (x1, y1) through (xn,yn) has
the form

The coefficients aj are generally obtained by computing a set of quantities known as divided
differences. The notation for divided difference is given by

In general,

2
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

Let us denote the divided differences as

The coefficients aj are related to the divided differences by aj = f1,j

If we let fk,j as the elements of an array at the jth sweep of the computation process, where j = 0, 1, 2,
…, n-1, the elements of the computation process may be illustrated in Figure 5.1 using four data pairs.

j=0 j=1 j=2 j=3

y1 = f1,0
f1,1
y2 = f2,0 f1,2
f2,1 f1,3
y3 = f3,0 f2,2

f3,1

y4 = f4,0

Figure [Link].1 Elements fk,j of the jth computation process

Example [Link].1 The elapsed time t(sec) for a car to accelerate to a speed v(mph) from an initial speed
of 30 mph is described by the data in Table [Link].1. Fit a polynomial to describe the function using
Newton’s form. Estimate the time for the car to accelerate to a speed of 50 mph.

Table [Link].1 Car’s speed versus time to accelerate

i 1 2 3 4

v(mph) 30 40 55 70

t(sec) 0.0 1.6 4.8 8.9

Solution: Since t = t(v)

3
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

When v = 50 mph

t = 3.59 sec

[Link] Neville’s Algorithm

Recall that the Newton interpolation involved two stages: the computation of coefficient and the
evaluation of the polynomial. Neville’s algorithm is a method that requires about the same effort as the
first stage of the Newton interpolation. It is used with n pairs of data (x1, y1) through (xn, yn) to obtain y(x)
directly for a given value of x.

The aim of the method is to evaluate a polynomial of degree (n-1) by a process that is similar to
the computation of the coefficients of the Newton polynomial. However, the quantities fk,j in Figure
[Link].1 are now

The interpolated value y(x) for a given x is the last computed quantity y(x) = f1, n-1.

4
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

Example [Link].1. From problem Example [Link].1, Estimate t when v is 50 mph using Neville’s
Algorithm.

Solution:

Therefore

t(50) = 3.63secs

5
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

7.1.2 SPLINE INTERPOLATION

The curves produced by polynomial interpolation pass through the n specified data points and
the curve may exhibit strong oscillations between the data points when polynomials of high degree are
used. As shown in Figure 7.1.2, suppose that the given data points are approximations to a straight line.
By forcing a fourth-degree polynomial to pass through five points, the curve that is produced will deviate
significantly from the line.

Figure [Link]. Polynomial Interpolation

This suggests that the interpolation might be unsatisfactory if we insist on increasing the order of the
polynomial. However, if we keep the order of the polynomial fixed and use it over different intervals, with
the length of the interval decreasing, then interpolation can be very accurate. This type of curve fitting is
known as spline interpolation. The spline interpolation is applied to n ordered pairs of data. We seek (n-
1) curves that connects points 1 and 2, points 2 and 3, … , and points (n-1) and n as shown in Figure
[Link].

b) first order splines a) Second order splines

c) Third order splines

Figure [Link] Spline Interpolation

6
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

[Link] Cubic Splines

The objective of the cubic spline is to derive third order polynomials for each interval between
data points. In addition, the two curves connecting points (k-1) and k and points k and (k+1) are required
to have the same slope at point k so that the resulting curve fit is both continuous and smooth.

7
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

Example [Link] Use natural spline to estimate the elapsed time t(sec) for a car to accelerate from
30 mph to 50 mph
i 1 2 3 4

v 30 40 55 70

t 0 1.6 4.8 8.9

t = 3.612 secs

7.2. Difference Formulas

One approach to numerical differentiation is to fit a curve with a simple form to the function, and
then differentiate the curve-fit function. Another approach is to obtain an approximation of the derivative
directly in terms of the function values. The resulting expression for the approximation is known as the
difference formula. Such formulas are easily useful for solving differential equations and may be derived
from curve fitting or from Taylor Series.

The Taylor series expansion for a function f(x) is

ℎ ′ ℎ2 ℎ3 ℎ𝑛
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + 𝑓 (𝑥) + 𝑓 ′′ (𝑥) + 𝑓 ′′′ (𝑥) + ⋯ + 𝑓 (𝑛) (𝑥)
1! 2! 3! 𝑛!

When the Taylor series is truncated immediately after the term containing the n th derivative, the
error that occurs is said to be of order hn+1 or O(hn+1) and is equal to R n+1 given by
ℎ𝑛+1 𝑓 (𝑛+1) (𝜀)
𝑅𝑛+1 = , 𝑤ℎ𝑒𝑟𝑒 𝑥 < 𝜀 < 𝑥 + ℎ
(𝑛 + 1)!

Simple difference formulas may be obtained by truncating the Taylor series after the first derivative term,
i.e.,
f (x + h) = f (x) ± hf '(x) + O(h2 )

Forward Difference

𝑓(𝑥+ℎ)−𝑓(𝑥)
Equation 7. 2. a 𝑓 ′ (𝑥) = + 𝑂(ℎ)
ℎ

Backward Difference

𝑓(𝑥)−𝑓(𝑥−ℎ)
Equation 7. 2. b 𝑓 ′ (𝑥) = + 𝑂(ℎ)
ℎ

8
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS
Eqs.(7.2.a,b) are first order formulas and have truncation errors of order h. Higher order
approximations may be derived by using additional function values. For example, we write expressions
for two function values by truncating the Taylor series after the second derivative. Thus, we have

ℎ2 ′′
𝑓(𝑥 + ℎ) = 𝑓(𝑥) + ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑂(ℎ3 )
2

4ℎ2 ′′
𝑓(𝑥 + 2ℎ) = 𝑓(𝑥) + 2ℎ𝑓 ′ (𝑥) + 𝑓 (𝑥) + 𝑂(ℎ3 )
2!

By Eliminating the second derivative term, we obtain a second-order forward difference formula given by

−3𝑓(𝑥) + 4𝑓(𝑥 + ℎ) − 𝑓(𝑥 + 2ℎ)

𝑓 ′ (𝑥) = + 𝑂(ℎ2 )
2ℎ
In similar manner, we obtain the second order backward difference formula
3𝑓(𝑥) − 4𝑓(𝑥 − ℎ) + 𝑓(𝑥 − 2ℎ)
𝑓 ′ (𝑥) = + 𝑂(ℎ2 )
2ℎ

Subtraction gives a second order difference formula for the first derivative.

Central Difference

𝑓(𝑥 + ℎ) − 𝑓(𝑥 − ℎ)
𝑓 ′ (𝑥) = + 𝑂(ℎ2 )
2ℎ

Addition of the expansions for f(x + h) and f(x - h) through the third derivative gives a central
difference formula for the 2nd derivative as

𝑓(𝑥 − ℎ) − 2𝑓(𝑥) + 𝑓(𝑥 + ℎ)

𝑓 ′′ (𝑥) = 2
+ 𝑂(ℎ2 )
ℎ

Example 7.2.1 A particle crosses three sensors which are spaced at intervals of 0.3 m. Measurements
of the time to travel from sensor 1 to sensor 2 and from sensor 2 to sensor 3 are 0.0550 and 0.0554
seconds, respectively. Estimate the speed and acceleration of the particle as it crosses the second
sensor.

Solution 1:

Let the subscripts denote sensors, x = position, t = time, v = speed, and a = acceleration.
𝑑𝑥
Speed at sensor 2 = 𝑣2 = at time 𝑡2
𝑑𝑡
𝑑𝑥
Acceleration at sensor 2 = 𝑎2 = at time 𝑡2
𝑑𝑡

9
 Department of Civil Engineering
CENUMES 313 – NUMERICAL SOLUTIONS TO CIVIL ENGINEERING PROBLEMS

10