Transcript chapter11

11. Numerical Differentiation and Integration
11.3 Better Numerical Integration, 11.4 Gaussian Quadrature, 11.5 MATLAB’s Methods
Natural Language Processing Lab
Dept. of Computer Science and Engineering, Korea Univertity
CHOI Won-Jong ([email protected])
Woo Yeon-Moon([email protected])
Kang Nam-Hee([email protected])
Contents



11.3 BETTER NUMERICAL INTEGRATION
 11.3.1 Composite Trapezoid Rule
 11.3.2 Composite Simpson’s Rule
 11.3.3 Extrapolation Methods for Quadrature
11.4 GAUSSIAN QUADRATURE
 11.4.1 Gaussian Quadrature on [-1, 1]
 11.4.2 Gaussian Quadrature on [a, b]
11.5 MATLAB’s Methods
2
11.3.1 Composite Trapezoid
Rule
11.3.2 Composite Simpson’s
Rule
CHOI WonJong ([email protected])
3
11.3.1 Composite Trapezoid
Rule
CHOI WonJong ([email protected])
4
11.3 BETTER NUMERICAL INTEGRATION

Composite integration(복합적분) : Applying one of the
lower order methods presented in the previous section
repeatedly on several sub intervals.
5
11.3.1 Composite Trapezoid Rule
If
we divide the interval of integration, [a, b], into two or more
subintervals and use the trapezoid rule on each subintervals, we
obtain the composite trapezoid rule.
ba
h
2
h
h
f
(
x
)
dx

f
(
x
)
dx

f
(
x
)
dx

[
f
(
a
)

f
(
x
)]

[ f ( x1 )  f (b)]
1
a
a
x1
2
2
h
ba
 [ f (a)  2 f ( x1 )  f (b)] 
[ f (a)  2 f ( x1 )  f (b)]
2
4
b
x1
b
6
11.3.1 Composite Trapezoid Rule
If
we divide the interval into n subintervals, we get

b
a
x1
f ( x)dx   f ( x )dx 
a

b
xn1
f ( x )dx
h
ba
n
h
h
 [ f (a )  f ( x1 )]   [ f ( xn 1 )  f (b)]
2
2
ba

[ f (a )  2 f ( x1 )   2 f ( xn 1 )  f (b)]
2n
MATLAB CODE
7
11.3.1 Composite Trapezoid Rule
Example
11.9
n=1
n=2
n=4
n=20
n=3
n=100
8
11.3.1 Composite Trapezoid Rule
Example
11.9
1
b
dx

[log
|
x
|

C
]
a
a x
21
2
1 x dx  [log | 2 | C ]  [log |1| C ]  log 1  0.69314718055995
b
9
11.3.2 Composite Simpson’s
Rule
CHOI WonJong ([email protected])
10
11.3.2 Composite Simpson’s Rule
Example
11.10
11
11.3.2 Composite Simpson’s Rule


Applying the same idea of subdivision of intervals to
Simpson’s rule and requiring that n be even gives the
composite Simpson rule.
[a,b]를 two subintervals [a,x2], [x2, b]로 나눈다면,
ba
ba
x2 
,h 
2
4

b
a
x2
b
a
x2
f ( x)dx   f ( x)dx   f ( x)dx
h
h
 [ f (a)  4 f ( x1 )  f ( x2 )]  [ f ( x2 )  4 f ( x3 )  f (b)]
3
3
12
11.3.2 Composite Simpson’s Rule
In
general, for n even, we have h=(b-a)/n, and Simpson’s rule is
ba
h
n

