関数を区間を考えたとき, との符号が異なる場合(),区間内に少なくとも1つの解が存在する. 区間の中点での関数値を求め,その値の符号がと同じならば, 解は内に存在する. このとき,解の存在する区間はになっている. 同様にして,の中点での関数値を使うことで,解の存在区間はとなる. この作業を繰り返すことで方程式の解を求める方法を2分法(bisection method)と呼ぶ. 図1は2分法で近似解を求める手順を示している. 現在の解の存在区間をとする. はで初期化される. 中点をとする. 図1ではであるので,解はに存在するので, とする. 次の中点をであり, 図1ではであるので,解はに存在する. これを反復して,と求めていく. 2分法の手順を以下に示す.
反復を終了するための収束条件は以下となる. 2分法のコード例を以下に示す. /*! * 2分法(bisection method) * @param[in] func 関数値を与える関数ポインタ * @param[in] xl,xr 探索範囲 * @param[out] x 解 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す) * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) * @return 1:成功,0:失敗 */ int bisection(double func(const double), double xl, double xr, double &x, int &max_iter, double &eps) { double f = func(xl); double fmid = func(xr); // 探索範囲の境界での値の符号が異なる場合のみ探索 if(f*fmid >= 0.0) return 0.0; int k; double dx = fabs(xr-xl), xmid; for(k = 0; k < max_iter; ++k){ xmid = 0.5*(xl+xr); // 中点 dx *= 0.5; // 中点での関数値を求める fmid = func(xmid); // 収束判定 if(dx < eps || fmid == 0.0){ x = xmid; max_iter = k; eps = dx; return 1; } // 新しい区間 if(f*fmid < 0){ xr = xmid; } else{ xl = xmid; f = fmid; } } max_iter = k; eps = dx; return 0; } |