Transcript CVEN 302, Chapter 11 - Texas A&M University

Chapter 18
Numerical Integration
of Functions
Numerical Integration
 Tabulated data – the accuracy of the integral
is limited by the number of data points
 Continuous function – we can generate as
many f(x) as required to attain the required
accuracy
 Richardson extrapolation and
Romberg integration
 Gauss Quadratures
 Round-off errors may limit the precision of lowerorder Newton-Cotes composite integration formula
 Use Romberg Integration for efficient integration
Romberg Integration
 More efficient methods to achieve
better accuracy have been developed
 Romberg integration - uses Richardson
extrapolation
 Idea behind Richardson extrapolation improve the estimate by eliminating the
leading term of truncation error at
coarser grid levels
Richardson Extrapolation
The exact integral can be represented as
I  I h  E h
This is true for any h = (ba)/n
I  I h1   E h1   I h2   E h2 
Use trapezoidal rule as an example

b  a 2
E
h f ( )
12
E h1  h1
 2
E h2  h2
2

Composite Trapezoidal Rule
4
• Evaluate the integral I   xe 2 x dx
0
h
 f (0 )  f ( 4 )  23847.66
2
h
n  2, h  2  I   f (0 )  2 f ( 2 )  f ( 4 )  12142.23
2
h
n  4 , h  1  I   f (0 )  2 f ( 1)  2 f ( 2 )
2
 2 f ( 3 )  f ( 4 )  7288.79
n  1, h  4  I 
h
 f (0 )  2 f (0.5 )  2 f ( 1)
2
 2 f ( 1.5 )  2 f ( 2 )  2 f ( 2.5 )  2 f ( 3 )
 2 f ( 3.5 )  f ( 4 )  5764.76
  357.12%
  132.75%
  39.71%
n  8, h  0.5  I 
h
 f (0 )  2 f (0.25)  2 f (0.5 )  
2
 2 f ( 3.5 )  2 f ( 3.75)  f ( 4 )  5355.95
  10.50%
n  16, h  0.25  I 
  2.66%
Richardson Extrapolation
 Truncation error for trapezoidal rule
 h1 
E h1   E h2   
 h2 
2
 Substitute into the exact integral
2
 h1 
I  I ( h1 )  E ( h2 )   I ( h2 )  E ( h2 )
 h2 
 Which leads to
I ( h1 )  I ( h2 )
E h2  
2
1  ( h1 / h2 )
Richardson Extrapolation
 Plugging back into I = I(h) + E(h)
I  I ( h2 ) 
1
h1 / h2 
2
1
I ( h2 )  I ( h1 )
 If h2 = h1/2, then
1
I ( h2 )  I ( h1 )
I  I ( h2 )  2
2 1
4
1
I  I ( h2 )  I ( h1 )
3
3
Richardson Extrapolation
 Combine two O(h2) estimates to get an O(h4) estimate
 Can also combine two O(h4) estimates to get an O(h6)
estimate
16
1
16
1
I
I ( h2 ) 
I ( h1 ) 
Im 
Il
15
15
15
15
 Combine two O(h6) estimates to get an O(h8) estimate
64
1
64
1
I
I ( h2 ) 
I ( h1 ) 
Im 
Il
63
63
63
63
 Im and Il are more and less accurate
estimates, respectively
Romberg Integration
 General form is called Romberg Integration
I j ,k 
4 k  1 I j  1, k  1  I j , k  1
4
k 1
1
 j: level of accuracy - j+1 is more accurate (more
segments)
 k: level of integration - k=1 is the original
