コンピュータ初級 (クラス2)

IT入門B2
ー方程式の反復解法ー
授業内容
•反復解法
•ニュートン法
•プログラムの作成
•演習
非線形方程式
連立一次方程式を解く場合には，ガウスの消去法等といった直接解法によ
り解の計算ができる．
非線形方程式を解く場合は直接解法が存在しないため，反復解法を用いて
解を数値的に計算する．
今回の授業では，反復解法の中でもよく知られたニュートン法を用いて方
程式の解を求めるプログラムを作成する．
反復解法
方程式
f ( x)  0,
f :R R
が与えられたとする．
(i )
適当な初期値 x ( 0 )から出発し，f ( x)  0のある解 x* に収束する列 x をつく
り，x (i )が x* に十分近づいたときに計算を打ち切り x (i ) を近似解とする．
 
収束の判定は，ある微小な値   0 に対し，例えば次のような条件
とする：
(1) x (i 1)  x (i )  
(2) x (i 1)  x (i )   x (i )
(3) f ( x (i ) )  
このとき，上の条件がすべて満たされても，
x (i )  x*  
が成り立つとは限らない．
ニュートン法
― 非線形方程式を解くための反復解法の一つ．
f : R  Rを1回連続微分可能な関数とし，方程式
f ( x)  0
を解くことを考える．
適当な初期点 x ( 0)  R を選び，これを方程式の最初の近似解とする．
ただし， f ( x(0) )  0 とする．
このとき f ( x(0) )  0 ならば，x ( 0 )を修正し
x (1)  x (0)  h
とし，方程式 f ( x)  0のより良い近似解となるようにする．
ニュートン法
f を x ( 0 ) のまわりでTaylor展開すると，
2
f ( x ( 0) )
f ( x ( 0) )
( 0)
( 0)


f ( x)  f ( x ) 
x  x 
x  x (0)   
1!
2!
となる．ここで，x  x (1)として1次の項で打ち切ると，
f ( x (1) )  f ( x ( 0 )  h)
 f ( x ( 0 ) )  f ( x ( 0 ) ) h
と近似できる．このとき，(右辺  0) になってほしいわけだから，
f ( x ( 0) )  f ( x ( 0) )h  0
f ( x (0) )
h
f ( x ( 0) )
と定まる．よって，
x
(1)
である．
x
( 0)
f ( x ( 0) )

f ( x ( 0) )
ニュートン法
ここで，
g ( x)  x 
f ( x)
f ( x)
とおくと，上式は
x(1)  g ( x(0) )
と書ける．
これを繰り返して一般化すると
x
( i 1)
f ( x (i ) )
 g(x )  x 
f ( x (i ) )
(i )
(i )
i  0,1,2,3,
となる．この反復解法をニュートン法と呼ぶ．
ニュートン法
x
( i 1)
f ( x (i ) )
 g(x )  x 
f ( x (i ) )
(i )
(i )
i  0,1,2,3,
 
ニュートン法により得られる点列 x (i ) が点 x*に収束したとする．
このとき
f ( x* )
x  g(x )  x 
f ( x* )
*
*
*
となる．すなわち
f ( x* )  0
したがって，x*が f ( x)  0 の解となっていることがわかる．
ニュートン法
f (x)
- 幾何学的な解釈 -
x (i 1)
x*


曲線 y  f (x) 上の点 x(i ) , f ( x(i ) ) における接線
y  f ( x(i ) )(x  x(i ) )  f ( x(i ) )
とx軸との交点のx座標
f ( x (i ) )
xx 
f ( x (i ) )
(i )
を求める操作の反復を行っている．
x (i )
ニュートン法
- 参考 r  R を定数として
g ( x)  x  rf ( x)
とする．
反復解法として
i  0,1,2,3,
x(i 1)  g ( x(i ) )  x(i )  rf ( x(i ) )
を用いる方法を簡易ニュートン法と呼ぶ．
r
1
f ( x ( 0) )
とする場合や，r を
 f ( x
( 0)
)0

