非線形方程式の数値計算（ニュートン法）（pdf

非線形方程式の数値計算
（ニュートン法）
電気情報工学科 ∗ コンピュータシミュレーション
2015 年 7 月 10 日 (金)
概要
非線形方程式の近似解を，数値計算により求める．特に，二分法とニュートン法
により近似解を求めるプログラムを作成する．
1 講義内容概略
数値計算による非線形方程式の近似解を求めてみる．今回はニュートン法によるプログ
ラムの作成を行う．特に以下の点について修得できることを目標とする．
• 二分法とニュートン法の計算原理をしっかりと理解し，それらを説明することがで
きる．
• アルゴリズムのフローチャートが書ける．
• C 言語でプログラムを作成し，実際に近似解を得ることができる．
2 ニュートン法
2.1
2.1.1
計算方法
漸化式
関数 f (x) のゼロ点 α に近い近似値 x0 から出発する．そして，関数 f (x) 上の点 (x0 , f (x0 ))
での接線が，x 軸と交わる点を次の近似解 x1 とする．そして，次の接線が x 軸と交わる点
を次の近似解 x2 とする．同じことを繰り返して x3 , x4 , · · · を求める (図 1)．この計算結果
∗
秋田工業高等専門学校
1
の数列 (x0 , x2 , x3 , x4 , · · · ) は初期値 x0 が適当であれば，真の解 α に収束することは図より
明らかであろう．
実際にこの数列を計算するためには，漸化式が必要である．図のようにすると，関数
f (x) 上の点 (xi , f (xi )) の接線を引き，それと x 軸と交点が xi+1 になる．まずは， xi+1 を求
めることにする．点 (xi , f (xi )) を通り，傾きが f 0 (xi ) の直線の方程式は，
y − f (xi ) = f 0 (xi )(x − xi )
(1)
である．図 1 より，y = 0 の時の x の値が xi+1 になることが分かる．したがって， xi+1 は，
xi+1 = xi −
f (xi )
f 0 (xi )
(2)
となる．これで，xi から xi+1 が計算できる．これをニュートン法の漸化式と言う．この漸
化式を用いれば，次々とより精度の良い近似解を求めることができる．
計算の終了は，
xi+1 − 1 ≤ ε
(3)
xi
の条件を満たした場合とするのが一般的である ε は計算精度を決める定数で，非常に小さ
な値をプログラマが設定する．これ以外にも計算の終了を決めることは可能なので，状況
に応じて，計算の打ち切り方法 (例えば break) を使えばいい．ニュートン法による計算収
束の様子を図 1 に示す．
150
125
100
75
50
25
-1
1
2
x3
4
3
x2
x1
5
x0
6
-25
図 1: f (x) = x3 − 3x2 + 9x − 8 の実数解をニュートン法で計算し，解の収束の様子を示して
いる．初期値 x0 = 5 から始まり，接線と x 軸の交点からより精度の高い回を求めている．
2
2.1.2
解への収束速度
ニュートン法を使う上で必要な式は，式 (2) のみである．計算に必要な式は分かったが，
数列がどのように真の解 α に収束するか考える．xi+1 と真値 α の差の絶対値，ようするに
誤差を計算する． f (α) = 0 を忘れないで，右辺 f (xi )/ f 0 (xi ) を α の周りでテイラー展開す
る．すると，
f (xi ) |α − xi+1 | = α − xi + 0
f (xi ) [
]
(
)
f (α)
f (α) f 00 (α)
(4)
= α − xi + 0
+ 1−
(xi − α) + O (α − xi )2 f (α)
f 02 (α)
(
)
= O (α − xi )2 となる．i + 1 番目の近似値は，i 番目に比べて 2 乗で精度が良くなるのである．これを，
二次収束と言い，非常に早く解に収束する．例えば，10−3 の精度で近似解が得られてい
るとすると，もう一回式 (2) を計算するだけで，10−6 程度の精度で近似解が得られるとい
うことである．一次収束の二分法よりも，収束が早いことが分かる．
2.1.3
解けない方程式
アルゴリズムから，二分法は解に必ず収束する．ただし，収束のスピードが遅いのが欠
点である．一方，ニュートン法は解に収束するとは限らない．初期値が悪いと解に収束し
ない場合がある．厳密にその条件を求めるのは大変なので，実例を示すことにする．
非線形方程式
3 tan−1 (x − 1) +
x
=0
4
(5)
を計算することを考える．これは，初期値が悪いと収束しない方程式の例である．例えば
初期値 x0 = 3 の場合，図 2 のように収束しない 1 ．これを初期値 x0 = 2.5 にすると図 3 の
ように収束する．
このようにニュートン法は解に収束しないで，振動する場合がある．こうなると，プロ
グラムは無限ループに入り，永遠に計算し続ける．これは資源の無駄遣いなので，慎むべ
きである．通常は，反復回数の上限を決めて，それを防ぐ．ニュートン法を使う場合は，
この反復回数の上限は必須である．
ニュートン法で収束する必要条件が分かればこの問題は解決する．しかし，それを探
すのは大変である．というか私には分からない．一方，十分条件は簡単にわかる．閉区間
[a, b] で， f (a) < 0, f (b) > 0 のような関数を考える．このとき，
• 常に f 0 (x) > 0, f 00 (x) > 0 で，初期値が x0 = b の場合
• 常に f 0 (x) > 0, f 00 (x) < 0 で，初期値が x0 = b の場合
発散だと思う．もしかしたら振動かも．私は興味がありませんので，発散か振動か誰か証明してくださ
い．
1
3
は，必ず収束する．ニュートン法の図から明らかである． f (a) > 0, f (b) < 0 の場合は，
f (x) に-1 倍すれば，先の十分条件を考えることができる．
実際には収束しない場合のほうが稀であるので，ニュートン法は非常に強力な非線形方
程式の解法である．ただ，反復回数を忘れないことが重要である．また，二分法と組み合
わせて使うことも考えられる．
4
5
x5
-15
x3
-10
-5
x1
2
x2
x0 5
x4
10
x1
x6
15
-3
-2
x5
x3
-1
1 x4
x2
2
x0
3
-2
-5
-4
図 2: ニュートン法で解が求まらない場合
2.2
図 3: ニュートン法で解が求まる場合
フローチャート
二分法同様，関数と計算を打ち切る条件はプログラム中に書くものとする．そうする
と，図 4 のようなニュートン法のフローチャートが考えられる．なお，漸化式を取扱う場
合，一次元配列を利用すると便利である．
始め
i = -1
x0 を入力
i ← i+1
xi+1 = xi - f(xi) / f’(xi)
解の精度検査
(xi+1- xi) / xi ≦
解 xi+1 を表示
no
yes
yes
反復回数検査
i ≦ imax
no
収束しないと
表示
終り
図 4: ニュートン法のフローチャート
4
2.3
プログラムの穴埋め問題
リスト 1: ニュートン法のプログラム（穴埋め）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# i n c l u d e < s t d i o . h>
# i n c l u d e <math . h>
# d e f i n e EPS 1 e −7
# d e f i n e IMAX 100
/ / 計算精度
/ / 計算回数の上限
/ / ユーザー定義関数のプロトタイプ宣言
void b i s e c t i o n ( ) ;
v o i d newton ( ) ;
double f ( double x ) ;
double fd ( double x ) ;
/ / ===================================
/ / メイン関数
/ / ===================================
i n t main ( ) {
newton ( ) ;
return 0;
}
/ / =====================================
/ / 解くべき非線形関数
/ / =====================================
double f ( double x ) {
double y ;
y = x + exp ( x ) + s i n ( x ) ;
return y ;
}
/ / ====================================
/ / 解くべき非線形関数の導関数
/ / ====================================
double fd ( double x ) {
double y ;
y = 1 + exp ( x ) + c o s ( x ) ;
return y ;
}
/ / ======================================
/ / ニュートン法による計算
/ / ======================================
v o i d newton ( ) {
d o u b l e x [IMAX ] ;
i n t i = −1;
FILE ∗ f p ;
5
51
52
53
54
55
56
57
58
59
60
61
62
63
64
f p = f o p e n ( ” newton . d a t ” , ”w” ) ;
p r i n t f ( ” ニュートン法による計算を行います． \ n ” ) ;
p r i n t f ( ” の初期値を入力してください： x ” ) ;
s c a n f ( ”%l f %∗c ” ,&x [ 0 ] ) ;
/ / ターミナルへの表示
/ / ターミナルへの表示
/ / 初期値の入力
/ / −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
/ / ここにニュートン法の計算とファイル出力のプログラムを書く
/ / −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
}
f c l o s e ( fp ) ;
p r i n t f ( ” ニュートン法による解：%0.15 f \ n ” , x [ i + 1 ] ) ; / / 繰返しごとの計算結果のターミナル表示
p r i n t f ( ” ニュートン法による計算回数：%d \ n ” , i ) ;
/ / 繰返しごとの計算回数のターミナル表示
3 練習問題
[練習 1] 下記の非線形方程式をニュートン法と二分法で解きなさい．
[練習 2] ニュートン法と二分法について，横軸に計算回数，縦軸に計算結果の解
を同時にプロットしたグラフを作成し，解の収束速度について考察しな
さい．また，計算結果をファイル出力し，グラフは gnuplot を用いること．
今回解く方程式:
f (x) = x + exp(x) + sin(x)
(6)
4 レポート課題
上記の練習問題をプログラムリストと結果を示したレポートを作成し，提出する．特
に，体裁がしっかりした報告書を提出すること．
提出期限：2015 年 7 月 17 日の講義終了時まで．
提出場所：坂本研究室入口の段ボールのレポート入れ．または講義中に手渡し．
6

Download Report