第10,11回 数値積分法

第10,11回 数値積分法
1.基本的事項
 定積分の復習
b
関数 f(x)
I = ∫ f ( x)dx = F (b) − F (a )
a
ここで F ( x) は f ( x) の原始関数.
すなわち F ′( x) = f ( x)
x
a
b
⎛π ⎞
⎛ π⎞
cos xdx = [sin x ]π−π/ 2/ 2 = sin ⎜ ⎟ − sin ⎜ − ⎟ = 1 − (−1) = 2
−π / 2
⎝2⎠
⎝ 2⎠
例I = ∫
π /2
·  原始関数が簡単な形式では書けない被積分関数が存在する
(原始関数が求まらない)
·  解析的には解けない積分が存在する
2
1.基本的事項
 数値積分とは?
·  被積分関数の代わりに補間多項式を用いる
b
∫a
f ( x)dx
被積分関数
b
∫a
補間多項式 Pn(x)
Pn ( x ) dx
補間多項式
(多項式は必ず
積分できる)
被積分関数 f(x)
x
a
b
 数値積分公式の種類
·  中点公式
·  台形公式
·  シンプソン公式(1/3,3/8)
3
1.基本的事項
 補間多項式を用いた数値積分
f1
f0
f3
f4
x0
a
x1
fn
Pn で補間
f2
関数 f
x2
x3
f6
f5
x4
x5
x6
...
xn
b
x
b
積分値 I = ∫ f ( x)dx
a
b
近似値 I n = ∫a Pn ( x)
n が大きくなってしまう
4
1.基本的事項
 補間多項式を用いた数値積分
fn
P3 で補間
f1 f
2
f3
f0
P3 で補間
f4
x0
x1
x2
関数 f
x3
f6
f5
x4
x5
x6
...
xn
x
I 3(1)
I 3(0)
b
積分値 I = ∫ f ( x)dx
a
( 0)
(1)
(k )
近似値 I 3 = I 3 + I 3 + … + I 3
複合公式:多項式の次数をあまり大きくしないですむ
5
2.ニュートン・コーツ公式
 等間隔分点 xi+1 – xi = h に対し,ニュートンの補間多項式を使う