b
a
h
f ( x)dx  [ f (a)  4 f ( x1 )  2 f ( x2 )  4 f ( x3 )  2 f ( x4 )   2 f ( xn2 )  4 f ( xn1 )  f (b)]
3
13
11.3.2 Composite Simpson’s Rule
Example
11.10
14
11.3.2 Composite Simpson’s Rule
Example
11.11 Length of Elliptical Orbit
3
x(r )  cos(r ), y (r )  sin(r )
4
L
b
a
( x ')  ( y ') dr  0.25
2
2
b
a
16sin 2 ( r )  9 cos 2 ( r ) dr
15
11.3.2 Composite Simpson’s Rule
Example
11.11 Length of Elliptical Orbit
3
x(r )  cos(r ), y (r )  sin(r )
4
L
b
a
( x ') 2  ( y ') 2 dr  0.25
b
a
16sin 2 ( r )  9 cos 2 ( r ) dr
days
0 10 20 30 40 50 60 70 80 90 100
r = [0.00 1.07 1.75 2.27 2.72 3.14 3.56 4.01 4.53 5.22 6.28]
Using
Composite Simpson’s Rule and the length between day 0
and 10 (n=20) is 0.88952. (Trapezoid=0.889567, Text=0.8556)
Using Composite Simpson’s Rule and the length between day 60
and 70 (n=20) is 0.382108. (Trapezoid=0.382109, Text=0.3702)
The former is 2.3279 times faster than the latter.
16
11.3.3 Extrapolation Methods
for Quadrature
Woo Yeon-Moon([email protected])
17
Richardson Expolation
Truncation error(절단 오차)

•
b

a
I  I ( f , h)  E ( f , h)
사다리꼴

h
2j
f ( x)dx  [ f (a)  2 f ( x1 )  ...  2 f ( xn 1 )  f (b)]   c j h
2
j 1
h
 [ f (a)  2 f ( x1 )  ...  2 f ( xn 1 )  f (b)]  ch 2
simpso
2
n
h
h
4
 [ f (a)  4 f ( x1 )  2 f ( x2 )]  [ f ( x2 )  4 f ( x3 )  f (b)]  ch
3
3
18
Richardson Expolation

To obtain an estimate that is more accurate
• using two or more subintervals (h를 줄임)
- 그러나, 세부구간의 수가 일정한 범위를 넘어서면
round-off error가 커지게 된다.
계
산
오
차
trapezoid
simpson
세부 구간의 수
 Richardson Extrapolation
간격이 다른 2개의 식을 구한 결과를 대수적으로 정리함으로써
보다 정확한 값을 산출
19
Richardson Extrapolation

Richardson Extrapolation using the trapezoid rule
I  IT (h1 )  ch  IT (h2 )  ch2
2
1
2
I (h2 )  I (h1 )
I  I (h2 ) 
 h 2 
1
   1
 h2 

(if h_2 = ½ h_1)
4 I (h2 )  I (h1 )
I
3
Simpson
rules
20
Example 11.12 Integral of 1/x
 start with one subinterval (h=1)
2
dx
1
1 1 1 3
1 x  I0  2 [ f (1)  f (2)]  2 [1  2 ]  4  0.75
 two subintervals (h=1/2)
1
1 1 2 1 17
I1  [ f (1)  2 f (1.5)  f (2)]  [ 
 ]
 0.7083
4
4 1 1.5 2 24
 to apply Richardson extrapolation
1
h
A  [4 A( )  A(h)]
3
2
1
I  [4 I1  I 0 ]  [4(0.7083)  0.7500] / 3  0.6944
3
 exact value of the integral is ln(2)=0.693147..
21
Example 11.12 Integral of 1/x

Form a table of the approximations
Ⅰ
h=1
0.7500
h=1/2

Ⅱ
0.7083
0.6944
E0  0.75  0.6944  0.0556
E1  0.7083  0.6944  0.0139
E0  (2) 2 E1
E (h)  O(h )
2

0.6944 ≠0.693147
22
Romberg Integration

Approximate an Error
O(h 2 )
 Trapezoid rules :
 Richardson extrapolation :

O(h 4 )
continued ( using simpson rules)
I  I S (h1 )  ch14  I S (h2 )  ch2 4
I (h2 )  I (h1 )
I  I (h2 ) 
 h 4 
1
   1
 h2 

16 I (h2 )  I (h1 )
I
15
23
Romberg Integration

Improving the result by Richardson extrapolation
E  c1h  c2 h  c3h  c4 h  c5h 
2

1st
4
2nd
6
8
10
3rd
Romberg integration : iterative procedure using
Richardson extrapolation
4k I (h / 2)  I (h)
I ( h) 
k
4 1

k means the improving level(=
degr ee of t he er r or
)
2
24
Example 11.12 Integral of 1/x using
Romberg Integration

