C.2.2 二分法の場合

bisection.c

/*
 * bisection.c -- 二分法 (bisection method) で方程式 f(x)=0 を解く
 *   コンパイル: gcc -o bisection bisection.c -lm
 *   実行: ./bisection
 *   入力は "0 1 1e-15" など。ここで 1e-15 は 1かける 10 の -15 乗を意味する。
 */

#include <stdio.h>
#include <math.h>

double f(double x)
{
  return cos(x) - x;
}

int main()
{
  int i, maxitr = 100;
  double a, b, c, eps;
  double fa, fb, fc;

  printf(" 探す区間の左端, 右端, 許容精度ε: ");
  scanf("%lf %lf %lf", &a, &b, &eps);

  fa = f(a);
  fb = f(b);
  if (fa * fb > 0.0) {
    printf(" f(α) f(β) > 0 なのであきらめます。\n");
    exit(0);
  }
  else {
    printf("f(%20.15f)=%9.2e, f(%20.15f)=%9.2e\n", a, fa, b, fb);
    for (i = 0; i < maxitr; i++) {
      c = (a + b) / 2;
      fc = f(c);
      if (fc == 0.0)
        break;
      else if (fa * fc <= 0.0) {
        /* 左側 [a,c] に根がある */
        b = c;
        fb = fc;
      }
      else {
        /* 左側 [a,c] には根がないかもしれない。[c,b] にあるはず */
        a = c;
        fa = fc;
      }
      printf ("f(%20.15f)=%9.2e, f(%20.15f)=%9.2e\n", a, fa, b, fb);
      if ((b - a) <= eps)
        break;
    }
    printf ("二分法による近似解=%20.15f\n", c);
  }
  return 0;
}

区間 $(0,1)$ 内に解があることがわかるから、二分法で許容精度 $10^{-15}$ を指示して解かせたのが以下の結果 (入力は 0 1 1e-15)。関数値が、区間の左端では正、右端では負になったまま区間が縮小 して行くのを理解しよう。

isc-xas06% gcc -o bisection bisection.c -lm
isc-xas06% ./bisection 
 探す区間の左端α, 右端β, 許容精度ε=0 1 1e-15
f(   5.000000000000000e-01)= 3.78e-01, f(   1.000000000000000e+00)=-4.60e-01
f(   5.000000000000000e-01)= 3.78e-01, f(   7.500000000000000e-01)=-1.83e-02
f(   6.250000000000000e-01)= 1.86e-01, f(   7.500000000000000e-01)=-1.83e-02
f(   6.875000000000000e-01)= 8.53e-02, f(   7.500000000000000e-01)=-1.83e-02
f(   7.187500000000000e-01)= 3.39e-02, f(   7.500000000000000e-01)=-1.83e-02
f(   7.343750000000000e-01)= 7.87e-03, f(   7.500000000000000e-01)=-1.83e-02
f(   7.343750000000000e-01)= 7.87e-03, f(   7.421875000000000e-01)=-5.20e-03
f(   7.382812500000000e-01)= 1.35e-03, f(   7.421875000000000e-01)=-5.20e-03
中略
f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151627e-01)=-3.44e-15
f(   7.390851332151591e-01)= 2.55e-15, f(   7.390851332151627e-01)=-3.44e-15
f(   7.390851332151591e-01)= 2.55e-15, f(   7.390851332151609e-01)=-4.44e-16
f(   7.390851332151600e-01)= 1.11e-15, f(   7.390851332151609e-01)=-4.44e-16
f(   7.390851332151600e-01)= 1.11e-15
isc-xas06%