(ラグランジュ補間でも同じ) Δf ( x0 )
Δ2 f ( x0 )
PN ( x) = f ( x0 ) + ( x − x0 )
+ ( x − x0 )( x − x1 )
+
2
h
2!h
+ ( x − x0 )( x − x1 )  ( x − x N −1 )
Δf ( xk ) = f ( xk +1 ) − f ( xk ),
ΔN f ( x0 )
n!h N
Δr f ( xk ) = Δr −1 f ( xk +1 ) − Δr −1 f ( xk ).
N=0
N=1
N=2
N=3
開型公式
閉型公式
6
2.ニュートン・コーツ公式
2.1 N = 0:中点公式(開型公式)
fn
fi
f0
N=0
x0
xi–h/2
xi
xi+h/2
xn
ひとつの区間の補間多項式P0(i ) = f ( xi ) = f i
(i )
ひとつの区間の積分 I 0
xi + h / 2 (i )
P ( x)dx
xi − h / 2 0
=∫
=∫
xi + h / 2
f dx
xi − h / 2 i
= h ⋅ fi
全体の積分 I 0 = I 0(0) + I 0(1) + … + I 0( n )
= h( f 0 + f1 + … + f n )
7
2.ニュートン・コーツ公式
2.2 N = 1:台形公式(閉型公式)
fi+1
fi
fn
f0
xi h xi+1
xn
Δf ( xi )
f −f
(i )
= f i + ( x − xi ) i +1 i
ひとつの区間の補間多項式 P1 = f i + ( x − xi )
h
h
xi +1 ( i )
xi +1 ⎛
f i +1 − f i ⎞
(i )
ひとつの区間の積分 I1 = ∫ P1 ( x)dx = ∫ ⎜ f i + ( x − xi )
⎟dx
xi
xi ⎝
h
⎠
N=1
x0
x − xi
,
→
h
dx = h ⋅ du
u≡
1
= ∫ h( f i + u ( f i +1 − f i ) )du
0
1
⎡
⎤
u
h
= h ⎢uf i + ( f i +1 − f i )⎥ = ( f i + f i +1 )
2
⎣
⎦0 2
2
8
2.ニュートン・コーツ公式
2.2 N = 1:台形公式(閉型公式)
fi+1
fi
fn
f0
N=1
x0
xi
xi+1
xn
h
ひとつの区間の積分 I1(i ) = ( f i + f i +1 )
2
( 0)
(1)
( 2)
( n −1)
全体の積分 I1 = I1 + I1 + I1 + … + I1
h
= ( f 0 + 2 f1 + 2 f 2 + … + 2 f n −1 + f n )
2
9
演習
1.2
定積分 ∫ ( x 4 + k )dx
(k=1,2,3) を,積分区間 [0, 1.2] を n (4, 6) 等分
0
して,
台形公式により求め,解析解と比較せよ.
f4
4等分の場合
f3
f0
0
f1
f2
1.2
10
2.ニュートン・コーツ公式
2.3 N = 2:シンプソン1/3公式
fi+1
fi
fi+2
fn
f0
N=2
x0
xi
2つの区間の補間多項式
P2(i )
2つの区間の積分
I 2(i )
h
xi+1
h
xi+2
xn
Δf i
Δ2 f i
= f i + ( x − xi )
+ ( x − xi )( x − xi +1 ) 2
h
2h
=∫
xi + 2
xi
P2(i ) ( x)dx
11
2.ニュートン・コーツ公式
Δf ( xi ) = f i +1 − f i
Δ2 f ( xi ) = Δ1 f ( xi +1 ) − Δ1 f ( xi )
= f i + 2 − f i +1 − ( f i +1 − f i )
= f i + 2 − 2 f i +1 + f i
12
2.ニュートン・コーツ公式
fi
fi+1
fi+2
fn
f0
x0
xi
xi+1
xi+2
xn
h
h
h
h
( f 0 + 4 f1 + f 2 ) ( f 2 + 4 f 3 + f 4 ) ( f 4 + 4 f 5 + f 6 ) ( f n − 2 + 4 f n −1 + f n )
3
3
3
3
全体の積分
13
演習
1.2
定積分∫ ( x 4 + k )dx
(k=1,2,3) を,積分区間 [0,1.2] を n (4, 6) 等分
0
して,
シンプソン1/3公式により求め,解析解と比較せよ.
f4
4等分の場合
f3
f0
0
f1
f2
1.2
14
2.ニュートン・コーツ公式
2.4 N = 3:シンプソン3/8公式
fi
fi+1
fi+2
fi+3
fn
f0
N=3
x0
3つの区間の補間多項式
xi
P3(i )
xi+1
xi+2
xi+3
Δf i
Δ2 f i
= f i + ( x − xi )
+ ( x − xi )( x − xi +1 ) 2
h
2h
+ ( x − xi )( x − xi +1 )( x − xi + 2 )
(i )
3つの区間の積分 I 3
xn
=∫
xi + 3
xi
Δ3 f i
2 ⋅ 3h 3
P3(i ) ( x)dx
15
2.ニュートン・コーツ公式
3つの区間の積分
Δ3 f ( xi ) = Δ2 f ( xi +1 ) − Δ2 f ( xi )
全体の積分
= Δ1 f ( xi + 2 ) − Δ1 f ( xi +1 ) − (Δ1 f ( xi +1 ) − Δ1 f ( xi ))
= f i +3 − 2 f i + 2 + f i +1 − ( f i + 2 − 2 f i +1 + f i )
= f i +3 − 3 f i + 2 + 3 f i +1 − f i
16
2.ニュートン・コーツ公式
 N = 0:中点公式(開型公式)
I 0 = h( f 0 + f1 + … + f n )
 N = 1:台形公式(閉型公式)
h
I1 = ( f 0 + 2 f1 + 2 f 2 + … + 2 f n −1 + f n )
2
 N = 2:シンプソン1/3公式
h
I 2 = ( f 0 + 4 f1 + 2 f 2 + 4 f 3 + … + 2 f n − 2 + 4 f n −1 + f n )
3
 N = 3:シンプソン3/8公式
