Partt 5

P
Chapter 18

Numerical
N i l IIntegration
t ti off
Functions

PowerPoints organized by Dr. Michael R. Gustafson II, Duke University

Revised by Prof. Jang, CAU
All images copyright © The McGraw-Hill Companies, Inc. Permission required for reproduction or display.
 Chapter Objectives
• U
Understanding
d di how
h Richardson
Ri h d extrapolation
l i
provides a means to create a more accurate
integral estimate by combining two less
accurate estimates.
• Understanding how Gauss quadrature provides
superior integral estimates by picking optimal
abscissas at which to evaluate the function.
function
• Knowing how to use MATLAB’s built-in
functions quad and quadl to integrate
functions.
 Richardson Extrapolation
• Although the composite Simpson 1/3 rule is adequate
for manyy problems,
p , there are more efficient methods
available: Richardson extrapolation, Gaussian
quadrature, and adaptive quadrature.
• Richard
Ri h d extrapolation
t l ti methods
th d use two
t estimates
ti t off an
integral to compute a third, more accurate
approximation.
pp
• If two O(h2) estimates I(h1) and I(h2) are calculated for
an integral using step sizes of h1 and h2, respectively,
an improved
i d O(h4) estimate
ti t may be b formed
f d using:
i
1
I  I(h2 )  I(h2 )  I (h1 )
(h1 / h2 ) 1
2
1
• For the special case where the interval is halved
(h2=h1/2), this becomes: 4 1
I  I(h2 )  I (h1 )
3 3
 Example 18.1
Q. Use Richardson extrapolation to evaluate the integral of the
g function from a = 0 to b = 0.8. Exact solution is
following
1.640533.
For n=1 and 2 I  4 (1.0688)  1 (0.1728)  1.367467
3 3
 Et = 1.640533-1.367467=0.273067 (t = 16.6%)
For n= 2 and 4 I  4 (1.4848)  1 (1.0688)  1.623467
3 3
 Et = 1.640533-1.623467
1 640533 1 623467 = 0.017067
0 017067 (
( t = 1.0%)
1 0%)
Results of single and composite integral using the trapezoidal rule
Segments n Step size h Integral t
1 0.8 0.1728 89.5%
2 04
0.4 1 0688
1.0688 34 9%
34.9%
4 0.2 1.4848 9.5%
Richardson Extrapolation (cont)
• For the cases where there are two O(h4)
estimates and the interval is halved (hm=hl/2),
an improved O(h6) estimate may be formed
using:
16 1
I  I m  Il
15 15

• For the cases where there are two O(h6)

estimates and the interval is halved (hm=hl/2),
an improved
p O(h8) estimate mayy be formed
using: 64 1
I Im  Il
63 63
 Example
p 18.2
Q. In Ex. 18.1 we used Richardson extrapolation to
compute two integral estimate of O(h4). To combine
these two estimates to compute
p the integral
g with
O(h6).
4 1
Sol )
Sol.) I  (1.0688)  (0.1728)  1.367467
3 3
4 1
I  (1.4848)  (1.0688)  1.623467
3 3
16 1
I (1.623467)  (1.367467)  1.640533
15 15
 Which is the exact value of the integral.
 The Romberg Integration
Algorithm
l ih
• Note that the weighting factors for the
Richardson extrapolation add up to 1 and that
as accuracy increases, the approximation using
th smaller
the ll step
t size
i isi given i greatert weight.
i ht
• In general, I j1,k1  I j,k1
k1
4
I 
4 k1 11
j,k

where Ij+1,k-1 and Ij,k-1 are the more and less

accurate integrals,
g respectively.
p y
• Ij,k is the new approximation, where k is the
level of integration (k=1: original trapezoidal
rule,
l k k=2:
2 O(h4) estimates,
ti t k=3:
k 3 O(h6) estimates)
ti t )
and j is used to determine which 4 1
approximation is more accurate. I 1, 2 
3
I 2,1 
3
I1,1
 Romberg Algorithm Iterations
• The chart below shows the process by
which lower level integrations are
combined
bi d tot produce
d more accurate
t
estimates: Romberg Integration
MATLAB Code for Romberg
 Gauss Quadrature
• Gauss quadrature describes a
class of techniques for evaluating
the area under a straight line by Error
joining any two points on a curve
rather than simply choosing the
endpoints.
• The key is to choose the line that
Error is decreased
balances the positive and
negative
ti errors. f (a )  f (b)
I  (b  a )
2

• The
h particular
i l Gauss
G quadrature
d
formulas described in this section
are called Gauss-Legendre
formulas
 Gauss-Legendre
Gauss Legendre Formulas
• The method of undetermined coefficients
offers
ff a third
hi d approach
h (G
(Gauss Quadrature).
Q d )
I  c 0 f (a )  c1 f (b)

• The trapezoidal rule yield exact results when

the function being integrated is a constant or
a straight
t i ht line.
li Two
T simple
i l equations
ti are y=1
1
and y=x.
(b a ) / 2
c 0  c1   1dx
(b  a ) / 2

c 0  c1  b  a
ba ba (ba ) / 2
 c0  c1  xdx
2 2  ( b  a ) / 2 ba ba
I f (a)  f (b)
ba ba 2 2
 c0  c1 0 Equivalent to the trapezoidal rule
2 2
 Two point Gauss-Legendre
