Transcript x - 九州工業大学

数値計算法I
第8回
数値積分
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
数値積分
• 関数f(x)が複雑な形をしていて、解析的に積分し
にくい時、数値計算により(定)積分を行う。
I   f (x)dx
b
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 
2
i1

fn 

簡単化のため、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での微分は
fi1  fi h  h  h (4)
fi 
 fi  fi 
fi L
h
2
6
24
2
3
6
台形公式による数値積分の誤差
Ii  
xi1
xi
f (x)dx   f (xi  z)dz
h
0

z2  z 3  z 4 (4)
   fi  zfi fi  fi  fi L
0
2
3!
4!
h

dz
h2
h3  h4  h5 (4)
 hfi  fi fi 
fi 
fi L
2
6
24
120
fi1  fi h  h2  h3 (4)
fi 
 fi  fi 
fi L を代入すると、
h
2
6
24
fi1  fi h3  h4
Ii  h
 fi 
fiL
2
12
24
台形の面積
この分が誤差
7
台形公式による数値積分の誤差
fi1  fi h3  h4
Ii  h
 fi 
fiL
2
12
24
一区間辺りの誤差Ei
h3 
Ei  fi
12
全部で、
h3 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
hM
12
12
n 12
12
8
台形公式による数値積分の誤差
ba 2
E
hM
12
台形公式の誤差は刻み幅の2乗に比例する。
9
台形公式による数値積分（例）
4
0 1 x2 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
h
10-1
100
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)  y2a2 (x)  y3a3 (x)
ここで、a1(x),a2(x),a3(x)は2次式
P2(x)が、(x1,y1),(x2,y2),(x3,y3)を通る時、
P2 (x1 )  y1  y1a1(x1 )  y2a2 (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)  y2a2 (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

(xi1  xi1 )(xi1  xi ) xi1
xi1
19
シンプソン公式

xi1
xi1
xi1
fi1
f (x)dx  (x  x )(x  x ) xi1 (x  xi )(x  xi1 )dx
i1
i
i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  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
fi1 xi1
f (x)dx  2 xi1 (x  xi )(x  xi1 )dx
2h
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 3
2 
(x

x
)(x

x
)dx

(z

h)(z

2h)dz


hz

2h
z

h
i
i1
3 2

xi1
0

0 3
xi1
2h
2h
 z3
4 3
2
(x

x
)(x

x
)dx

(z)(z

2h)dz


hz


h
i1
i1
3

xi1
0
3

0
xi1
2h
2h
z h 2
2 3
(x

x
)(x

x
)dx

(z)(z

h)dz


z

h
i1
i
3 2 
xi1
0

0 3
xi1
2h
3
21
シンプソン公式

xi1
xi1
fi1 2 3 fi  4 2  fi1 2 2
f (x)dx  2 h  2   h   2 h
2h 3
h  3  2h 3
h
nは偶数
  fi1  4 fi  fi1 
3
よって、x0(=a)からxn(=b)までの積分値は、
0から2の積分

b
a
2から4の積分
4から6の積分
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 
3
i1
i1

fn 

22
シンプソン公式の誤差
fi+1とfi-1をxiの近傍でテイラー展開する
h2
fi1  fi  hfi
2
h2
fi1  fi  hfi
2
3
4
5
6
h
h
h
h
fi  fi  fi (4)  fi (5)  fi (6) L
3!
4!
5!
6!
3
4
5
6
h
h
h
h
(4)
(5)
(6)


fi  fi  fi  fi  fi L
3!
4!
5!
6!
和をとって、
4
6
h
h
2 
(4)
(6)
fi1  fi1  2 fi  h fi  fi 
fi L
12
360
差をとって、
h3  h5 (5)
fi1  fi1  2hfi fi 
fi
3
60
23
シンプソン公式の誤差

xi1
xi1
f (x)dx   f (xi  z)dz   f (xi  z)dz   f (xi  z)dz
h
h
0
h
0
h
h
z2  z 3  z 4 (4) z5 (5) z6 (6)
  ( fi  zfi fi  fi  fi  fi  fi L )dz
2
3!
4!
5!
6!
0
0
z2  z 3  z 4 (4) z5 (5) z6 (6)
  ( fi  zfi fi  fi  fi  fi  fi L )dz
2
3!
4!
5!
6!
h
2 h3  2 h5 (4) 2 h7 (6)
 2hfi 
fi 
fi 
fi L
32
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 h3  2 h5 (4) 2 h7 (6)
xi1 f (x)dx  2hfi  3 2 fi  5 4! fi  7 6! fi L
h3  fi1  2 fi  fi1 h2 (4)  2 h5 (4) 2 h7 (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 h5 (4) n b  a h4 (4) b  a 4 (4) b  a 4
E
f 
f 
h f 
hM
2 90
2 n 90
180
180
25
シンプソン公式の誤差
n h5 (4) n b  a h4 (4) b  a 4 (4) b  a 4
E
f 
f 
h f 
hM
2 90
2 n 90
180
180
Mは4階微分の最大値
シンプソン公式の誤差は刻み幅の4乗に比例
台形公式は2乗
26
シンプソン公式の例

 /2
0
xcos xdx
解析解は

 /2
0
xcos xdx  xsin 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  xsin 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
解析解はπ
n
h
2
4
8
16
32
64
128
シンプソン公式の例
1
4
0 1 x2 dx
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
h
10-1
100
29