3 二分法(bisection method)

3.1 計算方法

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

$\displaystyle 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ずつ縮小していく様子がわかる.

計算の終了は,

$\displaystyle 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軸との交点)に収束していく様子が分かる.
\includegraphics[keepaspectratio, scale=0.7]{figure/function_solution/NibunMethod.eps}

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: 二分法のフローチャート
\includegraphics[keepaspectratio, scale=0.8]{figure/flow_chart/flow_nibun.eps}



ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成18年7月11日


no counter