CS321

Numerical Analysis
Lecture 4

Numerical Integration

Professor Jun Zhang

Department of Computer Science
University of Kentucky
Lexington, KY 40506‐0633

October 6, 2015
 Definite Integral

A definite integral has an interval for integration.  For a fixed integration 
interval, the result is a number


0
2
sin x dx  1

An indefinite integral does not have an integration interval.  The result  
of an indefinite integral (antiderivative) is a class of functions

 sin x dx   cos x  C
Numerical integration is for computing definite integrals

Fundamental Theorem of Calculus:
b
a F ' ( x ) dx  F (b)  F (a )
x
a F ' ( t ) dt  F ( x )  F (a )

2
Numerical Integration
 Partition

The definite integral of a function can be viewed as the area under a 
curve.  This point of view lends us means to compute definite integral

Let P be a partition of the interval of [a,b] as

P  a  x 0  x1  x 2    x n 1  x n  b

We have n subintervals as [xi,xi+1].  Let mi be the greatest lower bound of 

(a nonnegative function) f(x) on [xi,xi+1] as
m i  inf  f ( x ) : x i  x  x i 1 
and Mi as the least upper bound on the same subinterval

M i  sup f ( x ) : x i  x  x i 1 

4
Partition
 Lower and Upper Sums

The lower sums and upper sums of f corresponding to the given partition 
P is n 1
L( f ; P )   m i  x i 1  x i 
i 0
n 1
U( f ; P)   M i  x i 1  x i 
i 0

If we consider the definite integral of a nonnegative f as the area under 
the curve, we have
b
L( f ; P )  a f ( x ) dx  U ( f ; P )

for all partitions P

If f is continuous on [a,b], then the above inequality defines the definite 
integral.  The integral also exists if f is monotone (either increasing or 
decreasing) on [a,b]
6
Upper and Lower Bounds
 Riemann‐Integrable Functions

If the greatest lower bound equals the least upper bound for all partitions of 
[a,b], i.e.,
inf U ( f ; P )  sup L( f , P )
P P

Then f is said to be Riemann‐integrable

Every continuous function defined on a closed and bounded interval of the 
real line is Riemann‐integrable

We have
b
lim L( f ; Pn )  a f ( x ) dx  lim U ( f ; Pn )
n  n 

where P0, P1,.., Pn,… are a sequence of partitions such that the length of the 

largest subinterval in Pn converges to 0 as n → ∞

We can construct nested (refined) partitions
8
 Computation

1.) need a procedure to evaluate f(x)

2.) determine a partition (how many subintervals) of the interval [a,b]

3.) compute mi and Mi on each subinterval

4.) compute the sums L(f; P) and U(f; P)

5.) an approximate value is obtained

a f ( x ) dx  12 U  f ; P   L f ; P 
b

6.) the error of this approximation is bounded by
1
2
U  f ; P   L f ; P 

9
 Trapezoid Rule

A strategy that is better than estimating both the upper and the lower 
bounds of the area beneath the curve is to use trapezoids

The interval [a,b] is first divided into subintervals [xi,xi+1], 0 ≤ i ≤ n – 1.  A 

typical trapezoid has the subinterval [xi,xi+1] as its base, and the two 
vertical sides are f(xi) and f(xi+1).  The area is given by the base times the 
average height.  The basic trapezoid rule for the subinterval [xi,xi+1] is

 f ( x i )  f ( xi 1 )( xi1  xi )
xi 1
xi
f ( x) dx  1
2

The total area under the curve is
n 1
b
a f ( x ) dx  T ( f ; P )   Ai
i 0

1 n 1
  ( x i 1  x i )[ f ( x i )  f ( x i 1 )]
2 i 0
10
Trapezoid Rule
 Uniform Spacing

The lengths of the subintervals in a partition can be different.  For fast 
computation, a uniform partition of the interval may be advantageous

Let n be the number of subintervals, then h = (b – a)/n is the uniform 

interval spacing.  The nodal points are xi = a + i h, i = 0, 1,…, n.  Hence 
the composite trapezoid rule is
h n 1
T ( f ; P )   [ f ( x i )  f ( x i 1 )]
2 i 0

Note that the end point of an interval is the starting point of the next 
interval.  This fact can save almost half of the computation, 

