Document

数値計算法I
第8回
数値積分
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
数値積分
• 関数f(x)が複雑な形をしていて、解析的に積分し
にくい時、数値計算により(定)積分を行う。
b
I   f (x)dx
a
n
I   wi f (xi )
i0
aからbの間をn分割する。
f(xi)は各分割点での値
wiは各分割点での値に与える重み
重みの与え方によって、積分方法が異なる。
2
台形公式による数値積分
最も単純で直感的な方法
f (x)
h
x0
a
x1
x2
---- xn-1 xn
b
ba
刻み幅 h 
n
3
台形公式による数値積分
f (x1 ) f (x2 )
x0
a
x1
赤線で囲まれた台形の面積
x2
---- xn-1 xn
b
f (xi )  f (xi 1 )
h
2
で各区間の積分を近似
4
台形公式による数値積分
f (x0 )  f (x1 )
f (x1 )  f (x2 )
f (xn  2 )  f (xn 1 )
f (xn 1 )  f (xn )
h
h L
h
h
2
2
2
2
h
h
 f (x0 )  f (x1 )h  f (x2 )h  L  f (xn  2 )h  f (xn 1 )h  f (xn )
2
2
n 1
h
h
 f0  h fi  fn
2
2
i 1
I
n 1

h
  f0  2 fi  fn 

2
i 1
簡単化のため、f (xi ) を fi と書く
5
台形公式による数値積分の誤差
fi+1をxiの近傍でテイラー展開する
2
3
4
h  h  h (4 )
fi 1  fi  hfi 
fi 
fi 
fi  L
2
3!
4!
xiでの微分は
2
3
fi 1  fi h  h  h (4 )
fi 
 fi 
fi 
fi  L
h
2
6
24
6
台形公式による数値積分の誤差
Ii 

xi1
xi
h
f (x)dx   f (xi  z)dz
0

z 2  z 3  z 4 (4 )
   fi  zfi 
fi  fi 
fi  L
0 
2
3!
4!
h

 dz
h2
h 3  h 4  h 5 (4 )
 hfi 
fi 
fi 
fi 
fi  L
2
6
24
120
fi 1  fi h  h 2  h 3 (4 )
fi 
 fi 
fi 
fi  L を代入すると、
h
2
6
24
fi 1  fi h 3  h 4
Ii  h

fi 
fi L
2
12
24
台形の面積
この分が誤差
7
台形公式による数値積分の誤差
fi 1  fi h 3  h 4
Ii  h

fi 
fi L
2
12
24
一区間辺りの誤差Ei
h 3 
Ei 
fi
12
全部で、
h 3 n 
E   Ei   fi
12 i 1
i 1
n
区間中の2階微分の最大値をMとすると、
h3
h2
b  a h2
ba 2
E  nM  nh M  n
M
h M
12
12
n 12
12
8
台形公式による数値積分の誤差
ba 2
E
h M
12
台形公式の誤差は刻み幅の2乗に比例する。
9
台形公式による数値積分（例）
4
0 1  x 2 dx
1
n
h
積分値
誤差
2
0.500000000000
3.100000000000
-0.041592653590
4
0.250000000000
3.131176470588
-0.010416183002
8
0.125000000000
3.138988494491
-0.002604159099
16
0.062500000000
3.140941612041
-0.000651041548
32
0.031250000000
3.141429893175
-0.000162760415
64
0.015625000000
3.141551963486
-0.000040690104
128
0.007812500000
3.141582481064
-0.000010172526
256
0.003906250000
3.141590110458
-0.000002543132
512
0.001953125000
3.141592017807
-0.000000635783
10
台形公式による数値積分（例）
10-1
10-2
error
10-3
10-4
10-5
10-6
10-7 -3
10
10-2
10-1
100
h
11
台形公式による数値積分例
implicit real*8 (a-h,o-z)
parameter(nmax=1000)
c
pi=4*datan(1.d0)
b=1
a=0
do in=1,9
n=2**in
h=(b-a)/n
sum=0
do i=0,n-1
sum=sum+(f(a+i*h)+f(a+(i+1)*h))*h/2
end do
error=sum-pi
write(6,100)n,h,sum,error
100 format(1x,i5,3(1x,f20.12))
end do
stop
end
c
function f(x)
implicit real*8 (a-h,o-z)
f=4/(1+x*x)
return
end
12
シンプソン法
fi 1
fi
fi 1
h
xi 1
h
xi
xi 1
fi-1,fi,fi+1の3点を通る2次曲線で、xi-1からxi+1
の間のf(x)を近似する。
2次補間
13
2次補間
P2 (x)  y1a1 (x)  y2 a2 (x)  y3a3 (x)
ここで、a1(x),a2(x),a3(x)は2次式
P2(x)が、(x1,y1),(x2,y2),(x3,y3)を通る時、
P2 (x1 )  y1  y1a1 (x1 )  y2 a2 (x1 )  y3a3 (x1 )
より、
a1 (x1 )  1, a2 (x1 )  0, a3 (x1 )  0
同様に、
a1 (x2 )  0, a2 (x2 )  1, a3 (x2 )  0
a1 (x3 )  0, a2 (x3 )  0, a3 (x3 )  1
14
2次補間
a1 (x1 )  1, a2 (x1 )  0, a3 (x1 )  0
a1 (x2 )  0, a2 (x2 )  1, a3 (x2 )  0
a1 (x3 )  0, a2 (x3 )  0, a3 (x3 )  1
を満たす2次式は、
a1 (x)  A(x  x2 )(x  x3 )
a2 (x)  B(x  x1 )(x  x3 )
a3 (x)  C(x  x1 )(x  x2 )
の形をしている。
15
2次補間
a1 (x1 )  1 となるには、係数Aは
A(x1  x2 )(x1  x3 )  1
より、
同様に
1
A
(x1  x2 )(x1  x3 )
1
B
(x2  x1 )(x2  x3 )
1
C
(x3  x1 )(x3  x2 )
16
2次補間
P2 (x)  y1a1 (x)  y2 a2 (x)  y3a3 (x)
は
(x  x2 )(x  x3 )
(x  x1 )(x  x3 )
(x  x1 )(x  x2 )
P2 (x)  y1
 y2
 y3