1
のなんらかの近似とする場合が多い．
(0)
f ( x )
プログラムの作成
以下の例題を解くためのプログラムを作成してみよう．
例題
方程式
x  1.2e  x
の解をニュートン法により求める．
x  1 .2 e  x
 x  1 .2 e  x  0
より， f ( x)  x 1.2e x とすると
f ( x)  1  1.2e x
であり，
x
( i 1)
f ( x (i ) )
 g(x )  x 
f ( x (i ) )
(i )
(i )
で定まる反復により解を求める．
i  0,1,2,3,
プログラムの作成
・ f ( x), f ( x) の計算：
f ( x), f ( x) はそれぞれ
f ( x)  x  1.2e  x
f ( x)  1  1.2e  x
で与えられている．これらの値の計算には，f ( x), f ( x)をそれぞれ
関数として定義することにしよう．
f ( x)  x 1.2e xを以下のような関数として定義できる：
double func(double x)
{
double y;
y = x – 1.2*exp(-x);
return y;
}
f (x )も同様に関数名func_d()等として定義してみよう．
プログラムの作成
このとき，関数func()をmain関数等から呼び出すことでf xの値
を計算できる．
double func(double x)
{
double y;
y = x – 1.2*exp(-x);
return y;
}
int main(void)
{
double f, x;
x = 1.0;
f = func(x);
…
← fにf x)の値が代入される
プログラムの作成
f ( x (i ) )
i  0,1,2,3,
x
x 
(i )
f ( x )
ニュートン法では，適当な初期点 x ( 0 )を選び，上式にしたがってi の値を順
次変化させながら次の点を次々と計算する．
適当な停止条件を満たしたところで計算を止め，そのときの値を解とする．
( i 1)
(i )
・初期点の設定：
現在いる点を表す変数xをdouble型で宣言し，xに適当な値を代入．
int main(void)
{
double x;
…
x = 1.0;
…
・ニュートン法の反復：
while文，for文等の繰り返し処理を利用して反復を実現する．
プログラムの作成
・反復の停止条件：
ニュートン法の反復を停止する条件として
x ( i 1)  x ( i )  
を利用することにする．
現在いる点xから，ニュートン反復により計算される次の点を表す
変数をx_nとして宣言している場合，
fabs(x_n – x) < e
が成立したときに反復を停止すれば良い．
このとき，eとして予め微小な値を定義しておく必要がある．
← 106の指数表現
e = 1e-6;
while (…) {
/* ニュートン反復により次の点x_nを計算する部分 */
if (fabs(x_n - x) < e)
…
}
break;
プログラムの作成
・最大反復回数の設定：
ニュートン法では，初期点の選び方により解に収束しない場合がある．
その場合を考慮して，最大の反復回数を設定しておき，解に収束しな
かった場合は強制的に終了させる．
例えば，変数maxに最大の反復回数を代入しておき，繰り返し文の
終了条件として利用する．
max = 100;
i = 0;
while (i<max) {
…
i++;
}
max = 100;
for (i=0; i<max; i++) {
…
}
繰り返し文が終了したとき，iの値がmaxに等しい場合は解に収束して
おらず，それ以外の時は解に収束している．
プログラムの作成
例題を解くためのプログラムの流れ
int main(void)
{
double e, x, x_n;
int i, max;
x = 1.0;
e = 1e-6;
max = 100;
if (i < max) {
/* 解に収束 */
} else {
/* 解に収束していない
}
i = 0;
while (i < max) {
/* ニュートン反復により
次の点x_nを計算する部分 */
if (fabs(x_n - x) < e) break;
/*
}
変数の更新
*/
return 0;
}
*/
プログラムの作成
プログラム作成上の注意
・exp(), fabs()等の数学関数を利用するにはプログラムの最初に
#include <math.h>
と記述．
・また，コンパイル時に-lmオプションを付け、
gcc -lm …
とする．
・詳細は，第3回の講義資料を参照のこと．
演習
1
例題を解くためのプログラムを完成させて，例題の解を計算してみよう．
2
方程式
x  cos(x)  0
の解をニュートン法により計算してみよう．

Download Report