Newton-Raphson法

プログラミング論 I
2008年6月12日
非線形方程式の近似解
Newton-Raphson法
http://www.ns.kogakuin.ac.jp/~ct13140/Prog.2008
for文復習
• 回数が決まっている繰り返し処理(n回とす
る)は, 必ず for(i=0; i<n; i++) とす
るのが分かりやすい(ことが多い).
2
for文復習
• 100, 102, 104, 106...198 の50
個を表示
–例 A
for(i=100; i<200; i+=2){
printf("%d\n", i);
}
3
for文復習
• 100, 102, 104, 106...198 の50
個を表示
–例 B
for(i=0; i<50; i++){
x = 100 + i*2;
printf("%d\n", x);
}
4
練習 0
右の様に表示されるプログラムを
記述せよ.
倍精度浮動小数点の表示は,
double x;
printf("%lf\n", x);
6.000000
6.300000
6.600000
6.900000
7.200000
7.500000
7.800000
8.100000
8.400000
8.700000
5
解答 0
int i;
double x;
for(i=0; i<10; i++){
x = 6.0 + 0.3*i;
printf("%lf\n", x);
}
6
概要
• 非線形方程式の近似解
– Newton-Raphson法
• 補間
– n次多項式での近似
テイラー展開
– 関数の補間
ラグランジュの補間法
7
Newton-Raphson法
Newton-Raphson Method
8
Newton-Raphson法
•
•
•
•
反復法の一種.
必ず収束するとは限らない.
初期値が解に近ければ高速に収束する.
n回目の調査で使った点から1次近似直線
を引き,そのx切片をn+1回目の候補とする
• 初期値は1個の値.
– 2分法やはさみうち法の様に区間ではない.
9
Newton-Raphson法
方程式を f (x) = 0 とする
(1) 解に近い値 a0 を定める.
(2) y = f (x) に対して,( an , f (an) ) で接する
接線(1次近似直線)を引く.
接線は y – f (an) = f '(an) (x–an)
10
Newton-Raphson法
(3) 近似直線がx軸と交わる点(an+1,0)を求める.
すなわち,折線 y – f (an) = f ’(an) (x–an) が
y=0となる x を求め,その x がan+1である.
f ( an )
an 1  an 
f ' ( an )
近似直線と y = f (x) が近ければ,
an+1は解に近いことが期待される.
11
Newton-Raphson法
(4) |an-an+1| が十分に小さければ,終了.
小さく無ければ,(2)に戻る
12
Newton-Raphson法の例
• f (x) = 0.1x3+x-5 とし f (x) = 0 の解を探す
60
50
40
f(x)
30
20
10
0
0
1
2
3
4
-10
x
5
6
7
8
13
Newton-Raphson法の例
• f '(x) = 0.3x2+1
• 初期値a0=7とする.
• y = f (x) と ( 7, f (7) )で接する接線を引き,
接線がx軸と交わる点(a1, 0)を求め,
このa1に着目する.
14
Newton-Raphson法の例
ただし,( 7, f (7) )における接線は
y – f (7) = f ’(7)(x – 7)
接線のx切片は
f (7 )
a1  7 
≒ 4.69
f ' (7 )
15
近似直線 と y= f (x) が
近ければ,
近似直線のx切片は
解に近いと期待される.
60
50
40
( 7, f (7) )
f(x)
30
20
10
( 7, f (7) )で
接する接線
0
0
1
2
3
4
5
6
7
8
-10
x
次に着目する点
16
Newton-Raphson法の例
• y = f (x) と (a1, f (a1) )で接する接線を引き,
接線がx軸と交わる点(a2, 0)を求め,
次は,このa2に着目する.
17
60
50
40
f(x)
30
20
10
0
0
1
2
3
4
5
6
7
8
-10
x
18
6
5
4
f(x)
3
2
1
0
2
2.5
3
3.5
4
-1
-2
x
19
Neton-Raphson法で収束しない例
接点
接線
次の候補
次の候補
接線
接点
20
Newton-Raphson法の特徴
• 接線(1次近似直線)が,方程式と近ければ,高速
に収束していく.
• 初期値が解に近くないと,収束しない.
• 微分できないとNewton-Raphson法は使えない.
– 傾きがゼロの箇所においても使えない.
• 「 | f (an)|が十分に小さいこと」を収束条件するこ
ともある.
• 厳密には「|an-an+1|や f (m)が十分に小さいこと」
は,解の精度を保証していない.
21
復習:2分法
• 2分法
f (x)=ーx3ーx2+5x+20
として f (x)=0 の解を探す
(注意:f (x) は単調増加ではない)
ただし解が[-4,4]に存在は既知
50
40
30
f(x)
20
10
0
-4
-3
-2
-1 -10 0
1
2
3
4
-20
-30
x
22
復習:2分法
• 解は[-4,4]で,f (-4)>0, f (4)<0
中間値0に着目.
f (-4)>0, f (0)>0, f (4)<0 なので,解は[0,4]に
• 解は[0,4]で,f (0)>0, f (4)<0
中間値2に着目.
f (0)>0, f (2)>0, f (4)<0 なので,解は[2,4]に
• 解は[2,4]で,f (2)>0, f (4)<0
中間値3に着目.
f (2)>0, f (3)<0, f (4)<0 なので,解は[2,3]に
23
復習:2分法
min=??, max=??;
(max-min)が大きいなら繰り返す
mid=(min+max)/2;
f(min)*f(mid)>0 同符号なら
解は[mid,max]
min = mid;
f(min)*f(mid)<0 異符号なら
解は[min,mid]
max = mid;
繰り返しの先頭に戻る
24
補間
25
補間
• (1) n次多項式での近似
– テイラー展開
• (2) 関数の補間
– ラグランジュの補間法
26
補間とは
• (x0, f (x0)), (x1, f (x1),…(xn, f (xn))が既知の
とき,
それ(x0, x1,… xn)以外の点における関数の値
を求める(予測する)
27
補間の例
x=6
10
f(x)
8
(10, 6.6)
6
4
f (5), f (10)が既知.
f (6)が必要になったら?
近くの f (5)=2.4で代用する?
(5, 2.4)
2
0
0
5
10
15
x
20
25
30
28
f (5)と f (10)からf (6)を補間する
おそらく f (5)より,
一次近似の方が正確.
さらに正確な補間方法は?
7
6
(10,6.6)
f(x)
5
直線を引いて(一次近似),
x=6 における値を予想.
4
3
2
(5,2.4)
f (5)=2.4 で代用
1
4
6
8
x
10
12
29
(1) n次多項式での近似
・近似は補間と密接な関係.
・関数をn次多項式で近似する.
関数の式が既知であることが前提
・テイラー展開
30
一次式で近似
接線
1.5
y =f (a)+f '(a)(x-a)=g(x)
1
0.5
y
近似に
用いた点
0
0
-0.5
-1
y = f (x)
x

2

y = f (x) の接線 y = g(x)は,y =f (x)の一次近似線.
f (a) = g(a) かつ f '(a) = g'(a) となっている.
つまり,x=a において,
「yの値」と「傾き(一階微分)」が一致している.
31
二次式で近似
y  f (a)  f ' (a)( x  a)
1.5
y =f (a)+f '(a)(x-a)=g(x)
1
0.5
y
近似に
用いた点
0
0
-0.5
-1
x
y = f (x)

2
1
 f ' ' (a)( x  a) 2
2
 h( x )
x=aに近いほど
近似は正確

y = h(x)は,y =f (x)の二次近似曲線.
(y=f (x)を近似した二次曲線)
f (a) = h(a) かつ f '(a) = h'(a) かつ f ''(a) = h''(a).
x=a において 「y」と「一階微分」と「二階微分」が
が一致している.一次近似よりより正確な近似. 32
三次式で近似
1
1
2
i ( x)  f (a )  f ' (a )( x  a )  f ' ' (a )( x  a )  f ' ' ' (a )( x  a ) 3
2
6
1
i ' ( x)  f ' (a )  f ' ' (a )( x  a )  f ' ' ' (a )( x  a ) 2
ここは
2
xの式
i ' ' ( x)  f ' ' (a )  f ' ' ' (a )( x  a )
ここは
i ' ' ' ( x)  f ' ' ' (a)
定数である
y = i(x)は,y =f (x)の三次近似曲線.
f (a)=i (a),f '(a)= i (a), f ''(a)= i ''(a),f '''(a)=i '''(a)
よって,x=a において,yの値,1~3階微分が一致.
33
テイラー展開 Taylor expansion
• f (x) を以下の級数に展開する


k 0
f
(k )
(a)
k
( x  a)
k!
a=0の場合,
マクローリン展開
1
1
2
3
f (a )  f ' (a )( x  a )  f ' ' (a)( x  a )  f ' ' ' (a )( x  a )
2
6
1
1 (5)
4
5

f ' ' ' ' (a )( x  a ) 
f (a )( x  a)  ...
34
24
120
テイラー展開
• テイラー展開を有限個の級数で打ち切る
 n1 f ( k ) (a)

f ( x)   
( x  a) k   Rn
 k 0 k !



Rnはラグランジュの剰余項
• lim Rn =0のとき, f ( x) 
n 
 f ( k ) (a)

k 0
k!
( x  a) k
lim R
換言すると nn =0のとき,展開の次数を上げていくと
より正確な近似になっていく.
35
テイラー展開 近似線
1次
2
2次
1.5
1
0次
0.5
5次
6次
0
-0.5
-1
-1.5
-2
0

0.5
2
1
y=sin(x)
3次
4次
36
n次多項式による近似
• 基本的に近似点の近傍では次数を上げる
ほど,正確な近似になる
37
単語の説明
10
f(x)
8
6
4
この区間を
補うことを
補外,外挿
という
この区間を補うことを
補間,内挿という
2
0
0
5
10
15
x
20
25
30
38
(2) 関数の補間
関数の式が未知の場合も使える
39
補間
既知の f (xi)から,f (6)を補間する例を考える.
10
f(x)
8
6
4
2
0
0
5
10
15
x
20
25
30
40
0次補間
x=6に一番近い,x=5を採用する.
7
6
既知の点
f(x)
5
既知の点
4
3
f (6)を補間
2
f (6) = f (5)
1
4
6
8
x
10
12
41
1次補間
(5, f (5)), (10, f (10)) を通る1次関数(直線)で補間する.
7
近似直線を用いて
f (6)を補間
6
既知の点
f(x)
5
既知の点
4
3
f (5)  f (10)
f (6) 
(6  5)  f (a)
5  10
2
1
4
6
8
x
10
12
42
1次補間
• 2点(a , f (a))と(b , f (b))より,x=mにおける
yの値 f (m)を予測する
• 2点(a , f (a))と(b , f (b)) を通る直線の式
f (a)  f (b)
y  f (a) 
( x  a)
a b
これに,x=mを代入して,
f (a)  f (b)
f ( m) 
(m  a)  f (a)
a b
43
2次補間
3点を通る2次関数を求め,その2次関数(放物線)で補間する.
8
f(x)
6
(5, f (5)), (10, f (10)),
(15, f (15)) を通る2次曲線.
これを近似曲線とする
4
近似曲線を用いて
f (6)を補間
2
0
4
6
8
10
12
x
14
16
44
2次補間
•3点(x0, y0), (x1, y1), (x2, y2) を通る2次曲線の
求め方(の例)
2次曲線は y=ax2+bx+c で表せる.
(3個の未知数が決まれば,2次曲線が決まる)
(x, y)=(x0, y0), (x, y)=(x1, y1), (x, y)=(x2, y2)を代入す
れば,方程式が3本立つ.
よって,3元1次連立方程式を解けば2次曲線は求
まる.
後で,別の方法(ラグランジュの補間法)を紹介する 45
2次補間
• 求めた2次近似曲線 y=g (x) に,
x=m を代入し,g(m)を補間値として用いる.
• 1次補間(直線による予測)より,正確である
ことが期待できる.
46
練習 1
y=f (x)は,2点 (1, 1), (2,4) を通る.
f (1.5) と f (3) を,一次補間(補外)
により予測せよ 8
7
6
f(x)
5
4
3
2
1
0
0
1
2
x
3
47
解答 1
• 両点を通る直線
傾き
4 1
3
2 1
(1, 1)を通り,傾き3の式は
y  1  3( x  1)
y  3x  2
x = 1.5 のとき,y = 2.5
x = 3 のとき,y = 7
48
ラグランジュの補間法
既知の点n+1個を全て通る
n次式を求める
49
ラグランジュの補間法
• n+1個の点を通る,n次方程式P(x)で補間
通る点を(x0, f (x0)), (x1, f (x1)),…, (xn, f (xn)) として
以下の方法により,P(x)が求まる
P( x) 
n
 f ( xk ) Lk ( x)
k 0
x  xm
ただし Lk ( x)  m
0 xk  xm
n
m k
50
ラグランジュの補間法
x  xm
Lk ( x)  
m  0 xk  xm
n
m k
( x  x0 )  ( x  xk 1 )( x  xk 1 )  ( x  xn )

( xk  x0 )  ( xk  xk 1 )( xk  xk 1 )  ( xk  xn )
(xーxk)以外
(xkーxk)以外
Lk(x)の分子は x の n次多項式.分母は定数.
よって,Lk(x)はxのn次多項式.
f (xk) は定数.よって, f (xk)Lk(x)は x の n次多項式.
よって,P( x) 
n
 f ( xk ) Lk ( x) は x の n次多項式.
k 0
51
ラグランジュの補間法の例
• 4点(x0, f (x0)), (x1, f (x1)), (x2, f (x2)), (x3, f (x3))を
通る,3次曲線の方程式 y = P(x) は,
( x  x1 )( x  x2 )( x  x3 )
y  P( x)  f ( x0 )
( x0  x1 )( x0  x2 )( x0  x3 )
( x  x0 )( x  x2 )( x  x3 )
 f ( x1 )
( x1  x0 )( x1  x2 )( x1  x3 )
( x  x0 )( x  x1 )( x  x3 )
 f ( x2 )
( x2  x0 )( x2  x1 )( x2  x3 )
(xーx2)
以外
(x2ーx2)
以外
( x  x0 )( x  x1 )( x  x2 )
 f ( x3 )
( x3  x0 )( x3  x1 )( x3  x2 ) 52
ラグランジュの補間法の例
• y = P(x) は,(x2, f (x2))を通るか確認
( x  x1 )( x  x2 )( x  x3 )
y  P( x2 )  f ( x0 )
( x0  x1 )( x0  x2 )( x0  x3 )
y=P(x)は
(x2 , f (x2))
を通る
3次曲線
( x  x0 )( x  x2 )( x  x3 )
 f ( x1 )
( x1  x0 )( x1  x2 )( x1  x3 )
( x  x0 )( x  x1 )( x  x3 )
 f ( x2 )
( x2  x0 )( x2  x1 )( x2  x3 )
( x  x0 )( x  x1 )( x  x2 )
 f ( x3 )
( x3  x0 )( x3  x1 )( x3  x2 )
ゼロ
ゼロ
1
(x=x2なら
分子と分母
が一致)
ゼロ
= f (x0)×0 + f (x1)×0 + f (x2)×1 + f (x3)×0 = f (x2)53
ラグランジュの補間法の例
10
f(x)
8
4次曲線で
近似し,補間
6
4
2
0
0
5
10
15
x
20
25
30
54
ラグランジュ補間法の問題点
1.5
補間曲線は,
確かに
全ての点を
通過している
y
1
端の方は
全く近く
いない
1
2
x 1
を11点を用いて
ラグランジュ補間法で補間
y=1/(x2+1)
0.5
Lagrange補間
補間点
0
-5
-3
-1
-0.5
1
3
正解
5
ラグランジュ補間法
による予想55
ラグランジュの補間法の問題点
• 補間の点数が増えると,激しく振動する
– 点数が増えると,
(全ての補間点を通過しているのは事実だが)
補間のn多項式の次数が大きくなり振動する.
– 次数が高いほど好ましいのではない.
56
(3次)スプライン補間
1.5
この区間を
3次関数
で補間
区間ごとに,別々の
3次多項式を用意し
補間する
1
y=1/(x2+1)
この区間を
3次関数
で補間
補間点
0.5
0
-5
-3
-1
1
-0.5
3
5
57
練習 2
• 3点 (x0, y0), (x1, y1), (x2, y2) を通る,2次方
程式を書け.
– ただし,  記号と, 記号は使用しないこと
58
解答 2
( x  x1 )( x  x2 )
f ( x0 )
( x0  x1 )( x0  x2 )
( x  x0 )( x  x2 )
 f ( x1 )
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )
 f ( x2 )
( x2  x0 )( x2  x1 )
59