Transcript 数値計算法I - 九州工業大学

数値計算法I
第3回
ニュートン法
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
代数方程式
• f(x)=0 を満足するxを探す
f (x)  an x an1x
n
n1

a2 x a1x a0
2
n個の解をもつが、n≥5については、数値計算するしか
ない。

f (x)  e x sinx  2
このような場合は、数値計算に頼るしかない

•2分法
•ニュートン法
2
ニュートン法
• f(x)が連続であり、且つ微分が可能であれば、特
殊な場合を除いて必ず解が求まる。
•初期値x0を決める
•x0での接線を求める
•f(x0)から接線を引く
•接線がx軸と交差する点をx1とする
•x1での接線を求める
f(x0)
•f(xn)から接線を引く
•接線がx軸と交差する点をxn+1とする
•xn+1での接線を求める
x0 x1
x2
x3
3
ニュートン法
初期値をx0とする。
点(x0,f(x0))での接線の方程式は、
y  f (x0 )(x  x0 ) f (x0 )
x軸と交差する点x1は

0  f (x0 )(x1  x0 )  f (x0 )
f (x0 )(x1  x0 )   f (x0 )
f (x0 )
(x1  x0 )  
f (x0 )
f (x0 )
x1  x0 
f (x0 )
4
ニュートン法
i+1番目の値は、
f (xi )
xi1  xi 
f (xi )
値が真値に近づいていくと、f(xi)≈0なので、
xi+1はxiに比べて殆ど変化しなくなる。

|xi+1-xi|≤eまたは|f(xi+1)|<eになるまで計算を繰り返す
eは十分に小さな数（普通は自分で決める）
e:収束判定指数
5
ニュートン法
長所
• 2分法に比べて収束が早い。
• 重解でも可
短所
• まず微分しないといけない。
• 初期値によっては解が求まらない。
– f’(xi)=0となるxiに行ってしまうと発散する。
6
2分法
長所
• 愚直な方法ではあるが、殆どの場合に解が求まる。
– 電卓でも実行可能
短所
• 収束が遅い
• 偶数重解は求められない。
○
×
7
ニュートン法
f(x)=x3-6x2+7x+2=0
Y = x^ 3-6*x^ 2+7*x+2
4
case1
x0=0
2
y
0
case2
x0=2.5
-2
case3
x0=4.0
-4
-2
-1
0
1
2
3
4
5
6
x
真値は
2 5,2,2 5
8
ニュートン法
• f(x)=x3-6x2+7x+2=0
• f’(x)=3x2-12x+7
• case1において、
回数
1
2
3
4
5
f(x)
x
0.00000000 0.2000000E+01
-0.28571429 -0.5131195E+00
-0.23763999 -0.1573670E-01
-0.23606963 -0.1655032E-04
-0.23606798 -0.1837464E-10
9
ニュートン法
case2において、
回数
1
2
3
4
case3において
回数
1
2
3
4
5
x
2.50000000
1.94117647
2.00008159
2.00000000
x
4.00000000
4.28571429
4.23763999
4.23606963
4.23606798
f(x)
-0.2375000E+01
0.2939141E+00
-0.4079302E-03
0.1085354E-11
f(x)
-0.2000000E+01
0.5131195E+00
0.1573670E-01
0.1655032E-04
0.1838174E-10
10
ニュートン法
初期値を1とおくと、
回数
1
2
3
4
5
6
x
1.00000000
3.00000000
1.00000000
3.00000000
1.00000000
3.00000000
f(x)
0.4000000E+01
-0.4000000E+01
0.4000000E+01
-0.4000000E+01
0.4000000E+01
-0.4000000E+01
1と3の間を繰り返す。
11
ニュートン法（初期値への依存性）
Y = x^ 3-6*x^ 2+7*x+2
4
2
y
0
-2
-4
-2
-1
0
1
2
3
4
5
6
x
初期値が1または3の時は赤線のループを繰り返し
12
ニュートン法（初期値への依存性）
Y = x^ 3-6*x^ 2+7*x+2
4
2
y
0
-2
-4
-2
-1
0
1
2
3
4
5
6
x
5
2
初期値が
の時は、微係数がゼロなので、
3
x軸と交差しない。
13
ニュートン法の収束
f(x)=0の真値解をaとし、
f(xi)とf’(xi)をaのまわりでテーラー展開する。
(xi  a)2
f (xi )  f (a)  (xi  a) f (a) 
f (a)
2
0
(xi  a)2
f (xi )  f (a)  (xi  a) f (a) 
f (a)
2
i+1番目の誤差は

f (xi )
xi1  a  xi 
a
f (xi )
14
ニュートン法の収束
xi1  a  xi 
f (xi )
a
f (xi )
(xi  a)2
(xi  a) f (a) 
f (a) 
2
 xi  a 
f (a)  (xi  a) f (a) 
(x  a)
f (a)  i
f (a) 
2
 xi  a  (xi  a)
f (a)  (xi  a) f (a) 
キャンセル
(x  a) f (a)
1 i

2 f (a)
 xi  a  (xi  a)
(x  a) f (a)
1 i

f (a)
2次収束
(xi  a) f (a)
 (xi  a)
2 f (a)
f (a) 誤差xi+1-aは前のステップの
 (xi  a)2
2 f (a) 誤差xi-a の2乗に比例
15
ニュートン法の例
x2-cos(x)=0の解
回数
Y = x^ 2-cos(x)
2
1
2
3
4
5
6
1
y
0
-1
-2
-2
-1.5
-1
-0.5
0
x
0.5
1
1.5
x
f(x)
0.50000000 -0.6275826E+00
0.92420693 0.2516907E+00
0.82910576 0.1188097E-01
0.82414613 0.3292121E-04
0.82413231 0.2558270E-09
0.82413231 0.1110223E-15
2
16
ニュートン法のソース
c 計算を倍精度で実施
implicit real*8 (a-h,o-z)
c 初期値の入力
write(6,*)' input x0'
read(5,*)x0
c 繰り返し最大数
itermax=100
c 最小誤差
errlim=1.d-10
iter=0
xp=x0
c 繰り返しの開始
10 continue
iter=iter+1
c 新しいx値の計算
xn=xp-f(xp)/fd(xp)
c 出力
write(6,100)iter,xp,f(xp)
c 収束判定
if(dabs(f(xp)).lt.errlim)goto 30
if(iter.gt.itermax)goto 30
c x値の更新
xp=xn
goto 10
30 continue
c 回数、x値、f(x)の順で出力
c 整数は2桁を出力、実数は小数点以下8桁と小数点以下7桁の
c 回数で表示
100 format(i5,f19.8,e15.7)
stop
end
c f(x)を計算する関数サブルーチン
function f(x)
implicit real*8 (a-h,o-z)
f=x*x*x-6*x*x+7*x+2
return
end
c f'(x)を計算する関数サブルーチン
function fd(x)
implicit real*8 (a-h,o-z)
fd=3*x*x-12*x+7
return
end
17