CHAPTER 4: MATHEMATICAL MODELING WITH MATLAB
Lecture 4.1: Finite
difference approximations for numerical derivatives
Forward, backward, and central differences
for derivatives:
Problem: Given a set of data points near the point (x_{0},y_{0}):
…; (x_{2},y_{2});
(x_{1},y_{1}); (x_{0},y_{0}); (x_{1},y_{1});
(x_{2},y_{2}); …
Suppose
the data values represent values of a function y = f(x). Find a
numerical approximation for derivatives f'(x_{0}), f''(x_{0}),
… of the function y = f(x) at the point (x_{0},y_{0}).
Example: Linear electrical circuits
consist of resistors, capacitors, inductors, and voltage and current sources.
In a simple resistorinductor (RL) oneport network driven by a current source,
the voltage V = V(t) develops across the port terminals when a
current I = I(t) is applied to the input port, The voltage output
V(t) can be determined as a sum of the voltage drop across the
resistor R I(t) and of the voltage drop across the inductor L
I'(t). The derivative I'(t) is to be found from the input
current I(t) measured at different time instances:
Solution:
The first
derivative f'(x_{0}) of the function y = f(x)
at the point (x_{0},y_{0}) can be
approximated by the slope of the secant line that passes through two points
(linear piecewise interpolation). Depending on whether the points are taken to
the right of the point (x_{0},y_{0}) (future data), to the left of the point (x_{0},y_{0})
(past data), or to both sides, the slope of the secant line is called the
forward, backward or central difference approximations.

Forward difference approximation: The secant line passes the points (x_{0},y_{0})
and (x_{1},y_{1}). f'(x_{0}) _{} D_{forward}(f;x_{0})
= _{} Forward
differences are useful in solving initialvalue problems for differential
equations by singlestep predictorcorrector methods (such as Euler methods).
Given the values f'(x_{0}) and f(x_{0}),
the forward difference approximates the value f(x_{1}). 

Backward difference approximation: The secant line passes the points (x_{1},y_{1})
and (x_{0},y_{0}).
f'(x_{0}) _{} D_{backward}(f;x_{0})
= _{} Backward differences are useful for approximating the
derivatives if the data values are available in the past but not in the
future (such as secant methods for root finding and control problems). Given
the values f(x_{1}) and f(x_{0}),
the backward difference approximates the value f(x_{1}),
if it depends on f'(x_{0}). 

