Transcript CVEN 302, Chapter 11

Chapter 17
Numerical
Integration Formulas
y  f(x)
Integration
I
M
lim  f ( x )x  
max x  0 i 1
i
M
i
b
a
f ( x )dx
A   f ( xi ) xi  I
i 1
Graphical Representation of Integral
Integral = area
under the curve
Use of a grid to
approximate an integral
Use of strips to
approximate an integral
Numerical Integration
Survey of land
area of an
irregular lot
Cross-sectional area
and volume flowrate
in a river
Net force
against a
skyscraper
Pressure Force on a Dam
Water exerting pressure on the upstream face of a dam:
(a) side view showing force increasing linearly with depth;
(b) front view showing width of dam in meters.
p = gh = h
Integration
 Weighted sum of functional values at discrete
points
 Newton-Cotes closed or open formulae
-- evenly spaced points
 Approximate the function by Lagrange
interpolation polynomial
 Integration of a simple interpolation polynomial
 Guassian Quadratures
 Richardson extrapolation and Romberg
integration
Basic Numerical Integration
 Weighted sum of function values

b
a
n
f ( x )dx   c i f ( xi )
i 0
f(x)
x0
 c0 f ( x0 )  c1 f ( x1 )    c n f ( xn )
x1
xn-1
xn
x
Numerical Integration
• Idea is to do integral in small parts, like the way
you first learned integration - a summation
• Numerical methods just try to make it faster
and more accurate
12
10
8
6
4
2
0
3
5
7
9
11
13
15
Numerical integration
Newton-Cotes formulas
- based on idea
b
b
a
a
I   f ( x )dx   fn ( x )dx
 Approximate f(x) by a polynomial
fn ( x)  a0  a1 x    an1 x
n1
 an x
n
 fn (x) can be linear
 fn (x) can be quadratic
 fn (x) can also be cubic or other
higher-order polynomials
 Polynomial can be piecewise over the data
Numerical Integration
 Newton-Cotes Closed Formulae -- Use
both end points





Trapezoidal Rule : Linear
Simpson’s 1/3-Rule : Quadratic
Simpson’s 3/8-Rule : Cubic
Boole’s Rule : Fourth-order*
Higher-order methods*
 Newton-Cotes Open Formulae -- Use
only interior points


midpoint rule
Higher-order methods
Closed and Open Formulae
(a) End points are known
(b) Extrapolation
Trapezoidal Rule
• Straight-line approximation

b
a
1
f ( x )dx   c i f ( x i )  c0 f ( x0 )  c 1 f ( x1 )
i 0
h
  f ( x0 )  f ( x 1 )
2
f(x)
L(x)
x0
x1
x
Trapezoidal Rule
• Lagrange interpolation
L( x ) 
x  x0
x  x1
f ( x0 ) 
f ( x1 )
x0  x 1
x 1  x0
xa
dx
let a  x0 , b  x 1 ,  
, d 
; hba
ba
h
x  a    0

  L( )  ( 1   ) f ( a )  ( ) f ( b )
 x  b    1

b
a
b
1
a
0
f ( x )dx   L( x )dx  h L( )d
1
1
0
0
 f ( a ) h  ( 1   ) d   f ( b ) h   d
 f ( a ) h( 

2
2
1

)  f ( b)h
0
2
2
1
0
h
  f ( a )  f ( b )
2
Example:Trapezoidal Rule
• Evaluate the integral
• Exact solution

4
0

4
0
xe 2 x dx
4
 x 2x 1 2x 
xe dx   e  e 
4
2
0
2x
1
1 2x
 e ( 2 x  1)  5216.926477
4
0
• Trapezoidal Rule
I
4
0
4 0
 f (0 )  f ( 4 )  2(0  4 e 8 )  23847.66
xe dx 
2
5216.926  23847.66

 357.12%
5216.926
2x
Better Numerical Integration
 Composite integration
 Multiple applications of Newton-Cotes
formulae
 Composite Trapezoidal Rule
 Composite Simpson’s Rule
 Richardson Extrapolation
 Romberg integration
