方程式の反復解法（PDF）

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)
f ( x) = f ( x ) +
x−x +
x − x (0) + L
1!
2!
(0)
(
)
(
)
となる．ここで，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,L)
となる．この反復解法をニュートン法と呼ぶ．
ニュートン法
x
( i +1)
f ( x (i ) )
= g(x ) = x −
f ′( x ( i ) )
(i )
(i )
(i = 0,1,2,3,L)
{ }
ニュートン法により得られる点列 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,L)
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. 2e − x
⇔ x − 1. 2e − 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,L)
プログラムの作成
・ f ( x), f ′( x) の計算：
f ( x), f ′( x) はそれぞれ
f ( x ) = x − 1. 2 e − 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,L)
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
の解をニュートン法により計算してみよう．