trapezoidal rule estimate (O(h2)), k=2 is
improved (O(h4)), k=3 corresponds to O(h6), etc.
Romberg Integration
• Accelerated Trapezoidal Rule
I j ,k 
4 k  1 I j  1, k  1  I j , k  1
4
k 1
1
; k  2, 3,
Trapezoid
Simpson' s
Boole' s
k1
O( h2 )
k2
O ( h4 )
k3
O ( h6 )
k4
O ( h8 )
k5
O ( h10 )
h
I 1, 1
I 1, 2
I 1, 3
I 1, 4
I 1, 5
h/ 2
h/ 4
I 2,1
I 3,1
I 2,2
I 3, 2
I 2,3
I 3, 3
I 2,4
h/8
I 4,1
I 4,2
h / 16
I 5,1
4 I j  1, 1  I j , 1
16 I j  1, 2  I j , 2
64 I j  1, 3  I j , 3
256I j  1, 4  I j , 4
3
15
63
255
Romberg Integration
Accelerated Trapezoid Rule
4
I   xe dx  5216.926477
2x
0
Trapezoid Simpson' s
Boole' s
k1
k2
k3
k4
k5
O( h2 )
O ( h4 )
O ( h6 )
O ( h8 )
O ( h10 )
h4
23847.7
8240.41
5499.68
5224.84
5216.95
h2
12142.2
5670.98
5229.14
5217.01
h1
7288.79
5256.75
5217.20
h  0.5
5764.76
5219.68
h  0.25 5355.95

 2.66%  0.0527%  0.0053%  0.00168%  0.00050%
Romberg Integration
Accelerated trapezoidal Rule
» intg = romberg(‘example1’,0,pi,0.00001,2)
I =
0.0000
0.0000
-5.5122

0.0000
-5.1677
0
-3.8758
0
0
0
» intg = Romberg(‘example1’,0,pi,0.00001,3)
I =
0.0000
0.0000
-5.5122
-4.9221
0.0000
-5.1677
-4.9313
0
-3.8758
-4.9461
0
0
-4.6785
0
0
0
» intg = romberg(‘example1’,0,pi,0.00001,4)
I =
0.0000
0.0000
-5.5122
-4.9221
-4.9349
0.0000
-5.1677
-4.9313
-4.9348
0
-3.8758
-4.9461
-4.9348
0
0
-4.6785
-4.9355
0
0
0
-4.8712
0
0
0
0
» intg = romberg(‘example1’,0,pi,0.00001,6)
I =
0.0000
0.0000
-5.5122
-4.9221
-4.9349
0.0000
-5.1677
-4.9313
-4.9348
-4.9348
-3.8758
-4.9461
-4.9348
-4.9348
-4.9348
-4.6785
-4.9355
-4.9348
-4.9348
0
-4.8712
-4.9348
-4.9348
0
0
-4.9189
-4.9348
0
0
0

-4.9308
0
0
0
0
2
x sin(2 x )dx
-4.9348
-4.9348
0
0
0
0
0
-4.9348
0
0
0
0
0
CVEN 302-501
Homework No. 12
• Chapter 17
• Problem 17.3 (25), 17.5 (25)– Hand Calculation
• Chapter 18
• Problem 18.2 (25), 18.4 (25)
• Due on Monday, 11/17/2008 at the beginning of the
period
Gauss Quadrature
b
I   f ( x )dx
a
Assume
I  c0 f a  c1 f b
a and b are limits of integration
Trapezoidal rule should give exact results
for constant and linear functions
Trapezoidal rule
gives exact solution
for constant and
linear functions
Gauss Quadrature
 Now instead of trapezoidal, which has fixed end
points (a,b), let them float
 4 unknowns - x0 ,x1 ,c0 ,c1
1
I   f ( x )dx  c0 f ( x0 )  c1 f ( x1 )
1
 4 equations - constant, linear (had before in
trapezoidal rule), quadratic, cubic
 Integrate from -1 to 1 to simplify math
Trapezoidal vs. Gauss-Quadrature
Exact for constant and
linear functions
Exact for constant, linear,
quadratic and cubic functions
Gauss Quadrature
 The idea is that if we evaluate the
function at certain points (non-uniformly
distributed), and sum with certain
weights, we will get accurate integral
 Evaluation points and weights are
tabulated
Gauss-Legendre Quadrature
 Choose (c0 , c1 , x0 , x1) to yield highest
possible accuracy
Gauss Quadratures
 Newton-Cotes Formulas
 use evenly-spaced functional values
 Gauss Quadratures (Gauss-Legendre
formulas)
 change of variables so that the interval of