b  n 1 1 
a f ( x ) dx  h  f ( x i )  [ f ( x 0 )  f ( x n )]
 i 1 2 

12
Trapezoid Rule with Uniform Spacing
 Error Analysis (1)
If f”(x) exists and is continuous on [a,b], the error of the composite 
trapezoid rule T is
b a 2
a f ( x ) dx  T   b12 h f " ( )  O( h 2 )
for some ξ in (a,b)

Proof. We first prove the result for a=0, b=1 and h=1. That is

f ( x) dx  12  f (0)  f (1)   121 f " ( )

1
0

This simplified formula will be proved with the help of polynomial 
interpolation

Define a polynomial of degree one that interpolates f at 0 and 1

p( x)  f (0)  [ f (1)  f (0)]x

14
 Error Analysis (2)
It follows that
1 1
0
p ( x ) dx  f (0)  [ f (1)  f (0)]
2
1
 [ f (0)  f (1)]
2

Using the error formula for the polynomial interpolation, we have

f ( x )  p( x )  1
2 f " [ ( x )] x ( x  1)
Integrate it on both sides,
1 1 1 1
0
f ( x)dx  p( x)dx   f "[ ( x)]x( x  1)dx
0 2 0

15
 Error Analysis (3)
Using the Mean‐Value Theorem for Integrals
1 1
 f "
0
[ ( x )] x ( x  1) dx  f " [ ( s )]
0 x( x  1)dx
1
 f " ( )
6
So we have
1 1 1
0
f ( x)dx  [ f (0)  f (1)]  
2 12
f " ( )
We then do a change of variable, and let

g (t )  f [a  t (b  a)], x  a  (b  a)t
dx  (b  a)dt, g ' (t )  f '[a  t (b  a)](b  a)
g" (t )  f "[a  t (b  a)](b  a) 2
16
 Error Analysis (4)
Then using the result for the special case,
b 1
a
f ( x )dx  (b  a )  f [a  t (b  a )]dt
0
1
 (b  a )  g (t )dt
0

1 1 
 (b  a )  [ g (0)  g (1)]  g " ( ) 
2 12 
ba (b  a )3
 [ f ( a )  f (b)]  f " ( )
2 12
This is the error formula for the trapezoid rule with only one subinterval.

Let [a,b] be divided into n equal subintervals by points

x0 , x1 ,..., xn , with subinterval [ xi , xi 1 ]

17
 Error Analysis (5)
Let h be the interval length
xi 1 h 1
xi
f ( x ) dx  [ f ( xi )  f ( xi 1 )]  h 3 f " ( )
2 12
Sum over all subintervals to get the composite trapezoid rule

n1
h n1 h3 n1
f (x)dx   f (x)x   [ f (xi )  f (xi1)]  f "(i )
b xi1
a
i 0

xi 2 i0 12 i0
Note that h=(b-a)/n, we use Intermediate-Value Theorem of Continuous
Functions,

h 3 n 1 b  a 2  1 n 1  ba 2
  f " ( i )   h   f " ( i )   h f " ( )
12 i 0 12  n i 0  12

18
 Example 1
Show that
ah
a f ( x ) dx 

 f (a )  f (a  h)  h12 f " (a )  

3
h
2

Need to define
t
F ( t )  a f ( x ) dx
We can expand F(a+h) using Taylor series
h2 h3
F(a  h)  F(a)  hF' (a)  F"(a)  F"' (a) 
2 3!
Then by the Fundamental Theorem of Calculus, we have F’(t) = f(t).
Note that F(a) = 0, F”(t)=f’(t), F”’(t)=f”(t), and so on.

19
 Example 2
We have
ah h2 h3
(1) a  f ( x) dx hf (a) 
2
f ' (a) 
3!
f " (a)  
We can also apply the Taylor series directly on f(t) as
h2 h3
(2) f (a  h)  f (a )  hf ' (a )  f " (a)  f " ' (a)  
2 3!
Adding f(a) on both sides of (2) and multiplying it by h/2, we obtain
h h2 h3
(3) [ f (a)  f (a  h)]  hf (a)  f ' (a)  f "(a) 
2 2 4
Subtracting (3) from (1), we finally get

ah h 1 3

a
f (x) dx [ f (a)  f (a  h)]   h f "(a) 
2 12
20
 Estimate Grid Spacing