g
Formulas
I  c 0 f ( x 0 )  c1 f ( x1 )
- x0 and x1 are not fixed at the end
points ->
points. > four unknowns
unknowns. -> > four
equations are required.
It fits the integral g of these
functions y = constant, y = x, y =
x2, y = x13 1
c 0  c1   1dx
d 2 c 0 x 0  c1 x1   xdx
d 0
1 1

2 1 3
c 0 x 0  c1 x1   x dx  0
1
c 0 x  c x   x dx 
2 2 2 3 3
0 1 1 1 1
3
c0  c1  1
1 1
1 1 I  f ( )  f ( )
x0    0.5773503 x1   0.5773503 3 3
3 3
 Two point Gauss-Legendre
g
Formulas
•Integrals over other intervals require a change in variables
to set the limits from -1 to 1.
a  a 0 (1)  a1 ba ba
a0  a1 
b  a 0 (1)  a1 2 2

(b  a)  (b  a ) x d ba
x dx  dx d
2 2
 Example
p 18.3
Q. Use two point Gauss-Legendre formula to evaluate the integral of
the following function from x = 0 and 0.8.
0 8 The exact value is
1.640533. f ( x)  0.2  25 x  200 x 2  675 x 3  900 x 4  400 x 5
x  0.4  0.4 x d dx  0.4dx d
0.8
0 (0.2  25 x  200 x 2  675 x 3  900 x 4  400 x 5 )dx
1
  [0.2  25(0.4  0.4 x d )  200(0.4  0.4 x d ) 2  675(0.4  0.4 x d ) 3
-1

 900(0.4  0.4 x d ) 4  400(0.4  0.4 x d ) 5 ]0.4dx d

f xd 
xd  1 / 3  f xd  : 0.516741
xd  1 / 3  f xd  : 1.305837
  
I  f xd  1 / 3  f xd  1 / 3  1.822578  I  f(
1
) f(
1
)
1 1
I  f( ) f( ) 3 3
3 3
 Higher-Point
g Formula
I  c0 f x0   c1 f x1     cn1 f xn1 
<Weighting factors and functions used in Gauss-Legendre formulas>
Points Weighting factors Function arguments Truncation errors
c0 = 1.0000000 x0 = –0.577350269
2  f(4)()
c1 = 1.0000000 x1 = 0.577350269
c0 = 0.5555556 x0 = –0.774596669
3 c1 = 0.8888889 x1 = –0.0  f(6)()
c2 = 0.5555556 x2 = 0.774596669
c0 = 0.3478548 x0 = –0.861136312
c1 = 0.6521452 x1 = –0.339981044
4  f(8)()
c2 = 0.6521452 x2 = 0.339981044
c3 = 0.3478548 x3 = 0.861136312
c0 = 0.2369269 x0 = –0.906179846
c1 = 0.4786287 x1 = –0.538469310
5 c2 = 0.5688889 x2 = 0.0  f(10)()
c3 = 0.4786287 x3 = 0.538469310
c4 = 0.2369269 x4 = 0.906179846
c0 = 0.1713245 x0 = –0.932469514
c1 = 0.3607616 x1 = –0.661209386
c2 = 0.4679139 x2 = –0.238619186
6  f(12)()
c3 = 0.4679139 x3 = 0.238619186
c4 = 0.3607616 x4 = 0.661209386
c5 = 0.1713245 x5 = 0.932469514
 Example
p 18.4
Q. Use the three point Gauss-Legendre formula to estimate the
integral for the same function as in Ex
Ex. 18
18.3.
3 The exact value is
1.640533. f ( x )  0.2  25 x  200 x 2  675 x 3  900 x 4  400 x 5

Sol.)
I  0.5555556 f ( 0.7745967 )  0.8888889 f (0)  0.5555556 f (0.7745967 )
I  0.2813013  0.8732444  0.4859876  1.640533

– Because Gauss quadrature requires function evaluations at

nonuniformly spaced points within the integration interval, it is not
appropriate for cases where the function is unknown. -> Not suited
f engineering
for i i problems
bl that
h deal
d l with
i h tabulated
b l d data.
d
– Where the function is known, its efficiency can have advantages.
 Adaptive Quadrature
• M
Methods
th d such h as Simpson’s
Si ’ 1/3 rule
l has
h a
disadvantage in that it uses equally spaced
points - if a function has regions of abrupt
changes, small steps must be used over the
entire domain to achieve a certain accuracy. y
• Adaptive quadrature methods for
integrating
g g functions automaticallyy adjust
j
the step size so that small steps are taken
in regions of sharp variations and larger
steps are taken where the function changes
gradually.
 Adaptive Quadrature in MATLAB
• MATLAB has
h twot built-in
b ilt i ffunctions
ti ffor
implementing adaptive quadrature:
– quad:
q uses adaptive
p Simpson
p quadrature;
q ; possibly
p y more
efficient for low accuracies or nonsmooth functions
– quadl: uses Lobatto quadrature; possibly more efficient
for
o highg accu
accuracies
ac es and
a d smooth
s oot functions
u ct o s
• q = quad(fun, a, b, tol, trace, p1, p2, …)
– fun : function to be integrates
– a, b: integration bo
bounds
nds
– tol: desired absolute tolerance (default: 10-6)
– trace: flag
g to display
p y details or not
– p1, p2, …: extra parameters for fun
– quadl has the same arguments