b

a
Trapezoid rule
2
1
h
f ( x)dx   dx  [ f (a)  2 f ( x1 )  ...  2 f ( xn 1 )  f (b)]
x
2
1
 For k=0, I_0 = 0.75
 For k=1, I_1 = 0.7083
 For k=2, I_2 = 0.6941

To apply Richardson extrapolation A(h)  1  4 A( h )  A(h) 
Ⅰ
h=1
0.7500
h=1/2
0.7083
h=1/4
0.6970
h=1/8
0.6941
Ⅱ
3
2

0.6944
0.6933
0.6943
25
Example 11.12 Integral of 1/x using
Romberg Integration
second
level of extrapolation
1 
h

C (h)  16 B( )  B(h) 
15 
2

Ⅰ
h=1
0.7500
h=1/2
0.7083
h=1/4
0.6970
Ⅱ
Ⅲ
0.6944
0.6933
[16(0.6933)-0.6944]/15
26
Example 11.12 Integral of 1/x using
Romberg Integration
2

five levels of extrapolation to find values for
0.7500
0.6944
0.6932
0.6931
0.6931
0.7083
0.6933
0.6931
0.6931
0.6931
0.6970
0.6932
0.6931
0.6931
0.6941
0.6931
0.6931
0.6934
0.6931
1
1 x dx
0.6931
0.6932
27
Matlab function for Romberg Integration
28
11.4 Gaussian Quadrature
Kang Nam-Hee ([email protected])
29
11.4.1 Gaussian Quadrature on [-1,1]

Gaussian Quadrature Formular
 Get the definite integration of f(x) on [-1,1] using linear
combinations of coefficient ck and evaluated function
value f(xk) at the point xk
 Appropriate values of the points xk and ck depend on the
choice of n
 By choosing the quadrature point x1 ,… xn as the n zeros
of the nth-degree Gauss-Legendre polynomial, and by
using the appropriate coefficients, the integration
formular is exact for polynomials of degree up to 2n-1
30
11.4.1 Gaussian Quadrature on [-1,1]

Gaussian Quadrature Formular (cont.)
 n=2
 n=3
31
11.4.1 Gaussian Quadrature on [-1,1]

Example 11.13 integral of exp(-x2) Using G.Q

n
Xi
ci
2
3
±0.557753
0
±0.77459
±0.861136
±0.339981
1
8/9
5/9
0.34785
0.652145
4
Table 11.2 parameters of Gaussian quadrature
32
11.4.1 Gaussian Quadrature on [-1,1]

Gaussian-Legendre Polynomials
33
11.4.2 Gaussian Quadrature on [a,b]

Extends Gaussian Quadrature for f(t) on [a, b]
 by Transformation f(t) on [a, b] to f(x) on [-1,1]
 For the given integral
 change interval of t by using next formular
 so the interval
34
11.4.2 Gaussian Quadrature on [a,b]

Extends Gaussian Quadrature for f(t) on [a, b] (cont.)

f(t) rewrite for variable x
remark the factor (b-a)/2 (∵td convert to dx)
 Apply f(x) to the integral
35
11.4.2 Gaussian Quadrature on [a,b]

Example 11.14
integral of exp(-x2) on [0,2] using G.Q with n = 2
 Consider again the integral
 Transform f(t) on [0,2] to f(x) on [-1,1] using next formular
36
11.4.2 Gaussian Quadrature on [a,b]

Example 11.14 (cont)
 So we can get
 Apply Gaussian Quadrature to the integral with n = 2
37
11.4.2 Gaussian Quadrature on [a,b]

Matlab function for Gaussian Quadrature
38
11.5 MATLAB’s Methods
Woo Yeon-Moon([email protected])
39
11.5 MATLAB’s Methods






p=polyfit(x,y,n) – find the coefficients of the
polynomial of degree n
polyder(p) - calculates the derivative of
polynomials
diff(x) - x = [1 2 3 4 5];
y = diff(x)
y= 1 1 1 1
traps(x,y)
Q=quad(‘f’,xmin,xmax) (simpson rules)
Q=quad8(‘f’,xmin,xmax) (Newton-Cotes eight-panel
rule)
40