I3 =
3h
( f 0 + 3 f1 + 3 f 2 + 2 f 3 + … + 2 f n −3 + 3 f n − 2 + 3 f n −1 + f n )
8
いずれも次のような形をしている
I N = a0( N ) f 0 + a1( N ) f1 + … + an( N ) f n
17
演習
1.2
定積分∫ ( x 4 + k )dx
(k=1,2,3) を,積分区間 [0,1.2] を n (3, 6) 等分して,
0
シンプソン3/8公式により求め,解析解と比較せよ.
f6
6等分の場合
f0
0
f1
f2
f3
f5
f4
1.2
18
3.誤差の解析
関数 f(x) は十分に滑らかであるとする.ニュートン・コーツ公式の誤差につ
いてシンプソン1/3公式を例に考える.関数 f(x) をテーラー展開して
∞
1
f ( x) = ∑ ( x − x0 ) k f 0( k )
k = 0 k!
区間 [x0, x2] で積分すると
x2
∫x
0
x2 ⎛ ∞
x2
∞
⎞
⎡
⎤
1
1
f ( x)dx = ∫ ⎜⎜ ∑ ( x − x0 ) k f 0( k ) ⎟⎟dx = ⎢ ∑
( x − x0 ) k +1 f 0( k ) ⎥
x0
⎝ k = 0 k!
⎠
⎣k =0 (k + 1) × k!
⎦ x0
2 k +1 k +1 ( k )
=∑
h f0
k = 0 ( k + 1)!
∞
一方,シンプソン1/3公式の[x0, x2]の区間をテーラー展開すると
⎛ ∞ 1 k (k ) ⎞ ⎛ ∞ 2k k (k ) ⎞ ⎞
1
1 ⎛⎜
h( f 0 + 4 f1 + f 2 ) = h⎜ f 0 + 4⎜⎜ ∑ h f 0 ⎟⎟ + ⎜⎜ ∑ h f 0 ⎟⎟ ⎟⎟
3
3 ⎝
⎝ k = 0 k!
⎠ ⎝ k = 0 k!
⎠⎠
∞
2 k + 4 k +1 ( k )
= 2hf 0 + ∑
h f0
k =1 3 × k!
19
3.誤差の解析
両者を比較すると
1
f
(
x
)
dx
−
h( f0 + 4 f1 + f2 )
∫x0
3
∞
∞ k
2k +1 k +1 ( k )
2 + 4 k +1 ( k )
= ∑
h f0 − 2hf0 − ∑
h f0
k = 0 ( k + 1)!
k =1 3 × k !
E0 =
x2
⎛ 2k +1
2k + 4 ⎞ k +1 ( k )
= ∑⎜
−
⎟ h f0
3× k ! ⎠
k =1 ⎝ ( k + 1)!
1 5 (4)
h5 (4)
= − h f0 + K ≤
f (ξ ) ここで x0 ≤ ξ ≤ x2
90
90
全区間にわたる誤差は
∞
E = E0 + E2 + ...
b − a ⎛ 1 5 (4) ⎞
≤
⎜ hu ⎟
2h ⎝ 90
⎠
b − a 4 (4)
=
hu
ここで u = max{| f (4) ( x) |}
a ≤ x ≤b
180
20
3.誤差の解析
 ニュートン・コーツ公式の誤差をまとめると
b − a 2 (2)
hu
12
b − a 4 (4)
n=2
hu
180
b − a 4 (4)
n=3
hu
80
2(b − a ) 6 (6)
n=4
hu
945
55(b − a ) 6 (6)
n=5
hu
12096
n =1
左の表は,n が大きくて偶数の場
合が誤差が小さいことを示して
いる.これはあくまで関数 f(x) が
滑らかな場合であり,そうでない
場合は, n が小さくて奇数の場合
が誤差が小さいこともある.
21
計算機演習
 台形公式,シンプソン1/3公式,シンプソン3/8公式のプログラムを作れ.
 授業中に行った演習問題を上記のプログラムを用いて,台形公式,シ
ンプソン1/3公式,シンプソン3/8公式で解け.
dx
 定積分 ∫−1
を,積分区間 [–1,1] を n (6, 60, 120, 180, 240) 等分
3− x
して,上記のプログラムを用いて,台形公式,シンプソン1/3
公式,シンプソン3/8公式により 求め,解析解 log 2 = 0.69314718...と比較せよ. 1
22