Example. If the composite trapezoid rule is used to compute
1  x2

0
e dx
with an error of at most 0.5 × 10-4, what is the uniform grid spacing h?

From the graph of the second derivative (a decreasing function)
2
f " ( x )  (4 x 2  2)e  x
We find that
f " ( x )  f " ( 0)  2
We need
ba 2
 h f " ( )  16 h 2  0.5  10  4
12

It follows that h ≤ 0.01732.  The number of subintervals is n ≥ [1/h] = 58

21
 Recursive Trapezoid Idea

What can we  do if the initial partition of interval is not fine enough?

22
 Recursive Trapezoid Formula

Given a parameter n, dividing [a,b] into 2n equally spaced subintervals, we 

have
n 1
h
T ( f ; P )  h  f ( x i )   f ( x 0 )  f ( x n )
i 1 2
n 1
h
 h  f ( a  i h)   f ( a )  f ( b )
i 1 2
Note that n = 2n and h = (b – a)/2n
2 n 1
h
R ( n ,0 )  h  f ( a  i h )   f ( a )  f ( b )
i 1 2
Notice that R(n,0) can be viewed as dividing each subinterval of 
R(n – 1,0) into two equal sub‐subintervals.  If we already computed 
R(n – 1,0), how can we compute R(n,0) cheaply?

23
 Recursive Formula

If R(n – 1,0) is available, R(n,0) can be computed as

2 n1 1
1
R(n,0)  R(n  1,0)  h  f [a  (2k  1)h]
2 k 1

For n ≥ 1 using h = (b – a)/2n.  Initial starting value is

R(0,0)  12 (b  a )[ f (a )  f (b)]
The trick is to only sum the function values at every other grid points

Proof. Note that
2 n 1
R ( n ,0 )  h  f ( a  i h )  C
i 1

with C = h[f(a) + f(b)]/2 and

2 n  1 1
R( n  1,0)  2h  f (a  2 j h)  2C
j 1
24
 Recursive Formula

Hence, we have
2 n 1 2 n1 1
1
R (n,0)  R (n  1,0)  h  f (a  ih)  h  f (a  2 jh)
2 i 1 j 1

2 n1
 h f [a  (2k  1)h]
k 1

Each term in the first sum that corresponds to an even value of index is 
cancelled by a term in the second term. The final result has only the odd 
values of index

We can use the recursive trapezoid formula to compute a sequence of 
approximations to a definite integral using the trapezoid rule, without 
recomputing the values at points that have already been computed in the 
previous step

25
 Two Dimensional Integration
For one dimensional numerical integration on [0,1], using uniform space 
h = 1/n h n 1
i

1
0 f ( x ) dx 
2
[ f ( 0 )  2
i 1
f    f (1)]
n
n
i
  Ai f  
i0 n
For two dimensional integration on a unit square
n
1 1 1 i 
0 0 f ( x , y ) dx dy  0  Ai f  , y  dy
i 1 n 
n
1  i 
  i 0 n  dy
A f  , y
i 0  
n n
 i j
  i  j  n,n
A A f
i 0 j 0  
n n
 i j
   i j  n,n
A A f
i 0 j 0   26
 Romberg Algorithm

Recursive composite trapezoid method
R(0,0)  12 (b  a )[ f (a )  f (b )]
For h = 1/2n and n ≥ 1 2 n 1
1
R(n,0)  R(n  1,0)  h f [a  (2k  1)h]
2 k 1
Using Richardson extrapolation, we can have
R( i , j )  R( i , j  1)
 j
1
[ R( i , j  1)  R( i  1, j  1)]
4 1

For i ≥ j and j ≥ 1.  This is the Romberg algorithm, which may yield better 

approximate values for larger j

27
 Deriving Romberg Algorithm

Composite trapezoid rule on 2n-1 subintervals
b
a f ( x ) dx  R( n  1,0)  a 2 h 2  a 4 h 4  a 6 h 6  
with h = (b – a)/2n-1 and the coefficients ai depend on f but not on h

After one refinement and replacing n – 1 with n and h with h/2, we have

b
a f ( x ) dx 
R( n,0)  14 a 2 h 2  16
1
a 4 h 4  64
1
a6 h6  