Central difference approximation: The secant line passes the points (x_{1},y_{1})
and (x_{1},y_{1}). f'(x_{0}) _{} D_{central}(f;x_{0})
= _{} Central
differences are useful in solving boundaryvalue problems for differential
equations by finite difference methods. Approximating values of f'(x_{0})
that occurs in differential equations or boundary conditions, the central
difference relates unknown values f(x_{1}) and f(x_{1}) by an
linear algebraic equation. 
h = 0.1; x0 = 1; x_1 =
x0h; x1 = x0 + h;
% the three data points are located at equal distance h (the step
size)
y0 = exp(10*x0); y_1 = exp(10*x_1); y1 = exp(10*x1); yDexact =
10*exp(10*x0);
yDforward = (y1y0)/h; % simple form of forward difference for
equally spaced grid
yDbackward = (y0y_1)/h; %
simple form of backward difference
yDcentral = (y1y_1)/(2*h);
% simple form of central difference
fprintf('Exact = %6.2f\nForward = %6.2f\nBackward = %6.2f\nCentral = %6.2f',yDexact,yDforward,yDbackward,yDcentral);
Forward = 378476.76
Backward = 139233.82
Central = 258855.29
Errors of numerical differentiation:
Numerical
differentiation is inherently illconditioned process. Two factors determine
errors induced when the derivative f'(x_{0}) is replaced
by a difference approximations: truncation and rounding errors.
Consider
the equally spaced data points with constant step size: h = x_{1}
– x_{0 }= x_{0} –_{ } x_{1}. The theory based on the Taylor expansion
method shows the following truncation errors:
·
Forward difference approximation:
f'(x_{0}) – D_{forward}(f,x_{0})
= _{}f''(x), x _{} [x_{0},x_{1}]
The
truncation error of the forward difference approximation is proportional to h,
i.e. it has the order of O(h). The error is also proportional to
the second derivative of the function f(x) at an interior point x
of the forward difference interval.
·
Backward difference approximation:
f'(x_{0}) – D_{backward}(f,x_{0})
= _{} f''(x), x _{} [x_{1},x_{0}]
The
truncation error of the backward difference approximation is as bad as that of
the forward difference approximation. It also has the order of O(h)
and is also proportional to the second derivative of the function f(x)
at an interior point x of the backward difference interval.
·
Central difference approximation:
f'(x_{0}) – D_{central}(f,x_{0})
= _{}f'''(x), x _{} [x_{1},x_{1}]
The
truncation error of the central difference approximation is proportional to h^{2}
rather than h, i.e. it has the order of O(h^{2}).
The error is also proportional to the third derivative of the function f(x)
at an interior point x of the central difference interval. The
central difference approximation is an average of the forward and backward
differences. It produces a much more accurate approximation of the derivative
at a given small value of h, compared to the forward and backward
differences. If the data values are available both to the right and to the left
of the point (x_{0},y_{0}), the use of the
central difference approximation is the most preferable.
x0 = 1; x_1 = x0h; x1 = x0+h;
y0 = exp(10*x0); y_1 = exp(10*x_1); y1 = exp(10*x1); yDe =
10*exp(10*x0);
yDf = (y1y0)./h; yDb = (y0y_1)./h; yDc = (y1y_1)./(2*h);
eDf = abs(yDfyDe); eDb = abs(yDbyDe); eDc = abs(yDcyDe);
loglog(h,eDf,'g:',h,eDb,'b',h,eDc,'r');
If
the step size h between two points becomes smaller, the truncation
error of the difference approximation decreases. It decreases faster for
central difference approximations and slower for forward and backward
difference approximations. For example, if h is reduced by a
factor of 10, the truncation error of the central difference
approximation is reduced by a factor of 100, while the truncation
errors of the forward and backward differences are reduced only by a factor of 10.
When
h becomes too small, the difference approximations are taken for
almost equal values of f(x) at the two points. Any rounding error
of computations of f(x) is magnified by a factor of 1/h. As
a result, the rounding error grows with h for very small values
of h. An optimal step size h = h_{opt} can
be computed from minimization of the sum of truncation and rounding errors:
·
Forward difference approximation:
e_{forward} =  f'(x_{0})
– D_{forward}(f,x_{0}) <
M_{2} _{} + 2 _{},
where
M_{2} = max  f''(x). The minimum of error occurs for h
= h_{opt }=2_{}, when e_{forward }= 2_{}.
·
Central difference approximation:
e_{central} =  f'(x_{0})
– D_{central}(f,x_{0}) < M_{3} _{}+ 2 _{},
where
M_{3} = max  f'''(x). The minimum of error occurs for h
= h_{opt }=_{}, when e_{forward }= 3_{}.
M3 = 1000*exp(10*(x0+0.1));
hoptForward = 2*(eps/M2)^(1/2)
hoptCentral = (6*eps/M3)^(1/3)
eoptForward = 2*(eps*M2)^(1/2)
eoptCentral = 3*(eps^2*M3/6)^(1/3)
hoptCentral = 2.8127e008
eoptForward = 7.2924e005
eoptCentral = 2.3683e008
Example: Two central difference
approximations are applied to compute the derivative of the current I'(t):
green pluses are for h = 10 and blue dots are for h = 5.
The exact derivative I'(t) is shown by red solid curve.
Newton forward difference interpolation:
Problem: Given a set of (n+1) data points:
(x_{1},y_{1});
(x_{2},y_{2}); …; (x_{n},y_{n}); (x_{n+1},y_{n+1})
Find
a Newton polynomial of degree n, y = P_{n}(x), that passes through all (n+1) data
points. The Newton polynomial is a discrete Taylor polynomial of the form:
y = P_{n}(x) = c_{0}
+ c_{1} (xx_{1}) + c_{2}(xx_{1})(xx_{2})
+ c_{n }(x – x_{1})(x – x_{2})…(xx_{n})
where the coefficients c_{j} are
constructed from diagonal divided differences.
Comparison:
Lagrange
polynomial interpolation is convenient when the same grid points [x_{1},x_{2},…,x_{n+1}]
are
repeatevely used in several applications. The data values can be stored in
computer memory to reduce a number. The Lagrange interpolation is not useful
however when additional data points are added or removed to improve the
appearance of the interpolating curve. The data set has to be completely
recomputed every time when the data points are added or removed. The Newton
polynomial interpolation is more convenient for adding and deleting additional
data points. It is also more convenient for an algorithmic use of the Horner's
nested multiplication. Of course, the Newton interpolating polynomials
coincides with the Vandermond and Lagrange interpolating polynomials for a
given data set, since the interpolating polynomial y = P_{n}(x) of
order n connecting (n+1) data points are unique.
The difference between Vandermonde, Lagrange, and Newton interpolating
polynomials lies only in computational aspects.
Solution:
The
coefficients c_{j} are to be found from the conditions: P_{n}(x_{k})
= y_{k} that results in a linear system:
Ac
= y. The
Gauss elimination algorithm or induction methods can be used to prove that c_{j}
= d_{j,j}, where d_{j,j} are diagonal
entries in the table of forward divided differences:
grid points 
data values 
first differences 
second differences 
third differences 
fourth differences 
x_{1} 
y_{1} 




