Numerical Integration Formulas

Download as ppt, pdf, or txt
Download as ppt, pdf, or txt
You are on page 1of 60

Chapter 17

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

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


Waterexertingpressureontheupstreamfaceofadam:
(a)sideviewshowingforceincreasinglinearlywith
depth;
(b)frontviewshowingwidthofdaminmeters.

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

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

Approximate f(x) by a polynomial

f n ( x ) a0 a1 x an1 x n1 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


Simpsons 1/3-Rule : Quadratic
Simpsons 3/8-Rule : Cubic
Booles 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

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

Better Numerical Integration


Composite integration
Multiple applications of Newton-Cotes
formulae
Composite Trapezoidal Rule
Composite Simpsons Rule
Richardson Extrapolation
Romberg integration

Apply trapezoidal rule to multiple


segments over integration limits
Two segments

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

Composite Trapezoidal Rule

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)

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

Composite Trapezoidal Rule

x sin( 2 x )dx

function f = example1(x)
% a = 0, b = pi
f=x.^2.*sin(2*x);

Composite Trapezoidal Rule

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

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

Exact = -4. 9348

x 2 sin( 2 x )dx
n=8
I = -4.6785

Exact = -4. 9348

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
h
n 8, h 0.5 I f (0 ) 2 f (0.5 ) 2 f ( 1)
2
2 f ( 1.5 ) 2 f ( 2 ) 2 f ( 2.5 ) 2 f ( 3 )
n 1, h 4 I

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%

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

Simpsons 1/3-Rule

Approximate the function by a parabola

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

Composite Simpsons Rule


Piecewise Quadratic approximations
ba
h
n
f(x)

...
x0 h x1 h x2 h x3 h

x4

xn-2 xn-1

xn

Composite Simpsons 1/3 Rule

Applicable only if the number of segments is even

Composite Simpsons 1/3 Rule


Applicable only if the number of segments is even
I

x2

x0

x4

xn

x2

xn 2

f ( x )dx f ( x )dx

f ( x )dx

Substitute Simpsons 1/3 rule for each integral


f ( x0 ) 4 f ( x 1 ) f ( x 2 )
f ( x 2 ) 4 f ( x 3 ) 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

Simpsons 1/3 Rule


Truncation error (single application)
1 5 (4)
(b a)5 ( 4 )
ba
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 ) (4)
Ea
f
4
180n
5

Composite Simpsons 1/3 Rule


4

Evaluate the integral I xe 2 x dx


0
n = 2, h = 2
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%

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

Example: Simpsons Rules


Evaluate the integral

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

Composite Simpsons 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;

Simpsons 1/3 Rule

Composite Simpsons 1/3 Rule

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 Simpsons rule


with odd number of intervals

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

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
I

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

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
1.5000
2.0000
2.5000
Columns 8 through 11
3.6000
3.8000
3.9000
4.0000
y=x.*exp(2.*x)
y =
1.0e+004 *
Columns 1 through 7
0
0.0007
0.0030
0.0109
0.0371
Columns 8 through 11
0.4822
0.7593
0.9518
1.1924
integr = trapz(x,y)
integr =
5.3651e+003

3.0000

3.3000

0.1210

0.2426

Integral of Unevenly-Spaced Data

Trapezoidal rule
Could also be evaluated with Simpsons rule for higher accuracy

Composite Simpsons Rule with


Unequal Segments
4

Evaluate the integral I xe 2 x dx


0
h1 = 1.5, h2 = 0.5
3

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%

Newton-Cotes Open Formula


Midpoint Rule (One-point)

f ( x )dx ( b a ) f ( x m )

ab
(b a) 3
(b a) f (
)
f ( )
2
24
f(x)

xm

Two-point Newton-Cotes Open Formula


Approximate by a straight line

ba
(b a) 3
f ( x 1 ) f ( x 2 )
f ( x )dx
f ( )
2
108

f(x)

x0

x1

x2

x3

Three-point Newton-Cotes Open Formula


Approximate by a parabola

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

Exact if using single-segment Simpsons 1/3 rule


(because the function is quadratic in x and y)

You might also like