Subtracting the 1st equation from 4 times the 2nd equation

b
a
f ( x ) dx  R (n ,1)  14 a 4 h 4  165 a 6 h 6  
where for n ≥ 1
R( n,1)  R( n,0)  13 [ R( n,0)  R( n  1,0)]

28
 More Romberg Algorithm
We could apply the extrapolation idea repeatedly to get

b
a
f ( x) dx R(n,2)  413 a6 h 6  4215 a8 h8  
where

R (n,2)  R(n,1)  151 [ R(n,1)  R (n  1,1)]

This time, the truncation error is of sixth order

A few steps of extrapolation may generate very accurate 
approximations

Too many extrapolations may make the computation tedious 

29
 General Extrapolation

Extrapolation processes can be applied in more general cases where the error 
term can be represented as
E  a h  b h   c h  
with 0 < α < β < γ, we show how the first term of the error expansion is 
annihilated.  Let
L   ( h)  a h  b h   c h   (1)

Replacing h by h/2 yields

L   h2   a h2   b h2   c h2    ( 2)

  

Multiplying (2) by 2α
2 L  2  h2   a h  2 b h2   2 c h2    ( 3)
 

30
 General Extrapolation – Cont.

Subtracting the previous two equations, we can remove the hα term
 h
( 2  1) L  2      ( h)
 2
 ( 2    1)b h   ( 2   1)c h  
We can write the new approximation formula as
2  h 1
L       ( h)  bh   ch  
2 1  2 2 1
This approximation formula raises the order of truncation error from 
O(hα) to O(hβ) with α < β

Please read the book on p. 217 for a concrete example to show how the 
approximation accuracy is improved using extrapolation

31
 Basic Simpson’s Rule

Simple trapezoid rule uses two points for approximations
Can we get more accurate approximation using more points? 
32
 Basic Simpson’s Rule

A three point numerical integration rule using the middle point of the 
interval is known as the Simpson’s rule with different weights for each point

a
a  2h
f ( x ) dx  h
3
 f ( a )  4 f ( a  h)  f ( a  2 h )

Using Taylor’s expansion, we can find the error term of this approximation 
is
h5 ( 4)
 f ( )
90
For some point ξ in (a,a+2h).  This should be compared to the error term of 
the simple trapezoid rule O(h3)

It is desirable to subdivide the interval adaptively so that refinement is only 
placed in the area of large fluctuation of function value

33
 Basic Simpson’s Rule

Comparing composite trapezoid rule and Simpson’s rule

34
 Basic Simpson’s Rule

Using Taylor series for f(x) at a, we have 
1 2 1 3 1 4 (4)
f (ah)  f (a) hf'(a)  h f "(a)  h f "'(a)  h f (a)
2! 3! 4!
If we replace the interval size h by 2h, we obtain

4
4 2
f (a2h)  f (a)2hf'(a)2h2 f "(a) h3F"'(a) h4 f (4) (a)
3 4!
By combining these two expansions, we get

20 4 (4)
f (a)4f (ah) f (a2h) 6f (a)6hf'(a)4h f "(a)2h f "'(a) h f (a)
2 3

4!

35
 Basic Simpson’s Rule

On the other hand, define
x
F ( x)  a
f ( t ) dt
We expand F(a+2h) as
5
4 2 2
F(a2h)  F(a)2hF'(a)2h2F"(a) h3F"'(a) h4 f (4) (a) h5F(5) (a)
3 3 5!
Note that F’=f, F(a)=0, F”=f’, F”’=f”, we have
5
a2h 4 2 2

a
f (x)dx2hf(a)2h2 f '(a) h3 f "(a) h4 f "'(a) h5 f (4) (a)
3 3 54!
From the previous page, we have
h 4
[ f ( a )  4 f ( a  h)  f ( a  2h)]  2hf ( a )  2h 2 f ' ( a )  h 3 f " ( a )
3 3
2 4 20 5 ( 4 )
 h f " ' (a)  h f (a)  
3 3  4!
36
 Basic Simpson’s Rule (1)

Subtracting the previous two equations, we have

a2h h h5 (4)
a
f (x)dx [ f (a)  4 f (a  h)  f (a  2h)] f 
3 90
We have the Simpson’s rule as
a2h b a  a b
a
f (x)dx
6
[ f (a) 4 f 
 2 
  f (b)]