Apply trapezoidal rule to multiple
segments over integration limits
Two segments
Three segments
7
7
6
6
5
5
4
4
3
3
2
2
1
1
0
0
3
5
7
9
11
13
15
3
5
Four segments
7
9
11
13
15
Many segments
7
7
6
6
5
5
4
4
3
3
2
2
1
1
0
0
3
5
7
9
11
13
15
3
5
7
9
11
13
15
Multiple Applications of
Trapezoidal Rule
Composite Trapezoidal Rule

b
a
x1
x2
xn
x0
x1
xn  1
f ( x )dx   f ( x )dx   f ( x )dx    
f ( x )dx
h
 f ( x0 )  f ( x 1 )   h  f ( x 1 )  f ( x 2 )     h  f ( x n 1 )  f ( x n ) 
2
2
2
h
  f ( x0 )  2 f ( x1 )    2f ( x i )    2 f ( x n1 )  f ( x n )
2

f(x)
ba
h
n
x0
h
x1
h
x2
h
x3
h
x4
x
Trapezoidal Rule
 Truncation error (single application)
1
3
Et  
f ( )( b  a )
12
 Exact if the function is linear ( f = 0)
 Use multiple applications to reduce the
truncation error
( b  a) 3
Ea  
12n 3
n
 f ( )
i 1
( b  a) 3

f  ;
2
12n
i
1 n
f    f ( i )
n i 1
Approximate
error
Composite Trapezoidal Rule


0
2
x sin( 2 x )dx
function f = example1(x)
% a = 0, b = pi
f=x.^2.*sin(2*x);
Composite Trapezoidal Rule
»
»
»
I
a=0; b=pi; dx=(b-a)/100;
x=a:dx:b; y=example1(x);
I=trap('example1',a,b,1)
=
-3.7970e-015
» I=trap('example1',a,b,2)
I =
-1.4239e-015
» I=trap('example1',a,b,4)
I =
-3.8758
» I=trap('example1',a,b,8)
I =
-4.6785
» I=trap('example1',a,b,16)
I =
-4.8712
» I=trap('example1',a,b,32)
I =
-4.9189
» I=trap('example1',a,b,64)
I =
-4.9308
» I=trap('example1',a,b,128)
I =
-4.9338
» I=trap('example1',a,b,256)
I =
-4.9346
» I=trap('example1',a,b,512)
I =
-4.9347
» I=trap('example1',a,b,1024)
I =
-4.9348
» Q=quad8('example1',a,b)
Q =
-4.9348
MATLAB
function


0
x 2 sin( 2 x )dx
n=2
I = -1.4239 e-15
Exact = -4. 9348


0
x 2 sin( 2 x )dx
n=4
I = -3.8758
Exact = -4. 9348


0
x 2 sin( 2 x )dx
n=8
I = -4.6785
Exact = -4. 9348


0
x 2 sin( 2 x )dx
n = 16
I = -4.8712
Exact = -4. 9348
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%
Composite Trapezoidal Rule
»
»
»
»
»
»
»
»
x=0:0.04:4; y=example2(x);
x1=0:4:4; y1=example2(x1);
x2=0:2:4; y2=example2(x2);
x3=0:1:4; y3=example2(x3);
x4=0:0.5:4; y4=example2(x4);
H=plot(x,y,x1,y1,'g-*',x2,y2,'r-s',x3,y3,'c-o',x4,y4,'m-d');
set(H,'LineWidth',3,'MarkerSize',12);
xlabel('x'); ylabel('y'); title('f(x) = x exp(2x)');
» I=trap('example2',0,4,1)
I =
2.3848e+004
» I=trap('example2',0,4,2)
I =
1.2142e+004
» I=trap('example2',0,4,4)
I =
7.2888e+003
» I=trap('example2',0,4,8)
I =
5.7648e+003
» I=trap('example2',0,4,16)
I =
5.3559e+003
Composite Trapezoidal Rule
4
I   xe 2 x dx
0
Simpson’s 1/3-Rule
• Approximate the function by a parabola

b
a
2
f ( x )dx   c i f ( x i )  c0 f ( x0 )  c 1 f ( x 1 )  c 2 f ( x 2 )
i 0
h
  f ( x0 )  4 f ( x 1 )  f ( x 2 )
3
L(x)
f(x)
x0
h
x1
h
x2
x
Simpson’s 1/3-Rule
( x  x0 )( x  x 2 )
( x  x 1 )( x  x 2 )
L( x ) 
f ( x0 ) 
f ( x1 )
( x0  x 1 )( x0  x 2 )
( x 1  x0 )( x 1  x 2 )
( x  x0 )( x  x 1 )

