3 二分法(bisection method)

3.1 計算方法

閉区間$[a,\,b]$で連続な関数$f(x)$の値が、

$$f(a)f(b)<0$$ (2)

ならば、 $f(\alpha)=0$となる$\alpha$が区間$[a,\,b]$にある。

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

実際にこの方法で

$$x^3-3x^2+9x-8=0$$ (3)

を計算した結果を図2に示す。この図より、$f(a)$と$f(b)$の関係 の式(2)を満たす区間$[a,b]$が1/2ずつ縮小していく様子がわかる。この 方法の長所と短所は、以下の通りである。

長所

閉区間$[a,b]$に解があれば、必ず解に収束する。間違いなく解 を探すので、ロバスト(robust：強靭な)な解法と言われている。 次に示すニュートン法とは異なり、連続であればどんな形の関数 でも解に収束するので信頼性が高い方法と言える。さらに、解の 精度も分かり便利である。解の誤差は、区間の幅$\vert b-a\vert$以下で る。

短所

収束が遅い(図6)。一次収束である。

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

3.2 アルゴリズム

関数はあらかじめ、プログラム中に書くものとする。更に、計算を打ち切る条 件もプログラム中に書くものとする。そうすると、図 3のような2分法のフローチャートが考えられる。

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

3.3 プログラム

このプログラムを暗記する必要はない。テストでは、アルゴリズム上、重要な部分を虫食 いにして出題するつもりである(たぶん)。 =1

#include <stdio.h>
double func(double x);

/*=============================================================*/
/*       main function                                         */
/*=============================================================*/

int main(){
  double eps=1e-15;         /* precision of calculation */
  double a, b, c;
  double test;
  char temp;
  int i=0;

  do{

    printf("\ninitial value a = ");
    scanf("%lf%c", &a, &temp);

    printf("initial value b = ");
    scanf("%lf%c", &b, &temp);

    test=func(a)*func(b);

    if(test >= 0){
      printf("   bad initial value !!  f(a)*f(b)>0\n\n");
    }

  }while(test >= 0);

  if(b-a<0){
    c=a;
    a=b;
    b=c;
  }

  while(b-a>eps){

    c=(a+b)/2;

    if(func(c)*func(a)<0){
      b=c;
    }else{
      a=c;
    }

    i++;
    printf(" %d\t%20.15f\n",i,c);

  }

  printf("\nsolution x = %20.15f\n\n",c);

  return(0);
}

/*=============================================================*/
/*       define function                                       */
/*=============================================================*/

double func(double x){
  double y;

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

  return(y);
}

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