The error term of the Simpson’s rule is

1  b a  (4)
5

   f ()
90 2 

37
 Adaptive Simpson’s Algorithm

Reduce the size of the intervals to get more accurate approximations
38
 Adaptive Simpson’s Algorithm

Given an interval [a,b], we can use the basic Simpson’s rule to compute 
an approximation to the integral as
b
I  a f ( x ) dx  S (a , b )  E (a , b )
where the approximation part is

(b  a )  a  b 
S (a , b)  f ( a )  4 f    f ( b )
6   2  

and the error term is
1 b  a  5
(4)
E (a , b )   90 2
f ( )
For simplicity, we assume f(4)(x) remains constant on (a,b).  Let h = b – a, 
we have
I  S (1 )  E (1 )
For the first step approximation with
S (1)  S ( a , b )
39
 Adaptive Simpson – (II)
And
5
1  h
E (1 )     f (4)
90  2 

We then subdivide the interval [a,b] and apply the basic Simpson’s rule 
on the subintervals [a,c] and [c,b] respectively.  We have a new 
approximation on [a,b] as the sum of two separate approximations
I  S ( 2)  E ( 2)
where c = (a + b)/2 with
S ( 2 )  S (a , c )  S (c , b)
and 5 5
1  h / 2  ( 4 ) 1  h / 2  ( 4 ) 1 (1 )
E ( 2)     f    f  E
90  2  90  2  16
This is certainly a better approximation since the subintervals are smaller 
than the original interval
20
 Adaptive Simpson – (III)

Subtracting the two approximations yields
S ( 2 )  S (1)  E (1)  E ( 2 )  15E ( 2 )
Hence the numerical integration can be
I  S ( 2 )  E ( 2 )  S ( 2 )  15
1
( S ( 2 )  S (1 ) )
The error term is then computable and can be used for building the 
adaptive process
1 ( 2)
S  S (1)  
15

If this test shows that the error is larger than ε, the interval [a,b] can be 
split into two subintervals [a,c] and [c,b] with c = (a + b)/2.  The 
previously described procedure is replaced by ε/2 to make sure that the 
error sum is smaller than ε

41
 Adaptive Simpson’s Algorithm

Refine the intervals at the places where the function changes quickly
42
 Adaptive Simpson – (IV)
Numerical integration on subintervals
b c b
I   f ( x) dx   f ( x) dx   f ( x) dx  I L  I R
a a c

Let S  be the sum of  S L( 2 )  on [a,c ] and S R( 2 )  on [c,b], we have

I  S  I L  I R  S L( 2)  S R( 2 )
 I L  S L( 2 )  I R  S R( 2)
 151 S L( 2)  S L(1)  151 S R( 2 )  S R(1)
If we want to have
IS 

It is more than enough to have
1
S L( 2 )  S L(1)  
15 2

and
1
S R( 2 )  S R(1)  
15 2

43
 Adaptive Simpson’s Algorithm

One decision to make is to choose where to refine the interval
44
 Computational Procedure

The interval [a,b] is divided into four subintervals of equal length.  Two 
Simpson approximations are computed using two double‐width subintervals 
and four single‐width subintervals
h a  b 
S1  f ( a )  4 f    f ( b )
6   2  

h ac c  b 
S2  f ( a )  4 f    2 f ( c )  4 f    f ( b )
12   2   2  


If |S2-S1| ≤15ε, we have done, and set
1
S  16 S 2  S1 
15
Otherwise the interval [a,b] is divided in half and the recursive procedure is 
applied on the two subintervals [a,c] and [c,b], until either the error tolerance 
is satisfied or the maximum number of subdivisions is reached
45
 Adaptive Simpson’s Algorithm

Which sub‐interval or both to divide for refinement
46
 Gaussian Quadrature Formulas
A general numerical integration formula is
b
a f ( x ) dx  A0 f ( x0 )  A1 f ( x1 )    An f ( x n )
It suffices to know the nodes x0, x1,…, xn and the weights A0, A1,…, An.  
For important special functions, they are listed in some reference books

Suppose a set of nodes is given, how to find the weights.  This can be 
done using Lagrange interpolation polynomial as
n
p( x )   f ( x i )l i ( x )
i 0
With
n  x  xj 
li ( x)    
 