integration is [1,1]
 select functional values at non-uniformly
distributed points to achieve higher accuracy
Gauss Quadrature on [a, b]
 To go to [1,1] from other limits [a,b] use linear transformation
 Change from a  x  b to 1  xd  1
x  a0  a1 x d
a  a0  a1 ( 1)
a0  ( a  b ) / 2
 

 a1  ( b  a ) / 2
b  a0  a1 ( 1)
 Coordinate transformation
ab ba
x

xd
2
2
Gauss Quadrature on [a, b]
 Coordinate transformation from [a,b] to [1,1]
ba
ab
x
xd 
2
2
 x d  1  x  a

 xd  1  x  b
a

b
a
f ( x )dx 
t1

1
1
t2
b
1
ba
ab ba
f(
xd 
)(
)dx d   g ( x d )dx d
1
2
2
2
Gauss Quadrature on [1, 1]

1
1
n
f ( x )dx   c i f ( xi )  c0 f ( x0 )  c1 f ( x1 )    c n f ( xn )
i 1
n2:

1
1
f ( x )dx
 c 0 f ( x0 )  c 1 f ( x 1 )
-1
x0
x1
1
• Choose (c0 , c1 , x0 , x1) such that the method
yields “exact integral” for f(x) = x0, x1, x2, x3
Gauss Quadrature on [1, 1]

n2:
1
1
f ( x )dx  c0 f ( x0 )  c1 f ( x1 )
 Exact integral for f = x0, x1, x2, x3
 Four equations for four unknowns
f


 f

f


 f
1
 1   1dx  2  c0  c 1
1
1
 x   xdx  0  c0 x0  c 1 x 1
1
2
 x   x dx   c0 x02  c 1 x12
1
3
1
2
2
1
 x   x 3 dx  0  c0 x03  c 1 x 13
3
1
1
I   f ( x )dx  f ( 
1
1
3
c0  1

c 1  1

1
  x0 
3


1
x

 1
3

) f (
1
3
)
Gauss Quadrature on [1, 1]
n 3:
-1
x0

1
1
f ( x )dx  c0 f ( x0 )  c1 f ( x1 )  c 2 f ( x2 )
x1
x2
1
 Choose (c0, c1, c2, x0, x1, x2) such that the method
yields “exact integral” for f(x) = x0, x1, x2, x3,x4, x5
Gauss Quadrature on [1, 1]
• Exact integral for f = x0, x1, x2, x3, x4, x5
f


f

f


f


f

f

1
 1   1dx  2  c0  c 1  c 2
1
1
 x   xdx  0  c0 x0  c 1 x1  c 2 x 2
1
1
 x   x 2 dx 
2
1
2
 c0 x02  c 1 x12  c 2 x 22
3
1
 x   x 3 dx  0  c0 x03  c 1 x13  c 2 x 23
3
1
1
 x   x 4 dx  0  c0 x04  c 1 x14  c 2 x 24
4
1
c0  5/9

c 1  8/9
c 2  5/9
 
 x0   3/5
 x1  0

 x 2  3/5
1
 x   x 5 dx  0  c0 x05  c 1 x15  c 2 x 25
5
1
I
1
1
5
3
8
5
3
f ( x )dx  f ( 
) f (0 ) f (
)
9
5
9
9
5
Example: Gauss Quadrature
Evaluate
4
I   xe2 x dx  5216.926477
0
Coordinate transformation
ba
ab
x
xd 
 2 xd  2; dx  2dxd
2
2
4
1
I   xe dx   ( 4 xd  4 )e
2x
4 xd  4
1
0
1
dxd   g ( xd )dxd
1
Two-point formula
I

1
1
g ( x d )dxd  g (
1
)  g(
1
)  (4 
4
4
)e
3
3
3
 9.167657324 3468.376279 3477.543936
4
3
 (4 
4
4
)e
3
(   33.34% )
4
3
Example: Gauss Quadrature
 Three-point formula
5
8
5
g (  0.6 )  g (0 )  g ( 0.6 )
1
9
9
9
5
8
5
4  0.6
4
 ( 4  4 0.6 )e
 ( 4 )e  ( 4  4 0.6 )e 4  0.6