x_{2} 
y_{2} 
f[x_{1},x_{2}] 



x_{3} 
y_{3} 
f[x_{2},x_{3}] 
f[x_{1},x_{2},x_{3}] 


x_{4} 
y_{4} 
f[x_{3},x_{4}] 
f[x_{2},x_{3},x_{4}] 
f[x_{1},x_{2},x_{3},x_{4}] 

x_{5} 
y_{5} 
f[x_{4},x_{5}] 
f[x_{3},x_{4},x_{5}] 
f[x_{2},x_{3},x_{4},x_{5}] 
f[x_{1},x_{2},x_{3},x_{4},x_{5}] 
Zero
divided differences f[x_{k}]_{ }are simply the values of the function y =
f(x) at the point (x_{k},y_{k}). First
divided differences f[x_{k},x_{k+1}] are forward
difference approximation for derivatives of the function y = f(x)
at (x_{k},y_{k}):
f[x_{k},x_{k+1}]
= _{}
Second,
third, and higherorder forward divided difference are constructed by using the
recursive rule:
f[x_{k},x_{k+1},…,x_{k+m}]
= _{}
With
the use of divided differences c_{j1} = f[x_{1},x_{2},…,x_{j}],
the Newton interpolating polynomial y = P_{n}(x) has
the explicit representation:
P_{n}(x) = _{}f[x_{1},x_{2},…,x_{j}] _{}
% example of explicit computation of
coefficients of Newton polynomials
x = [ 1,0.75,0.5,0.25,0]; n = length(x)1;
y = [ 14.58,6.15,1.82,0.23,0.00];
A = ones(n+1,1); % the coefficient matrix for
linear system Ac = y
for j = 1 : n
A = [ A,
A(:,j).*(x'x(j))];
end
A, c = A\(y'); c = c'
1.0000 0.2500 0 0 0
1.0000 0.5000
0.1250 0 0
1.0000 0.7500
0.3750 0.0938 0
1.0000 1.0000
0.7500 0.3750 0.0938
c = 14.5800 33.7200 32.8000 14.5067 0.2133
function [yi,c]
= NewtonInter(x,y,xi)
% Newton interpolation algorithm
% x,y  rowvectors of (n+1) data values (x,y)
% xi  a rowvector of xvalues, where interpolation is to be
found
% yi  a rowvector of interpolated yvalues
% c  coefficients of Newton interpolating polynomial
n = length(x)  1; % the degree of interpolation polynomial
ni = length(xi); % the number of xvalues, where interpolation is
to be found
D = zeros(n+1,n+1); % the matrix for Newton divided differences
D(:,1) = y'; % zeroorder divided differences
for k = 1 : n
D(k+1:n+1,k+1) =
(D(k+1:n+1,k)D(k:n,k))./(x(k+1:n+1)x(1:nk+1))';
end
c = diag(D);
% computation of values of the Newton interpolating polynomial at
values of xi
% the algorithm uses the Horner's rule for polynomial evaluation
yi = c(n+1)*ones(1,ni); % initialization of the vector yi as the
coefficient of the highest degree
for k = 1 : n
yi =
c(n+1k)+yi.*(xix(n+1k)); % nested multiplication
end
x = [
1,0.75,0.5,0.25,0]; y = [ 14.58,6.15,1.82,0.23,0.00];
xInt = 1 : 0.01 : 0;
[yInt,c] = NewtonInter(x,y,xInt);
c = c', plot(xInt,yInt,'g',x,y,'b*');
c = 14.5800 33.7200
32.8000 14.5067 0.2133
MATLAB
finite differences:
Let
the data points x_{1},x_{2},…,x_{n},x_{n+1}
be equally spaced with constant step size h = x_{2} – x_{1}.
Then, the divided differences can be rewritten as:
f[x_{k},x_{k+1},…,x_{k+m}]
= _{},
where
_{} is the mth forward
difference of a function y = f(x) at the point (x_{k},y_{k}).
_{} = y_{1} – y_{0}
_{} = y_{2 }– 2 y_{1} + y_{0}
_{} = y_{3 }–3 y_{2} + 3 y_{1} – y_{0}
_{} = y_{4 }– 4 y_{3} + 6 y_{2 } 4 y_{1} – y_{0}
_{} = y_{5} – 5 y_{4}
+ 10 y_{3} – 10 y_{2} + 5 y_{1} – y_{0}
Derivatives
of the interpolation polynomials y = P_{n}(x) approximate
derivatives of the function y = f(x). Matching the nth derivative
of the polynomial y = P_{n}(x) with f^{(n)}(x_{0}),
we find the forward difference approximation for higherorder
derivatives:
f^{(n)}(x_{0})
_{}_{}
·
diff(y): computes firstorder forward difference for a given vector y
·
diff(y,n): computes nth order forward difference for a given vector
y
·
gradient(u): computes horizontal and vertical forward differences for a given matrix
to approximate the x and yderivatives of u(x,y):
grad(u) = [u_{x},u_{y}]
·
del2(u): computes a discrete Laplacian for a given matrix u: _{}u = u_{xx} + u_{yy}
n = 50; x =
linspace(0,2*pi,n); h = x(2)x(1); y = sin(x); % function
y1 = diff(y)/h; y2 = diff(y,2)/h^2; y1ex = cos(x); y2ex =
sin(x);% derivatives
plot(x(1:n1),y1,'b',x,y1ex,':r',x(1:n2),y2,'g',x,y2ex,':r');
Hierarchies of higherorder difference
approximations:
Newton
interpolating polynomials and divided difference tables can be constructed for
backward differences, since the order of data points x_{1},x_{2},…,x_{n},x_{n+1}
is arbitrary. By arranging for data points in descenting order, the Newton
polynomial represents the backward differences. It is trickier to construct
tables for central differences. Since central differences are the most accurate
approximations, special algorithms are designed to automate derivation of
coefficients of central difference approximations for higherorder derivatives.
·
hierarchy of forward differences:
Let
the data points are equally spaced with constant step size h. Fix
a point (x_{0},y_{0}) and present a forward
difference approximation for the nth derivative as the inner
product multiplication:
f^{(n)}(x_{0})
_{}_{}D_{n}*y
where
y = [ y_{0},y_{1},y_{2}, …]
n = 100; A = diag(ones(n1,1),1)diag(ones(n,1));
m = 8; B = A;
for k = 1 : m1
D(k,:) = B(1,1:m);
B = B*A;
end
D = D
1 2
1 0 0 0 0
0
1 3
3 1 0
0 0 0
1 4
6 4 1 0 0
0
1 5
10 10 5
1 0 0
1 6
15 20 15
6 1 0
1 7 21 35 35 21 7 1
·
hierarchy of central differences
Let
the data points are equally spaced with constant step size h. Fix
a point (x_{0},y_{0}) and present a central
difference approximation for the nth derivative as the inner
product multiplication:
f^{(2m1)}(x_{0})
_{}_{}D_{2m1}*y;
f^{(2m)}(x_{0}) _{}_{}D_{2m}*y
where
y = [ …,y_{2},y_{1},y_{0},y_{1},y_{2},
…]
n = 100; A = diag(ones(n1,1),1)diag(ones(n1,1),1);
A2 = diag(ones(n1,1),1)+diag(ones(n1,1),1)2*diag(ones(n,1));
m = 4; B = A; C = A2; k = 1;
while (k < (2*m2) )
D(k,:) = B(m,1:2*m1);
D(k+1,:) = C(m,1:2*m1);
B = B*A2; C = C*A2;
k = k+2;
end
D = D
0 0
1 2 1 0 0
0 1
2 0 2 1 0
0 1
4 6 4
1 0
1 4
5 0 5
4 1
1 6 15 20 15 6 1
Richardson extrapolation for higherorder
differences:
Recursive
difference formulas for derivatives can be obtained by canceling the truncation
error at each order of numerical approximation. This method is called the
Richardson extrapolation. It can be used only if the data values are equally
spaced with constant step size h.
·
Recursive forward differences:
The
hierarchy of forward differences for first and higherorder derivatives has the
truncation error of order O(h). Denote the forward difference
approximation for a derivative f^{(n)}(x_{0}) as^{
}D_{1}(h). Compute the approximation with two step sizes h
and 2h:
f^{(n)}(x_{0})
= D_{1}(h) + _{}h; f^{(n)}(x_{0})
= D_{1}(2h) + 2_{}h
where
_{} is unknown
coefficient for the truncation error. By cancelling the truncation error of
order O(h), we define a new forward difference approximation for
the same derivative:
f^{(n)}(x_{0})
= 2D_{1}(h) – D_{1}(2h) = D_{2}(h)
The
new forward difference approximation D_{2}(h) for the
same derivative is more accurate since the truncation error is of order O(h^{2}).
The
firstorder approximation D_{1}(h) is a twopoint divided
difference, while the secondorder approximation D_{2}(h)
is a threepoint divided difference:
D_{1}(h) = _{}, D_{2}(h) = _{}
The
firstorder approximation D_{1}(h) is a threepoint
divided difference, while the secondorder approximation D_{2}(h)
is a fourpoint divided difference:
D_{1}(h) = _{}, D_{2}(h) = _{}
The
process can be continued to find a higherorder forwarddifference
approximation D_{m}(h) with the truncation error O(h^{m}).
The recursive formula for Richardson forward difference extrapolation:
D_{m+1}(h) =D_{m}(h)
+ _{}
·
Recursive central differences:
The
hierarchy of central differences for first and higherorder derivatives has the
truncation error of order O(h^{2}). Denote the central
difference approximation for a derivative f^{(n)}(x_{0}) as^{
}D_{1}(h). Compute the approximation with two step sizes h
and 2h:
f^{(n)}(x_{0})
= D_{1}(h) + _{}h^{2}; f^{(n)}(x_{0})
= D_{1}(2h) + 4_{}h^{2}
where
_{} is unknown
coefficient for the truncation error. By cancelling the truncation error of
order O(h^{2}), we define a new central difference
approximation for the same derivative:
f^{(n)}(x_{0})
= _{}
The
new central difference approximation D_{2}(h) for the
same derivative is more accurate since the truncation error is of order O(h^{4}).
The
firstorder approximation D_{1}(h) is a threepoint
divided difference, while the secondorder approximation D_{2}(h)
is a fivepoint divided difference:
D_{1}(h) = _{}, D_{2}(h) =
_{}
The
firstorder approximation D_{1}(h) is a threepoint
divided difference, while the secondorder approximation D_{2}(h)
is a fivepoint divided difference:
D_{1}(h) = _{}, D_{2}(h) =
_{}
The
process can be continued to find a higherorder centraldifference
approximation D_{m}(h) with the truncation error O(h^{2m}).
The recursive formula for Richardson forward difference extrapolation:
D_{m+1}(h) =D_{m}(h)
+ _{}
·
Numerical algorithm
In
order to compute the central difference approximation of a derivative f^{(n)}(x_{0})
up to order m, central difference approximations of lower
order are to be computed with larger step sizes: h, 2h, 4h, 8h, …,
(m1)h.
These
approximations can be arranged in a table of recursive derivatives:
step size 
D_{1} 
D_{2} 
D_{3} 
D_{4} 
D_{5} 
h 
D_{1}(h) 




2h 
D_{1}(2h) 
D_{2}(h) 



4h 
D_{1}(4h) 
D_{2}(2h) 
D_{3}(h) 


8h_{} 
D_{1}(8h) 
D_{2}(4h) 
D_{3}(2h) 
D_{4}(h) 

16h 
D_{1}(16h) 
D_{2}(8h) 
D_{3}(4h) 
D_{4}(2h) 
D_{5}(h) 
The
diagonal entries are values of higherorder central difference approximations
for the derivative f^{(n)}(x_{0})
where
x_{0} is the central point surrounded by the points x
= x_{0 } 2^{k1}h and x = x_{0 }+ 2^{k1}h
for k = 1,2,…,m.
The
higherorder approximation D_{k}(h) has the truncation
error O(h^{2k}). If h is small, the
truncation error rapidly decreases with larger k. However, the
rounding error grows with larger value of k. An optimal order m
exists at the minimum of the sum of the truncation and rounding errors.
Example:
The figure
below presents the central difference approximations D_{1}(h) and
D_{2}(h) for the derivative of the current I(t) with
the step size h = 10. Blue circles are found by fivepoint
central differences D_{2}(h), green pluses are obtained
by threepoint central differences D_{1}(h), and the
exact derivative I'(t) is shown by a red solid curve. Fivepoint
difference approximation D_{2}(h) is clearly more
accurate than the threepoint difference D_{1}(h).