Numerical Integration Formulas
Numerical Integration Formulas
Numerical Integration Formulas
Numerical
Integration Formulas
y f(x)
Integration
I
lim f ( x )x
max x 0 i 1
f ( x)dx
A f ( xi ) xi I
i 1
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
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
f ( x )dx c i f ( x i )
i 0
f(x)
x0
c0 f ( x0 ) c 1 f ( x 1 ) c n f ( x n )
x1
xn-1
xn
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
0
3
11
13
15
Numerical integration
Newton-Cotes formulas
- based on idea
f ( x )dx
f n ( x )dx
f n ( x ) a0 a1 x an1 x n1 an x n
Numerical Integration
Newton-Cotes Closed Formulae -- Use
both end points
midpoint rule
Higher-order methods
(b) Extrapolation
Trapezoidal Rule
Straight-line approximation
f ( x )dx c i f ( x i ) c0 f ( x0 ) c 1 f ( x 1 )
i 0
h
f ( x0 ) f ( x 1 )
2
f(x)
L(x)
x0
x1
Trapezoidal Rule
Lagrange interpolation
L( x )
x x0
x x1
f ( x0 )
f ( x1 )
x0 x 1
x 1 x0
xa
dx
let a x0 , b x 1 ,
, d
; hba
ba
h
x a 0
L( ) ( 1 ) f ( a ) ( ) f ( b )
x b 1
f ( x )dx L( x )dx h L( )d
1
f ( a ) h ( 1 )d f ( b ) h d
1
f ( a ) h( ) f ( b ) h
2 0
2
2
h
f ( a ) f ( b )
2
Example:Trapezoidal Rule
Evaluate the integral
Exact solution
x 2x 1 2x
xe dx e e
4
2
xe 2 x dx
4
2x
1 2x
e ( 2 x 1) 5216.926477
4
0
Trapezoidal Rule
4
4 0
2x
f (0 )
I xe dx
0
f ( 4 ) 2(0 4 e 8 ) 23847.66
2
5216.926 23847.66
357.12%
5216.926
Three segments
11
13
15
Four segments
11
13
15
Many segments
0
3
11
13
15
11
13
15
Multiple Applications of
Trapezoidal Rule
x1
x2
xn
x0
x1
xn 1
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)
ba
h
n
x0
x1
x2
x3
x4
Trapezoidal Rule
Truncation error (single application)
1
Et
f ( )( b a ) 3
12
Exact if the function is linear ( f = 0)
Use multiple applications to reduce the
truncation error
(b a) 3
Ea
12 n 3
f ( )
i1
(b a) 3
f ;
2
12 n
1 n
f f ( i )
n i 1
Approximate
error
x sin( 2 x )dx
function f = example1(x)
% a = 0, b = pi
f=x.^2.*sin(2*x);
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
x 2 sin( 2 x )dx
n=2
I = -1.4239 e-15
Exact = -4. 9348
x 2 sin( 2 x )dx
n=4
I = -3.8758
x 2 sin( 2 x )dx
n=8
I = -4.6785
x 2 sin( 2 x )dx
n = 16
I = -4.8712
2 f ( 3.5 ) f ( 4 ) 5764.76
h
n 16 , h 0.25 I f (0 ) 2 f (0.25 ) 2 f (0.5 )
2
2 f ( 3.5 ) 2 f ( 3.75 ) f ( 4 ) 5355.95
357.12%
132.75%
39.71%
10.50%
2.66%
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
I xe 2 x dx
0
Simpsons 1/3-Rule
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
x1
x2
Simpsons 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 x1 )
ab
2
x x1
ba
dx
h
,
, d
2
h
h
x x0 1
x x1 0
x x 1
2
let
x 0 a, x 2 b, x 1
L( )
( 1)
( 1)
f ( x0 ) ( 1 2 ) f ( x 1 )
f ( x2 )
2
2
Simpsons 1/3-Rule
L( )
( 1)
( 1)
f ( x0 ) ( 1 2 ) f ( x 1 )
f ( x2 )
2
2
h 1
f ( x )dx h L( )d f ( x0 ) ( 1)d
1
2 1
1
h 1
2
f ( x 1 ) h ( 1 )d f ( x 2 ) ( 1)d
0
2 1
f ( x0 )
( ) f ( x 1 ) h( )
2 3
2 1
3 1
3
f ( x2 ) ( )
2 3
2 1
h
f ( x )dx f ( x0 ) 4 f ( x 1 ) f ( x 2 )
3
...
x0 h x1 h x2 h x3 h
x4
xn-2 xn-1
xn
x2
x0
x4
xn
x2
xn 2
f ( x )dx f ( x )dx
f ( x )dx
I 2h
(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
(b a ) (4)
Ea
f
4
180n
5
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%
Simpsons 3/8-Rule
Approximate by a cubic polynomial
f ( x )dx c i f ( x i ) c0 f ( x0 ) c 1 f ( x 1 ) 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
f(x)
x1
x2
x3
Simpsons 3/8-Rule
L( x )
( x x 1 )( x x 2 )( x x 3 )
( x x0 )( x x 2 )( x x 3 )
f ( x0 )
f ( x1 )
( x0 x 1 )( x0 x 2 )( x0 x 3 )
( x 1 x0 )( x 1 x 2 )( x 1 x 3 )
( x x0 )( x x 1 )( x x 3 )
( x x0 )( x x 1 )( x x 2 )
f ( x2 )
f ( x3 )
( x 2 x0 )( x 2 x 1 )( x 2 x 3 )
( x 3 x0 )( x 3 x 1 )( x 3 x 2 )
f(x)dx
ba
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 )
ba
Et
h f ( )
f ( ) ; h
80
6480
3
Simpsons 1/3-Rule
I
xe 2 x dx
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
Simpsons 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
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
Hybrid Simpsons
1/3 & 3/8 rules
Newton-Cotes Closed
Integration Formulae
ba
h
n
n Name
1
2
3
4
5
Formula
Truncation Error
f ( x0 ) f ( x 1 )
1 3
Trapezoida l 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 ( x 1 ) 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 ( x 1 ) 50 f ( x 2 ) 50 f ( x 3 ) 75 f ( x 4 ) 19 f ( x5 )
275 7 ( 6 )
( b a)
h f ( )
288
12096
3.5
f ( x )dx f ( x )dx
f ( x )dx
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
2
1
0.5
0 2 e 4 2 e 4 3e 6
3e 6 3.5 e 7
2
2
2
0.5
3.5 e 7 4 e 8 5971.58
14.45%
2
3.0000
3.3000
0.1210
0.2426
Trapezoidal rule
Could also be evaluated with Simpsons rule for higher accuracy
I f ( x )dx f ( x )dx
0
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
3 e 6 4 ( 3 .5 e 7 ) 4 e 8
3
3
5413.23 3.76%
f ( x )dx ( b a ) f ( x m )
ab
(b a) 3
(b a) f (
)
f ( )
2
24
f(x)
xm
ba
(b a) 3
f ( x 1 ) f ( x 2 )
f ( x )dx
f ( )
2
108
f(x)
x0
x1
x2
x3
ba
2 f ( x 1 ) f ( x 2 ) 2 f ( x 3 )
f ( x )dx
3
7 ( b a) 5
f ( )
23040
f(x)
x0
x1
x2
x3
x4
Newton-Cotes Open
Integration Formulae
ba
h
n
n Formula
Truncation Error
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 ( x1 ) f ( x 2 ) f ( x 3 ) 11 f ( x 4 )
95 5 ( 4 )
5 ( b a)
h f ( )
24
144
11 f ( x1 ) 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
Areaunderthefunctionsurface
d
b
f ( x , y )dx dy
a
f ( x , y )dy dx
c
f ( x , y )dxdy
Double Integral
T(x, y) = 2xy + 2x x2 2y2 + 40
Two-segment trapezoidal rule