(x1  x2 )(x1  x3 )
(x2  x1 )(x2  x3 )
(x3  x1 )(x3  x2 )
と書ける。
fi-1,fi,fi+1の3点を通る2次曲線は
P2 (x)  fi 1
(x  xi )(x  xi 1 )
(x  xi 1 )(x  xi 1 )
(x  xi 1 )(x  xi )
 fi
 fi 1
(xi 1  xi )(xi 1  xi 1 )
(xi  xi 1 )(xi  xi 1 )
(xi 1  xi 1 )(xi 1  xi )
と書ける。
17
シンプソン公式
xi-1からxi+1の積分

xi1
xi1
f (x)dx
f(x)をP2(x)で近似して積分する。
fi 1
fi
fi 1
h
xi 1
h
xi
xi 1
18
シンプソン公式

xi1
xi1
f (x)dx 

xi1
xi1
P2 (x)dx
(x  xi )(x  xi 1 )
  fi 1
dx
xi1
(xi 1  xi )(xi 1  xi 1 )
xi1
(x  xi 1 )(x  xi 1 )
  fi
dx
xi1
(xi  xi 1 )(xi  xi 1 )
xi1
(x  xi 1 )(x  xi )
  fi 1
dx
xi1
(xi 1  xi 1 )(xi 1  xi )
xi1
fi 1

(x  xi )(x  xi 1 )dx

x
(xi 1  xi )(xi 1  xi 1 ) i1
xi1
fi

(x  xi 1 )(x  xi 1 )dx

x
(xi  xi 1 )(xi  xi 1 ) i1
xi1
fi 1

(x  xi 1 )(x  xi )dx

x
(xi 1  xi 1 )(xi 1  xi ) i1
xi1
19
シンプソン公式

xi1
xi1
xi1
fi 1
(x  xi )(x  xi 1 )dx
f (x)dx 

x
i1
(xi 1  xi )(xi 1  xi 1 )
xi1
fi

(x  xi 1 )(x  xi 1 )dx

x
(xi  xi 1 )(xi  xi 1 ) i1
xi1
fi 1

(x  xi 1 )(x  xi )dx

x
(xi 1  xi 1 )(xi 1  xi ) i1
xi 1  xi  h, xi 1  xi 1  2h 等を使って、
xi1
fi 1

(x  xi )(x  xi 1 )dx

x
(h)(2h) i1
xi1
fi

(x  xi 1 )(x  xi 1 )dx

x
(h)(h) i1
xi1
fi 1

(x  xi 1 )(x  xi )dx

x
(2h)(h) i1
20
シンプソン公式




xi1
xi1
xi1
xi1
xi1
xi1
xi1
xi1
fi 1 xi1
f (x)dx  2 x (x  xi )(x  xi 1 )dx
2h i1
fi xi1
 2  (x  xi 1 )(x  xi 1 )dx
h xi1
fi 1 xi1
 2  (x  xi 1 )(x  xi )dx
2h xi1
2h
 z3 3 2
2
2 
(x  xi )(x  xi 1 )dx   (z  h)(z  2h)dz    hz  2h z   h 3
3
3 2
0
0
2h
2h
 z3
4
2
(x  xi 1 )(x  xi 1 )dx   (z)(z  2h)dz    hz    h 3
3
3
0
0
2h
(x  xi 1 )(x  xi )dx 
2h
z
h 2
2 3
(z)(z

h)dz


z

h
3 2 
0
3

0
2h
3
21
シンプソン公式

xi1
xi1
f (x)dx 
fi 1 2 3 fi  4 2  fi 1 2 2
h  2  h   2 h
2
2h 3
h  3  2h 3
h
nは偶数
  fi 1  4 fi  fi 1 
3
よって、x0(=a)からxn(=b)までの積分値は、

