Numerical Dierentiation and Integration
Problem Statement: If the values of a function f are given at a few points, say, x0 , x1 ,
, xn , we attempt to estimate a derivative f (c) or an integral ab f (x)dx.
Basics of Numerical Dierentiation
Richardson Extrapolation
Basics of Numerical Integration
Quadrature Formulas
Trapezoidal Rule
1
Simpsons 3
Rule
3
Simpsons 8
Rule
Booles Rule
Basics of Numerical Dierentiation
f (x + h) f (x)
f (x) = limith0 (1)
h
1
f (x) [f (x + h) f (x)] (2)
h
h2
f (x + h) = f (x) + hf (x) + f () (3)
2
1 h
f (x) = [f (x + h) f (x)] f () (4)
h 2
1
f (x) [f (x + h) f (x h)] (5)
2h
h2 h3
f (x + h) = f (x) + hf (x) + f (x) + f (3) () (6)
2 3!
h2
h3 (3)
f (x h) = f (x) hf (x) + f (x) f ( ) (7)
2 3!
1 h2
f (x) = [f (x + h) f (x h)] [f (3) () + f (3) ( ))] (8)
2h 12
Examples:
1. f (x) = cos(x), evaluate f (x) at x = /4 with h = 0.01, h = 0.005
2. g(x) = ln(1 + x), evaluate g (x) at x = 1 with h = 0.01, h = 0.005.
3. t(x) = tan1 x, evaluate t (x) at x = 2 with h = 0.01, h = 0.005.
True Solutions
1. f (/4) = sin(/4) = 12 = 0.707106781
1
2. g (1) = 1+1 = 0.500000000
1
3. t ( 2) = 1+( 2)2
= 13 = 0.333333333
Approximation by Taylor Expansion
Theorem 1: Assume that f C 3 [a, b] and x h, x, x + h [a, b] with h > 0. Then
f (x + h) f (x h)
f (x) = D0 (h) (9)
2h
Furthermore, c [a, b] such that
f (x + h) f (x h)
f (x) + E1 (f, h) (10)
2h
2 f (3) (c)
where E1 (f, h) = h 6
= O(h2 ) is called the truncation error.
Theorem 2: Assume that f C 5 [a, b] and x 2h, x h, x [a, b] with h > 0. Then
f (x + 2h) + 8f (x + h) 8f (x h) + f (x 2h)
f (x) = D1 (h) (11)
12h
Furthermore, c [a, b] such that
f2 + 8f1 8f1 + f2
f (x) + E2 (f, h) (12)
12h
4 f (5) (c)
where E2 (f, h) = h 30
= O(h4 )
h By Theorem 1 By Theorem 2 Richardson
0.1 -0.716161095 -0.717353703
0.01 -0.717344150 -0.717356108
0.001 -0.717356000 -0.717356167
0.0001 -0.717360000 -0.717360833
f (0.8) -0.717356091 -0.717356091 -0.717356091
Table 1: Approximating the derivative of f (x) = cos(x) at x = 0.8
Richardsons Extrapolation
Recall that
f (x0 ) D0 (h) + Ch2 , f (x0 ) D0 (2h) + 4Ch2 (13)
Then
4D0 (h) D0 (2h) f2 + 8f1 8f1 + f2
f (x0 ) = D1 (h) (14)
3 12h
Similarly,
f2 + 8f1 8f1 + f2 h4 f (5) ()
f (x0 ) = + D1 (h) + Ch4 (15)
12h 30
f4 + 8f2 8f2 + f4 h4 f (5) ( )
f (x0 ) = + D1 (2h) + 16Ch4 (16)
12h 30
Then
16D1 (h) D1 (2h)
f (x0 ) (17)
15
Richardsons Extrapolation Theorem: Let Dk1 (h) and Dk1 (2h) be two approxi-
mations of order O(h2k ) for f (x0 ), such that
f (x0 ) = Dk1 (h) + c1 h2k + c2 h2k+2 + (18)
f (x0 ) = Dk1 (2h) + 22k c1 h2k + 22k+2 c2 h2k+2 + (19)
Then an improved approximation has the form
4k Dk1 (h) Dk1 (2h)
f (x0 ) = Dk (h) + O(h2k+2 ) = + O(h2k+2) (20)
4 1
k
Dierentiation Approximation Algorithms
Algorithm: Dierentiation Using Limits
Generate the numerical sequence
f (x + 2j h) f (x 2j h)
f (x) Dj = f or j = 0, 1, , n
2 (2j h)
until |Dn+1 Dn | |Dn Dn1 | or |Dn Dn1 | < tol, a user-specied tolerance,
which attempts to nd the best approximation f (x) Dn .
Algorithm: Dierentiation Using Extrapolation
To approximate f (x) by generating a table of D(j, k) for k j using f (x) D(n, n)
as the nal answer. The D(j, k) entries are stored in a lower- matrix. The rst column
is
f (x + 2j h) f (x 2j h)
D(j, 0) =
2 (2j h)
and the elements in row j are
D(j, k 1) D(j 1, k 1)
D(j, k) = D(j, k 1) + f or 1 k j.
4k 1
Matlab Codes for Richardson Extrapolation
%
% Script file: richardson.m
% Example: format long
% richardson(@cos,0.8,0.00000001,0.00000001)
% richardson(@sinh,1.0,0.00001,0.00001)
%
% Richardson Extrapolation for numerical differentiation
% P.333 of John H. Mathews, Kurtis D. Fink
% Input f is the function input as a string f
% delta is the tolerance error
% toler is the tolerance for the relative error
% Output D is the matrix of approximate derivatives
% err is the error bound
% relerr is the relative error bound
% n is the coordinate of the best approximation
%
function [D, err, relerr, n]=richardson(f,x,delta,toler)
err=1.0;
relerr=1.0;
h=1.0;
j=1;
D(1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
while (relerr>toler & err>delta & j<12)
h=h/2;
D(j+1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
for k=1:j
D(j+1,k+1)=D(j+1,k)+(D(j+1,k)-D(j,k))/(4^k-1);
end
err=abs(D(j+1,j+1)-D(j,j));
relerr=2*err/(abs(D(j+1,j+1))+abs(D(j,j))+eps);
j=j+1;
end
[n n]=size(D);
Basics of Numerical Integration
1 2 /2 1
0
2
ex dx = 2
2 x2
0 1 4
dx = 2
0 sin(x)dx = 2
4 14
1 xdx = 3
x t2 /2
erf (x) = 2 e dt, 0 < x < , erf (1) =?
0
x
(x) = 0 1 + cos2 tdt, 0 < x < , () =?
x sint
(x) = 0 t
dt, 0 < x < , (1) =?
x t2
(x) = 0 1+ 4(4t2 )
dt, 0 < x < 2, (2) =?
x t3
(x) = 0 et 1 dt, 0 < x < , (5) = 4.8998922
Quadrature Formulas
The general approach to numerically compute the denite integral ab f (x)dx is by evalu-
ating f (x) at a nite number of sample and nd an interpolating polynomial to approx-
imate the integrand f (x).
Denition: Suppose that a = x0 < x1 < x2 < < xn = b. A formula of the form
n
Q[f ] = i f (xi ) = 0 f (x0 ) + 1 f (x1 ) + + n f (xn ) (21)
i=0
such that b
f (x)dx = Q[f ] + E[f ]
a
is called a quadrature formula. The term E[f ] is called the truncation error for inte-
gration. The values {xj }nj=0 are called the quadrature nodes, and {j }nj=0 are called the
weights.
Closed Newton-Cotes Quadrature Formula
Assume that xi = x0 + i h are equally spaced nodes and fi = f (xi ). The rst four
closed Newton-Cotes quadrature formula are
x1
h
(trapezoidal rule) f (x)dx (f0 + f1 ) (22)
x0 2
x2
1 h
(Simpson s rule) f (x)dx (f0 + 4f1 + f2 ) (23)
3 x0 3
x3
3 3h
(Simpson s rule) f (x)dx (f0 + 3f1 + 3f2 + f3 ) (24)
8 x0 8
x4
2h
(Boole s rule) f (x)dx (7f0 + 32f1 + 12f2 + 32f3 + 7f4 ) (25)
x0 45
Approximating Integrals by Trapezoidal and
Simpsons Rules
Example: The arc length of a curve f (x) = 23 x3/2 between x [0, 1] can be computed by
1
2
= 1 + xdx = (2 2 1) = 1.21895142
0 3
Trapezoidal Rule
h
[f (x0 ) + 2f (x1 ) + 2f (x2 ) + + 2f (xn1 ) + f (xn )]
2
where 0 = x0 < x1 < x2 < < xn1 < xn = 1, h = xi = xi+1 xi 0 i n 1
1.21894654, when n = 50; 1.21895020, when n = 100; = 1.21895142
Simpsons 1/3 Rule
h
[f0 + 4f1 + 2f2 + 4f3 + 2f4 + + 2f2m2 + 4f2m1 + f2m ]
3
where 0 = x0 < x1 < x2 < cdots < xn1 < x2m = 1, h = xi = xi+1 xi 0 i
2m 1
1.21895133, when n = 12; 1.21895140, when n = 20; = 1.21895142
Matlab Codes for Simpsons 1/3 Rule
%
% Simpson.m - Simpsons 1/3 Rule for \sqrt(1+x)
%
format long
n=20;
h=1/n;
x0=0; x1=x0+h; x2=x1+h;
s=0;
for i=0:2:n-2,
f0=sqrt(1+x0);
f1=sqrt(1+x1);
f2=sqrt(1+x2);
s=s+f0+4*f1+f2;
x0=x2;
x1=x2+h;
x2=x2+h+h;
end
s=h*s/3.0;
Simpson Approximated Arc length is, s