f ( x2 )
( x 2  x0 )( x 2  x 1 )
ab
2
x  x1
ba
dx
h
, 
, d 
2
h
h
 x  x0    1

 x  x1    0
x  x    1
2

let
x0  a, x 2  b, x 1 
L( ) 
 (  1)
2
f ( x0 )  ( 1   2 ) f ( x 1 ) 
 (  1)
2
f ( x2 )
Simpson’s 1/3-Rule
L( ) 

 (  1)
2
f ( x0 )  ( 1   2 ) f ( x 1 ) 
 (  1)
2
f ( x2 )
h 1
f ( x )dx  h L( )dξ  f ( x0 )  ξ ( ξ  1)dξ
1
2 1
1
h 1
2
 f ( x1 )h ( 1  ξ )dξ  f ( x 2 )  ξ ( ξ  1)dξ
0
2 1
1
b
a
1
1
ξ
ξ
h ξ
 f ( x0 ) (  )  f ( x1 )h( ξ  )
3 1
2 1
2 3
3
2
3
1
ξ
h ξ
 f ( x2 ) (  )
2 1
2 3
3

b
a
2
h
f ( x )dx   f ( x0 )  4 f ( x1 )  f ( x 2 )
3
Composite Simpson’s Rule
Piecewise Quadratic approximations
ba
h
n
f(x)
…...
x0 h x1 h x2 h x3 h
x4
xn-2 xn-1
xn
x
Composite Simpson’s 1/3 Rule
 Applicable only if the number of segments is even
Composite Simpson’s 1/3 Rule
 Applicable only if the number of segments is even
x2
x4
xn
x0
x2
xn 2
I   f ( x )dx   f ( x )dx    
f ( x )dx
 Substitute Simpson’s 1/3 rule for each integral
f ( x0 )  4 f ( x 1 )  f ( x 2 )
f ( x2 )  4 f ( x3 )  f ( x4 )
 2h
6
6
f ( x n 2 )  4 f ( x n 1 )  f ( x n )
   2h
6
I  2h
 For uniform spacing (equal segments)
n 1
n 2

( b  a) 
I
 f ( x0 )  4  f ( x i )  2  f ( x j )  f ( x n ) 
3n 
i  1, 3 , 5
j 2 , 4 ,6

Simpson’s 1/3 Rule
 Truncation error (single application)
1 5 (4)
( b  a) 5 ( 4 )
ba
Et  
h f ( )  
f ( ) ; h 
90
2880
2
 Exact up to cubic polynomial ( f (4)= 0)
 Approximate error for (n/2) multiple
applications
(b  a )5 (4)
Ea  
f
4
180n
Composite Simpson’s 1/3 Rule
Evaluate the integral I   xe 2 x dx
0
• n = 2, h = 2
4
h
I   f (0 )  4 f ( 2 )  f ( 4 )
3
2
 0  4 ( 2 e 4 )  4 e 8  8240.411    57.96%
3


• n = 4, h = 1
h
I   f (0 )  4 f ( 1)  2 f ( 2 )  4 f ( 3)  f ( 4 )
3
1
 0  4 ( e 2 )  2( 2 e 4 )  4 ( 3e 6 )  4 e 8
3
 5670.975    8.70%


Simpson’s 3/8-Rule
 Approximate by a cubic polynomial

b
a
3
f ( x )dx   c i f ( x i )  c0 f ( x0 )  c 1 f ( x1 )  c 2 f ( x 2 )  c 3 f ( x 3 )
i 0

3h
 f ( x0 )  3 f ( x 1 )  3 f ( x 2 )  f ( x 3 ) 
8
L(x)
x0
h
f(x)
x1
h
x2
h
x3
x
Simpson’s 3/8-Rule
L( x ) 

( x  x1 )( x  x 2 )( x  x 3 )
( x  x0 )( x  x 2 )( x  x 3 )
f ( x0 ) 
f ( x1 )
( x0  x1 )( x0  x 2 )( x0  x 3 )
( x1  x0 )( x1  x 2 )( x1  x 3 )
( x  x0 )( x  x1 )( x  x 3 )
( x  x0 )( x  x1 )( x  x 2 )
f ( x2 ) 
f ( x3 )
( x 2  x0 )( x 2  x1 )( x 2  x 3 )
( x 3  x0 )( x 3  x1 )( x 3  x 2 )

