2.6 制御 (2) while 文, if 文 -- 2分法による方程式の解法

bisection.cpp

   1 /*
   2  * bisection.cpp -- 二分法 (bisection method) で方程式 f(x)=0 を解く
   3  *   コンパイル: g++ -o bisection bisection.cpp -lm
   4  *   実行: ./bisection
   5  */
   6 
   7 #include <iostream>
   8 #include <cstdlib> // exit()
   9 #include <iomanip>
  10 #include <cmath>
  11 
  12 using namespace std;
  13 
  14 int main(void)
  15 {
  16     int i, maxitr = 100;
  17     double alpha, beta, a, b, c, eps;
  18     double fa, fb, fc;
  19     double f(double);
  20 
  21     cout.setf(ios::scientific, ios::floatfield);
  22 
  23     cout << " 探す区間の左端α, 右端β, 許容精度ε=";
  24     cin >> alpha >> beta >> eps;
  25 
  26     a = alpha; b = beta;
  27     fa = f(a); fb = f(b);
  28 
  29     if (fa * fb > 0.0) {
  30         cout << " f(α) f(β) > 0 なのであきらめます。" << endl;
  31         exit(0);
  32     }
  33     else {
  34         for (i = 0; i < maxitr; i++) {
  35             c = (a + b) / 2; fc = f(c);
  36             if (fc == 0.0)
  37                 break;
  38             else if (fa * fc <= 0.0) {
  39                 /* 左側 [a,c] に根がある */
  40                 b = c; fb = fc;
  41             } else {
  42                 /* 左側 [a,c] には根がないかもしれない。[c,b] にあるはず */
  43                 a = c; fa = fc;
  44             }
  45             cout << "f(" << setw(24) << setprecision(15) << a << ")="
  46                          << setw(9) << setprecision(2) << fa
  47                  << ", f(" << setw(24) << setprecision(15) << b << ")="
  48                            << setw(9) << setprecision(2) << fb << endl;
  49             if ((b - a) <= eps)
  50                 break;
  51         }
  52         cout << "f(" << setw(24) << setprecision(15) << c << ")="
  53              << setw(9) << setprecision(2) << fc << endl;
  54     }
  55     return 0;
  56 }
  57 
  58 double f(double x)
  59 {
  60     return cos(x) - x;
  61 }

コンパイルと実行結果

   1 oyabun% g++ -i bisection bisection.cpp
   2 oyabun% ./bisection
   3  探す区間の左端α, 右端β, 許容精度ε=0 1 1e-14
   4 f(   5.000000000000000e-01)= 3.78e-01, f(   1.000000000000000e+00)=-4.60e-01
   5 f(   5.000000000000000e-01)= 3.78e-01, f(   7.500000000000000e-01)=-1.83e-02
   6 f(   6.250000000000000e-01)= 1.86e-01, f(   7.500000000000000e-01)=-1.83e-02
   7 f(   6.875000000000000e-01)= 8.53e-02, f(   7.500000000000000e-01)=-1.83e-02
   8 f(   7.187500000000000e-01)= 3.39e-02, f(   7.500000000000000e-01)=-1.83e-02
   9 f(   7.343750000000000e-01)= 7.87e-03, f(   7.500000000000000e-01)=-1.83e-02
  10 f(   7.343750000000000e-01)= 7.87e-03, f(   7.421875000000000e-01)=-5.20e-03
  11 f(   7.382812500000000e-01)= 1.35e-03, f(   7.421875000000000e-01)=-5.20e-03
  12 f(   7.382812500000000e-01)= 1.35e-03, f(   7.402343750000000e-01)=-1.92e-03
  13 f(   7.382812500000000e-01)= 1.35e-03, f(   7.392578125000000e-01)=-2.89e-04
  14 f(   7.387695312500000e-01)= 5.28e-04, f(   7.392578125000000e-01)=-2.89e-04
  15 f(   7.390136718750000e-01)= 1.20e-04, f(   7.392578125000000e-01)=-2.89e-04
  16 f(   7.390136718750000e-01)= 1.20e-04, f(   7.391357421875000e-01)=-8.47e-05
  17 f(   7.390747070312500e-01)= 1.74e-05, f(   7.391357421875000e-01)=-8.47e-05
  18 f(   7.390747070312500e-01)= 1.74e-05, f(   7.391052246093750e-01)=-3.36e-05
  19 f(   7.390747070312500e-01)= 1.74e-05, f(   7.390899658203125e-01)=-8.09e-06
  20 f(   7.390823364257812e-01)= 4.68e-06, f(   7.390899658203125e-01)=-8.09e-06
  21 f(   7.390823364257812e-01)= 4.68e-06, f(   7.390861511230469e-01)=-1.70e-06
  22 f(   7.390842437744141e-01)= 1.49e-06, f(   7.390861511230469e-01)=-1.70e-06
  23 f(   7.390842437744141e-01)= 1.49e-06, f(   7.390851974487305e-01)=-1.08e-07
  24 f(   7.390847206115723e-01)= 6.91e-07, f(   7.390851974487305e-01)=-1.08e-07
  25 f(   7.390849590301514e-01)= 2.92e-07, f(   7.390851974487305e-01)=-1.08e-07
  26 f(   7.390850782394409e-01)= 9.20e-08, f(   7.390851974487305e-01)=-1.08e-07
  27 f(   7.390850782394409e-01)= 9.20e-08, f(   7.390851378440857e-01)=-7.75e-09
  28 f(   7.390851080417633e-01)= 4.21e-08, f(   7.390851378440857e-01)=-7.75e-09
  29 f(   7.390851229429245e-01)= 1.72e-08, f(   7.390851378440857e-01)=-7.75e-09
  30 f(   7.390851303935051e-01)= 4.72e-09, f(   7.390851378440857e-01)=-7.75e-09
  31 f(   7.390851303935051e-01)= 4.72e-09, f(   7.390851341187954e-01)=-1.51e-09
  32 f(   7.390851322561502e-01)= 1.61e-09, f(   7.390851341187954e-01)=-1.51e-09
  33 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851341187954e-01)=-1.51e-09
  34 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851336531341e-01)=-7.33e-10
  35 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851334203035e-01)=-3.43e-10
  36 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851333038881e-01)=-1.48e-10
  37 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851332456805e-01)=-5.11e-11
  38 f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851332165767e-01)=-2.37e-12
  39 f(   7.390851332020247e-01)= 2.20e-11, f(   7.390851332165767e-01)=-2.37e-12
  40 f(   7.390851332093007e-01)= 9.81e-12, f(   7.390851332165767e-01)=-2.37e-12
  41 f(   7.390851332129387e-01)= 3.72e-12, f(   7.390851332165767e-01)=-2.37e-12
  42 f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332165767e-01)=-2.37e-12
  43 f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332156672e-01)=-8.48e-13
  44 f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332152124e-01)=-8.66e-14
  45 f(   7.390851332149850e-01)= 2.94e-13, f(   7.390851332152124e-01)=-8.66e-14
  46 f(   7.390851332150987e-01)= 1.04e-13, f(   7.390851332152124e-01)=-8.66e-14
  47 f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332152124e-01)=-8.66e-14
  48 f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151840e-01)=-3.91e-14
  49 f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151698e-01)=-1.53e-14
  50 f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151627e-01)=-3.44e-15
  51 f(   7.390851332151627e-01)=-3.44e-15
  52 oyabun%

桂田 祐史
2018-11-04