プログラム演習 (4) 三次方程式の解を二分法とニュートン法を用いて求めるプログラム演習 (1) では二次方程式の解を簡単なアルゴリズムで求めたが、今回はもう少しだけ賢いアルゴリズムで求めることとする。今回採り上げるのは、二分法とニュートン法である。なお、本演習では、解くべき問題を「f (x) = 0を満足する x の値を求める（ただし f (x) は連続関数）」というものとする。
はじめに二分法について説明する。
ある２つの x について、その関数値 f (x) の符号が互いに異なるようなものを探し、これを x1 と x2 とする。このようなときには、f (x) = 0 の解のひとつは必ず区間 [x1 , x2 ] に存在するということが言える。そこで、x1 と x2 の中点をx3 とし、そこでの関数値 f (x3 ) を計算する。下の図のように、f (x3 ) と f (x2 ) の符号が同一である場合には、f (x) = 0の解のひとつは [x1 , x3 ] に存在することが言え、解の存在（推定）範囲を 1/2 に狭めることが出来る。逆に f (x3 ) とf (x2 ) の符号が異なる（f (x3 ) と f (x1 ) の符号が同一である）場合には、f (x) = 0 の解のひとつは [x3 , x2 ] に存在することが言える。この操作を繰り返すことによって解の存在範囲を小さくすること、つまり精度の高い解の推定が行えることになる。これが二分法の考え方である。

Fig. 1: 二分法の考え方


問題１：二分法に基づいて x3 − 4x2 − x + 4 = 0 の解を求めるプログラムを作成し、解を数値的に求め、解析解と比較せよ。


ゼロから自分でプログラムを作成するのが難しい場合は、以下のサンプルプログラムを参考にして作成するとよい（そのままコピーしてはダメです）。

Listing 1: 二分法のサンプルプログラム
12345678910111213141516171819202122232425262728293031#include < s t d l i b . h>#include <cmath>#include <i o s t r e a m >using namespace s t d ;i n t main ( ){double a = 1 . ;double b= = 4.;
double c = =1;
double d=4;
double x i n i t = = 3.;
double d x i n i t = 1 0 . ;i n t i t e r m a x =20;
double t h r e s h o l d =1e = 6;
double x1=x i n i t ;double f x 1=a * ( x1 * x1 * x1)+b * ( x1 * x1)+ c * x1+d ;double x2=x1+d x i n i t ;double f x 2=a * ( x2 * x2 * x2)+b * ( x2 * x2)+ c * x2+d ;i f ( f x 1 * f x 2 > 0. ) {cout<<”# The s i g n s o f f x 1 and f x 2 a r e i n d e n t i c a l . \ n” ;exit (0);};i f ( abs ( f x 1 )< t h r e s h o l d ) {cout<<x1<<” \n” ;

1

32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
4950515253
54
55
56

exit (0);};i f ( abs ( f x 2 )< t h r e s h o l d ) {cout<<x2<<” \n” ;exit (0);};
double x3 , f x 3 ;f o r ( i n t i =0; i <i t e r m a x ; i ++){
x3=(x1+x2 ) * 0 . 5 ;f x 3=a * ( x3 * x3 * x3)+b * ( x3 * x3)+ c * x3+d ;i f ( abs ( f x 3 )< t h r e s h o l d ) break ;i f ( f x 1 * f x 3 < 0. ) {x2=x3 ;
f x 2=f x 3 ;} else {
x1=x3 ;
f x 1=f x 3 ;};};cout<<x3<<” \n” ;};return 0 ;このプログラムは、ax3 + bx2 + cx + d = 0 の解を求めるためのものである。

x の最初の推定値を x init とし、その推定値における関数値と、推定値に dx init を加えた値における関数値の符号が異なる場合に、二分法により解を探索する。この条件が満たされない場合には強制終了となるため、x init とdx init を適切な値に変更する必要がある。また、二分法により求まる解もこれらの値に依存するため、複数の解を探索する場合は、これらの値を変更して対応することになる。二分法の繰り返し回数は変数 iter max により設定している。なお、計算途中に関数値がゼロに近い値になった場合（具体的には、変数 threshold よりも絶対値が小さくなった場合）、繰り返し計算を終了し、解を出力させている。絶対値を計算するために abs 関数を使っているが、環境によっては返り値が整数型になってしまう場合がある。このような場合は解析解を再現することが出来ないため、千葉まで相談すること。
次に、二分法よりも効率的に解を探索することができるニュートン法について説明する。
下図のように、ある値 x1 とその値における関数値 f (x1 ) を考える。ここで、x1 において f (x) の接線を引き、その接線が f (x) = 0 と交わる x を解の推定値とするのがニュートン法である。この解の推定値 x2 は、図に記載されているように、以下の式で求めることができる。

x2 = x1 − ∆x = x1 −f (x1 )f ′ (x1 )(1)解の推定値が真値から大きく離れている場合は、一次微係数に基づく推定は精度が悪いため、推定された解と真値にはまだ差が残る。そのため、この操作を繰り返し行う必要がある。繰り返しの過程で推定値が真値に近づくにつれて推定の精度は向上し、最終的に真値へ収束していくことが期待される。繰り返し計算は、前回の推定値と新しい推定値の相対差を調べ、その絶対値がある値未満になった場合に終了させるのが一般的である。なお、初期推定値によっては繰り返しによって推定値が発散する場合もあるため、繰り返しの最大回数を設定しておくとよいだろう（このような設定をしておかないと、いつまで経っても計算が終わらないという状況に陥る）。

Fig. 2: ニュートン法の考え方
二分法では、初期推定値として２点を設定する必要があったが、ニュートン法では 1 点のみ設定すればよいので、使い勝手が良いと言えるだろう。また、一般的に、真値への収束もニュートン法のほうが速い。

2

問題２：問題１の解をニュートン法により求め、解析解と比較せよ。


ニュートン法はニュートン · ラプソン法（Newton-Raphson method）とも呼ばれる方程式の逐次近似解法であり、以下の定理で保証される方法である（以下は文献1 の p. 239 から引用）。


区間 [a, b] において f (x) が 2 回微分可能で、f ′′ (x) > 0、また f (a) と f (b) とは符号が異なるとする。このとき区間 (a, b) 間にただ 1 つ f (x) = 0 の根（解）がある。


このとき、f (x0 ) > 0 である区間 [a, b] 内の x0 を初期値としてニュートン法を行うと、計算結果は解に漸近することが保証される。一方、この区間で f ′′ (x) < 0 であるときは、f (x0 ) < 0 である区間 [a, b] 内の x0 を初期値とすると、同様に計算結果は解に漸近する。

参考文献[1] 東京大学教養学部統計学教室編、「基礎統計学 III、自然科学の統計学」、東京大学出版会.

1 tokyo3