b
a
4から6の積分
0から2の積分 2から4の積分
h
h
h
f (x)dx   f0  4 f1  f2    f2  4 f3  f4    f4  4 f5  f6 
3
3
3
h
h
L  fn  4  4 fn 3  fn  2    fn  2  4 fn 1  fn 
3
3

b
a
n /2
n /21

h
f (x)dx   f0  4  f2n1  2  f2n  fn 

3
i 1
i 1
22
シンプソン公式の誤差
fi+1とfi-1をxiの近傍でテイラー展開する
h2
fi 1  fi  hfi 
2
h2
fi 1  fi  hfi 
2
3
h
fi 
3!
3
h
fi 
3!
4
5
6
h
h
h
fi 
fi (4 ) 
fi (5) 
fi (6) L
4!
5!
6!
4
5
6
h
h
h
(4 )
(5)
(6)

fi 
fi 
fi 
fi L
4!
5!
6!
和をとって、
4
6
h
h
2 
(4 )
(6)
fi 1  fi 1  2 fi  h fi 
fi 
fi  L
12
360
差をとって、
h 3  h 5 (5)
fi 1  fi 1  2hfi 
fi 
fi
3
60
23
シンプソン公式の誤差

xi1
xi1
f (x)dx 

h
h
f (xi  z)dz 

h
0
0
f (xi  z)dz   f (xi  z)dz
h
h
z 2  z 3  z 4 (4 ) z 5 (5) z 6 (6)
  ( fi  zfi 
fi  fi 
fi 
fi 
fi L )dz
2
3!
4!
5!
6!
0
0
z 2  z 3  z 4 (4 ) z 5 (5) z 6 (6)
  ( fi  zfi 
fi  fi 
fi 
fi 
fi L )dz
2
3!
4!
5!
6!
h
2 h 3  2 h 5 (4 ) 2 h 7 (6)
 2hfi 
fi 
fi 
fi  L
3 2
5 4!
74 6!
6
h (4 ) h
(6)

fi 1  fi 1  2 fi  h fi 
fi 
fi  L
12
360
より、
2
4
f

2
f

f
h
h
(4 )
(6)
i
i 1
fi  i 1

f

f
 L を代入
i
i
2
h
12
360
24
2
シンプソン公式の誤差
2 h 3  2 h 5 (4 ) 2 h 7 (6)
xi1 f (x)dx  2hfi  3 2 fi  5 4! fi  7 6! fi  L
h 3  fi 1  2 fi  fi 1 h 2 (4 )  2 h 5 (4 ) 2 h 7 (6)
 2hfi  

fi  
fi 
fi  L
2
3
h
12
7 6!
 5 4!
xi1
h
1  5 (4 )
 1
  fi 1  4 fi  fi 1   
  h fi  L
 60 36 
3
h
1 5 (4 )
  fi 1  4 fi  fi 1  
h fi  L
3
90
シンプソン公式
誤差
x0からxnまで誤差を足すと、全誤差Eは
n h 5 (4 ) n b  a h 4 (4 ) b  a 4 (4 ) b  a 4
E
f 
f 
h f 
h M
2 90
2 n 90
180
180
25
シンプソン公式の誤差
n h 5 (4 ) n b  a h 4 (4 ) b  a 4 (4 ) b  a 4
E
f 
f 
h f 
h M
2 90
2 n 90
180
180
Mは4階微分の最大値
シンプソン公式の誤差は刻み幅の4乗に比例
台形公式は2乗
26
シンプソン公式の例

 /2
0
x cos xdx
解析解は

 /2
0
x cos xdx  x sin x 0  
 /2
 /2
0
n
2
4
8
16
32
64
128
256
512
h
0.785398163397
0.392699081699
0.196349540849
0.098174770425
0.049087385212
0.024543692606
0.012271846303
0.006135923152
0.003067961576
sin xdx  x sin x 0   cos x 0 
 /2
積分値
0.581572016637
0.571416499264
0.570834320361
0.570798689642
0.570796474290
0.570796336010
0.570796327371
0.570796326831
0.570796326797
 /2

2
1
誤差
0.010775689842
0.000620172469
0.000037993566
0.000002362847
0.000000147495
0.000000009216
0.000000000576
0.000000000036
27
0.000000000002
シンプソン公式の例
4
0 1  x 2 dx
1
解析解はπ
n
h
2
4
8
16
32
64
128
0.500000000000
0.250000000000
0.125000000000
0.062500000000
0.031250000000
0.015625000000
0.007812500000
積分値
3.133333333333
3.141568627451
3.141592502459
3.141592651225
3.141592653553
3.141592653589
3.141592653590
誤差
-0.008259320256
-0.000024026139
-0.000000151131
-0.000000002365
-0.000000000037
-0.000000000001
0.000000000000
28
シンプソン公式の例
10-1
台形公式
シンプソン法
10-3
error
10-5
10-7
10-9
10-11
10-13 -3
10
10-2
10-1
h
100
29

Download Report