3 二分法(bisection method)

3.1 計算方法

二分法の原理は非常に単純であるが，場合によっては非常に強力な方法である．これ は，閉区間$[a,\,b]$で連続な関数$f(x)$の値が，

$$f(a)f(b)\leq 0$$ (10)

ならば， $f(\alpha)=0$となる$\alpha$が区間$[a,\,b]$にある--ということを使う．このことは，中間値の定理から保証されるが，常識的に考えてあたりまえである．た だ，連続な区間で適用できることを忘れてはならない．

実際の数値計算は， $f(a)f(b)\leq 0$であるような2点 $a,\,b(a\leq b)$から出発する．そして， 区間$[a,\,b]$を2分する点$c=(a+b)/2$に対して，$f(c)$を計算を行う． $f(c)f(a)\leq 0$な らば$b$を$c$と置き換え， $f(c)f(a)\geq 0$ならば$a$を$c$と置き換える．絶えず，区間 $[a,b]$の間に解があるようにするのである．この操作を繰り返して，区間の幅$\vert b-a\vert$が 与えられた値 $\varepsilon$ よりも小さくなったならば，計算を終了する．

解へ収束は収束率1/2の一次収束となる．1回の計算で，解の精度は1/2になるからである． 10回計算を行えば，はじめに比べて解の精度は $(1/2)^{10}=1/1024$となる．非常に良い精 度の近似解を得るためには，かなりの計算が必要である．人間と電卓では手間がかかって 大変である．このような単純反復作業はコンピューターの得意分野である．

実際にこの方法で式(4)を計算した結果を図 2に示す．この図より，$f(a)$と$f(b)$の関係の式 (10)を満たす区間$[a,b]$が1/2ずつ縮小していく様子がわかる．

計算の終了は，

$$b-a\leq\varepsilon$$ (11)

の条件を満たした場合とするのが一般的である．ここで， $\varepsilon$は近似解の精度 である．これを変えることにより，任意の精度で近似解を求めることができる．アルゴリ ズムから，$b-a$が大体の計算精度を表すことは自明である．

図 2: $f(x)=x^3-3x^2+9x-8$の実数解を二分法で解散し，その解の収束の 様子を示している．初期値は $a=-1,\,b=11$として，最初の解$c=x_0=5$が求 まり，順次より精度の良い $x_1,\,x_2,\,x_3,\cdots$が求まる．それが，解析解 $x=1.1659\cdots$(x軸との交点)に収束していく様子が分かる．

3.2 フローチャート

関数$f(x)$はあらかじめ，プログラム中に書くものとする．C言語の関数--サブルーチン-- を利用するのが良いだろう．具体的には，

double f(double x){
  double y;

  y=x*x*x-3*x*x+9*x-8;
  
  return y;
}

と書く．もちろんこれを使うためには，プロトタイプ宣言が必要である．また，計算を打 ち切る条件もプログラム中に書くものとする．通常，これは

#define EPS 1e-10

とプログラムの先頭に書く．このEPSが求めるべき解の精度を表す．マクロ定数なの で，普通，大文字--C言語の習慣--を使う．

図3のような二分法のフローチャートの通りにすれば，目的の 動作をするプログラムができる．

図 3: 二分法のフローチャート

ホームページ: Yamamoto's laboratory
著者: 山本昌志