b
a
f(x)dx  
b
a
ba
L(x)dx ; h 
3
3h
 f ( x0 )  3 f ( x 1 )  3 f ( x 2 )  f ( x 3 )

8
 Truncation error
3 5 (4)
( b  a) 5 ( 4 )
ba
Et  
h f ( )  
f ( ) ; h 
80
6480
3
Example: Simpson’s Rules
 Evaluate the integral

Simpson’s 1/3-Rule
4
I   xe 2 x dx 
0

4
0
xe 2 x dx
h
 f ( 0 )  4 f ( 2 )  f ( 4 )
3


2
0  4 ( 2 e 4 )  4 e 8  8240.411
3
5216.926  8240.411

 57.96%
5216.926


Simpson’s 3/8-Rule
4
I   xe 2 x dx 
0
3h 
4
8

f
(
0
)

3
f
(
)

3
f
(
)

f
(
4
)

8 
3
3

3( 4/3 )
0  3( 19.18922)  3( 552.33933)  11923.832  6819.209
8
5216.926  6819.209

 30.71%
5216.926

Composite Simpson’s 1/3 Rule
function I = Simp(f, a, b, n)
% integral of f using composite Simpson rule
% n must be even
h = (b - a)/n;
S = feval(f,a);
for i = 1 : 2 : n-1
x(i) = a + h*i;
S = S + 4*feval(f, x(i));
end
for i = 2 : 2 : n-2
x(i) = a + h*i;
S = S + 2*feval(f, x(i));
end
S = S + feval(f, b); I = h*S/3;
Simpson’s 1/3 Rule
Composite Simpson’s 1/3 Rule
»
»
»
»
»
»
»
»
»
»
»
»
»
I
»
I
»
I
»
I
»
Q
x=0:0.04:4; y=example(x);
x1=0:2:4; y1=example(x1);
c=Lagrange_coef(x1,y1); p1=Lagrange_eval(x,x1,c);
H=plot(x,y,x1,y1,'r*',x,p1,'r');
xlabel('x'); ylabel('y'); title('f(x) = x*exp(2x)');
set(H,'LineWidth',3,'MarkerSize',12);
x2=0:1:4; y2=example(x2);
c=Lagrange_coef(x2,y2); p2=Lagrange_eval(x,x2,c);
H=plot(x,y,x2,y2,'r*',x,p2,'r');
xlabel('x'); ylabel('y'); title('f(x) = x*exp(2x)');
set(H,'LineWidth',3,'MarkerSize',12);
I=Simp('example',0,4,2)
=
8.2404e+003
I=Simp('example',0,4,4)
=
5.6710e+003
I=Simp('example',0,4,8)
=
5.2568e+003
I=Simp('example',0,4,16)
=
5.2197e+003
Q=Quad8('example',0,4)
=
5.2169e+003
n=2
n=4
n=8
n = 16
MATLAB fun
Multiple applications of Simpson’s rule
with odd number of intervals
Hybrid Simpson’s
1/3 & 3/8 rules
Newton-Cotes Closed
Integration Formulae
ba
h
n
n Name
1
2
3
4
5
Formula
TruncationError
f ( x0 )  f ( x 1 )
1 3
Trapezoidal rule
( b  a)

h f ( )
2
12
f ( x0 )  4 f ( x 1 )  f ( x 2 )
1 5 (4)
Simpson' s 1/3 rule ( b  a )

h f ( )
6
90
f ( x0 )  3 f ( x 1 )  3 f ( x 2 )  f ( x 3 )
3 5 (4)
Simpson's 3/8 rule ( b  a )

h f ( )
8
80
7 f ( x0 )  32 f ( x1 )  12 f ( x 2 )  32 f ( x 3 )  7 f ( x 4 )
8 7 (6 )
Boole' s rule
( b  a)

h f ( )
90
945
19 f ( x0 )  75 f ( x1 )  50 f ( x 2 )  50 f ( x 3 )  75 f ( x 4 )  19 f ( x5 )
275 7 ( 6 )
( b  a)