j o , j i  xi  x j 

b
If p is a good approximate to 
b
f, we anticipate   p( x ) dx is a good 
a
approximate to a f ( x ) dx

47
 Gaussian Quadrature

We integrate over p(x) as
b b
a f ( x ) dx  a p( x ) dx
n n
  f ( xi )  li ( x) dx   Ai f ( xi )
b

a
i 0 i 0
where we can compute
b
Ai  a l i ( x ) dx
Note that the polynomial interpolation is exact for a polynomial of 
degree at most n.  It follows that the integration will be exact for such 
polynomials

If the nodes can be chosen carefully, it is possible to increase the order 
of polynomial with the exact integration remarkably.  This was 
discussed by Karl Gauss

48
 An Example

Determine a quadrature formula when the interval is [‐2,2]  and the 
nodes are ‐1, 0, and 1.
We first need to compute the cardinal functions
 x  x j   x  0  x  1  1
l0 ( x )    
 (1)  0   1  1   2 x( x  1)

j 1, 2  0x  x    
j 

 x  x j   x  1  x  1 
l1 ( x)    
    1 x
2

j  0 , 2 x1  x j 
  0  (1)  0  1 
 x  x j   x  (1)  x  0  1
l2 ( x )    
 1  (1)  1  0   2 x( x  1)
 x
j  0 ,1 2  x    
j 

The weights are computed by integrating these functions
2
2 1 2 1 1 1  8
A0   l0 ( x) dx   x( x  1)dx   x 3  x 2  
2 2 2 2 3 2  2 3

49
 An Example

Similarly, we have
2
2 2  1  4
A1   l1 ( x) dx   (1  x )dx   x  x 3   
2
2 2
 3  2 3
2
2 1 2 1 1 1  8
A2   l2 ( x) dx   x( x  1)dx   x 3  x 2  
2 2 2 2 3 2  2 3
So the quadrature formula defined on the interval [‐2,2] and using the 
node ‐1, 0, 1, is
2 8 4 8
2 f (x)dx 3 f (1)  3 f (0)  3 f (1)
It can be verified that this formula gives exact values for the three 
functions 
f (x) 1, x, x2

50
 Gaussian Quadrature Theorem

Let q be a nontrivial polynomial of degree n + 1 such that

b k
a x q( x ) dx  0 (0  k  n)
Let x0, x1,…, xn be zeroes of q.  Then we define the formula
n
b b
a f ( x ) dx   Ai f ( x i ), Ai  a l i ( x ) dx
i 0
With these xi’s as nodes, the approximation will be exact for all 
polynomials of degree at most 2n+1.  All these nodes lie in the open 
interval (a,b)

We can first figure out the quadrature forn
1
1 f ( t ) dt   Ai f ( t i )
i 0
Then use the transformation t = [2x – (b – a)]/(b – a) for a Gaussian 
quadrature on the general interval [a,b]

51
 Gaussian Theorem Proof
The transformed integral is
b
a f ( x ) dx 
1
2
( b  a )
1

1
f 12 (b  a )t  12 (b  a ) dt
Proof of Gaussian Quadrature Theorem:
Let f be any polynomial of degree at most (2n + 1).  Dividing f by q with 
quotient p and a remainder r
f  pqr

Both p and r are polynomials of degree at most n

By hypothesis, we have
b
a q( x ) p( x ) dx  0
Since xi are roots of q, we have
f ( x i )  p( x i )q( x i )  r ( x i )  r ( x i )
52
 Gaussian Theorem Proof (II)
b
Since the degree of r is at most n, the integration a r ( x ) dx is exact
b b b
a f ( x ) dx  a p( x )q( x ) dx  a r ( x ) dx
n n
b
 a r ( x ) dx   Ai r ( x i )   Ai f ( x i )
i 0 i 0
Gaussian Quadrature Theorem guarantees high accuracy numerical 
integration with just a few nodes.  However, finding these nodes is not an 
easy task.  The roots of Legendre polynomials are the nodes for Gaussian 
quadrature on the interval [-1,1].  With q0(x) = 1, q1(x) = x, we have for
n≥2

 2n  1   n  1
qn ( x )    xq n 1 ( x )   q n  2 ( x )
 n   n 

53