9
9
9
5
8
5
 ( 2.221191545)  ( 218.3926001)  (8589.142689)
9
9
9
 4967.106689
(   4.79% )
1
I   g ( xd )dxd 
 Four-point formula
I   g ( xd )dxd  0.34785g ( 0.861136)  g (0.861136)
1
1
 0.652145g ( 0.339981)  g (0.339981)
 5197.54375
(  0.37%)
Points Weighting Factors Function Arguments Truncation Error
c0  1.0000000
x0  0 .577350269
2
 f (4) ( )
c 1  1.0000000
x1  0 .577350269
3
c0  0 .5555556
x0  0 .774596669
 f (6) ( )
c 1  0 .8888889
x1  0 .000000000
c 2  0 .5555556
x 2  0 .774596669
4
c0  0 .3478548
x0  0.861136312
 f (8) ( )
5
6
c1
c2
c3
c0
c1
 0 .6521452
 0 .6521452
 0 .3478548
 0 .2369269
 0 .4786287
x1
x2
x3
x0
x1
 0.339981044
 0.339981044
 0.861136312
 0.906179846
 0.538469310
c2
c3
c4
c0
c1
c2
 0 .5688889
 0 .4786287
 0 .2369269
 0 .1713245
 0 .3607616
 0 .4679139
x2
x3
x4
x0
x1
x2
 0.000000000
 0.538469310
 0.906179846
 0.932469514
 0.661209386
 0.238619186
c 3  0 .4679139
c 4  0 .3607616
c5  0 .1713245
x 2  0.238619186
x 4  0.661209386
x5  0.932469514
 f (10) ( )
 f (12) ( )
Gauss-Legendre
Formulas
Gauss Quadrature
Gauss Quadrature
» I=Gauss_quad(‘example1’,0,pi,2);
t =
-0.5774
-0.7746
-0.8611
-0.9062
0.5774
0
-0.3400
-0.5385
0
0.7746
0.3400
0
0
0
0.8611
0.5385
0
0
0
0.9062
c =
1.0000
0.5556
0.3479
0.2369
1.0000
0.8889
0.6521
0.4786
0
0.5556
0.6521
0.5689
0
0
0.3479
0.4786
0
0
0
0.2369
tt =
-0.5774
0.5774

cd =
2
x
sin( 2 x )dx
1
0
1
» I
I =
-8.6878
» Q=quad8(‘example1’,0,pi)
Q =
-4.9348
Exact Q = -4.9348

k=2
» I=Gauss_quad(‘example1’,0,pi,5);
t =
-0.5774
-0.7746
-0.8611
-0.9062
0.5774
0
-0.3400
-0.5385
0
0.7746
0.3400
0
0
0
0.8611
0.5385
0
0
0
0.9062
c =
1.0000
0.5556
0.3479
0.2369
1.0000
0.8889
0.6521
0.4786
0
0.5556
0.6521
0.5689
0
0
0.3479
0.4786
0
0
0
0.2369
tt =
-0.9062
-0.5385

0
2
0.5385
0
0.9062
cd =
0.2369
0.4786
0.5689
0.4786
0.2369
» I
I =
Exact Q = -4.9348
-4.9333

x sin( 2 x )dx
Gauss Quadrature
k=5
Adaptive Quadrature
 Composite Simpson’s 1/3 rule requires the use
of equally spaced points
 Use adaptive refinement in regions of relatively
abrupt changes
 Estimate truncation error between two levels
of refinement
 Automatically adjust the step size so that small
steps are taken in regions of sharp variations
while larger steps are used elsewhere
 MATLAB functions: quad and quadl
MATLAB Integration Methods
trapz(x,y)
* Composite trapezoid rule
q = quad(‘func’,xmin,xmax)
* Adaptive Simpson’s rule (p 439),
more efficient for low accuracies or nonsmooth functions
q =quadl(‘func’,xmin, xmax)
* Labatto quadrature – more efficient for high
accuracies and smooth functions (p439)