h f ( )
288
12096
Composite Trapezoidal Rule with
Unequal Segments
4
 Evaluate the integral I  0 xe 2 x dx
 h1 = 2, h2 = 1, h3 = 0.5, h4 = 0.5
2
3
3.5
2
3
I   f ( x )dx   f ( x )dx  
0
f ( x )dx  
4
3.5
f ( x )dx
h1
 f (0 )  f ( 2 )  h2  f ( 2 )  f ( 3)
2
2
h3
h4
  f ( 3 )  f ( 3.5 )   f ( 3.5 )  f ( 4 )
2
2
0.5 6
1
2
3e  3.5 e 7
 0  2 e 4  2 e 4  3e 6 
2
2
2
0.5
   14.45%
3.5 e 7  4 e 8  5971.58

2









Trapezoidal Rule for Unequally Spaced Data
MATLAB Function: trapz
 Z = trapz(x,y)
» x=[0 1 1.5 2.0 2.5 3.0 3.3 3.6 3.8 3.9 4.0]
x =
Columns 1 through 7
0
1.0000
Columns 8 through 11
3.6000
3.8000
» y=x.*exp(2.*x)
y =
1.0e+004 *
Columns 1 through 7
0
0.0007
Columns 8 through 11
0.4822
0.7593
» integr = trapz(x,y)
integr =
5.3651e+003
1.5000
2.0000
3.9000
4.0000
0.0030
0.0109
0.9518
1.1924
2.5000
3.0000
3.3000
0.0371
0.1210
0.2426
Integral of Unevenly-Spaced Data
 Trapezoidal rule
 Could also be evaluated with Simpson’s rule for higher accuracy
Composite Simpson’s Rule with
Unequal Segments
4
• Evaluate the integral I  0 xe 2 x dx
• h1 = 1.5, h2 = 0.5
3
4
I   f ( x )dx   f ( x )dx
0
3
h1
  f (0 )  4 f ( 1.5 )  f ( 3 )
3
h2
  f ( 3 )  4 f ( 3.5 )  f ( 4 )
3
1.5
0.5
3
6

0  4 ( 1.5 e )  3e 
3e 6  4 ( 3.5 e 7 )  4 e 8
3
3
 5413.23    3.76%




Newton-Cotes Open Formula
Midpoint Rule (One-point)

b
a
f ( x )dx  ( b  a) f ( xm )
ab
( b  a) 3
 ( b  a) f (
)
f ( )
2
24
f(x)
a
xm
b
x
Two-point Newton-Cotes Open Formula
 Approximate by a straight line

b
a
ba
( b  a) 3
 f ( x 1 )  f ( x 2 ) 
f ( x )dx 
f ( )
2
108
f(x)
x0
h
x1
h
x2
h
x3
x
Three-point Newton-Cotes Open Formula
 Approximate by a parabola

b
a
ba
2 f ( x1 )  f ( x2 )  2 f ( x3 )
f ( x )dx 
3
7 (b  a)5

f ( )
23040
f(x)
x0
h
x1
h
x2
h
x3
h
x4
x
Newton-Cotes Open
Integration Formulae
ba
h
n
n Formula
TruncationError
1 3
2 ( b  a) f ( x1 )
h f ( )
3
f ( x1 )  f ( x 2 )
3 3
3 (b  a)
h f ( )
2
4
2 f ( x1 )  f ( x 2 )  2 f ( x 3 )
14 5 ( 4 )
4 (b  a)
h f ( )
3
45
11 f ( x 1 )  f ( x 2 )  f ( x 3 )  11 f ( x 4 )
95 5 ( 4 )
5 (b  a)
h f ( )
24
144
11 f ( x 1 )  14 f ( x 2 )  26 f ( x 3 )  14 f ( x 4 )  11 f ( x5 ) 41 7 ( 6 )
6 (b  a)
h f ( )
20
140
Double Integral
 Area under the function surface

d
c
 b f ( x , y )dx dy 
 a


b
a
 d f ( x , y )dy dx 
 c

d
b
c
a

f ( x , y )dxdy
Double Integral
 T(x, y) = 2xy + 2x – x2 – 2y2 + 40
 Two-segment trapezoidal rule
 Exact if using single-segment Simpson’s 1/3 rule
(because the function is quadratic in x and y)