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

bisection.cpp

/*
 * bisection.cpp -- 二分法 (bisection method) で方程式 f(x)=0 を解く
 *   コンパイル: g++ -o bisection bisection.cpp -lm
 *   実行: ./bisection
 */

#include <iostream>
#include <cstdlib> // exit()
#include <iomanip>
#include <cmath>

using namespace std;

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

    cout.setf(ios::scientific, ios::floatfield);

    cout << " 探す区間の左端α, 右端β, 許容精度ε=";
    cin >> alpha >> beta >> eps;

    a = alpha; b = beta;
    fa = f(a); fb = f(b);

    if (fa * fb > 0.0) {
        cout << " f(α) f(β) > 0 なのであきらめます。" << endl;
        exit(0);
    }
    else {
        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;
            }
            cout << "f(" << setw(24) << setprecision(15) << a << ")="
                         << setw(9) << setprecision(2) << fa
                 << ", f(" << setw(24) << setprecision(15) << b << ")="
                           << setw(9) << setprecision(2) << fb << endl;
            if ((b - a) <= eps)
                break;
        }
        cout << "f(" << setw(24) << setprecision(15) << c << ")="
             << setw(9) << setprecision(2) << fc << endl;
    }
    return 0;
}

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

このプログラムの実行結果は以下のようになる。

コンパイルと実行結果

oyabun% g++ -i bisection bisection.cpp
oyabun% ./bisection
 探す区間の左端α, 右端β, 許容精度ε=0 1 1e-14
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.382812500000000e-01)= 1.35e-03, f(   7.402343750000000e-01)=-1.92e-03
f(   7.382812500000000e-01)= 1.35e-03, f(   7.392578125000000e-01)=-2.89e-04
f(   7.387695312500000e-01)= 5.28e-04, f(   7.392578125000000e-01)=-2.89e-04
f(   7.390136718750000e-01)= 1.20e-04, f(   7.392578125000000e-01)=-2.89e-04
f(   7.390136718750000e-01)= 1.20e-04, f(   7.391357421875000e-01)=-8.47e-05
f(   7.390747070312500e-01)= 1.74e-05, f(   7.391357421875000e-01)=-8.47e-05
f(   7.390747070312500e-01)= 1.74e-05, f(   7.391052246093750e-01)=-3.36e-05
f(   7.390747070312500e-01)= 1.74e-05, f(   7.390899658203125e-01)=-8.09e-06
f(   7.390823364257812e-01)= 4.68e-06, f(   7.390899658203125e-01)=-8.09e-06
f(   7.390823364257812e-01)= 4.68e-06, f(   7.390861511230469e-01)=-1.70e-06
f(   7.390842437744141e-01)= 1.49e-06, f(   7.390861511230469e-01)=-1.70e-06
f(   7.390842437744141e-01)= 1.49e-06, f(   7.390851974487305e-01)=-1.08e-07
f(   7.390847206115723e-01)= 6.91e-07, f(   7.390851974487305e-01)=-1.08e-07
f(   7.390849590301514e-01)= 2.92e-07, f(   7.390851974487305e-01)=-1.08e-07
f(   7.390850782394409e-01)= 9.20e-08, f(   7.390851974487305e-01)=-1.08e-07
f(   7.390850782394409e-01)= 9.20e-08, f(   7.390851378440857e-01)=-7.75e-09
f(   7.390851080417633e-01)= 4.21e-08, f(   7.390851378440857e-01)=-7.75e-09
f(   7.390851229429245e-01)= 1.72e-08, f(   7.390851378440857e-01)=-7.75e-09
f(   7.390851303935051e-01)= 4.72e-09, f(   7.390851378440857e-01)=-7.75e-09
f(   7.390851303935051e-01)= 4.72e-09, f(   7.390851341187954e-01)=-1.51e-09
f(   7.390851322561502e-01)= 1.61e-09, f(   7.390851341187954e-01)=-1.51e-09
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851341187954e-01)=-1.51e-09
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851336531341e-01)=-7.33e-10
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851334203035e-01)=-3.43e-10
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851333038881e-01)=-1.48e-10
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851332456805e-01)=-5.11e-11
f(   7.390851331874728e-01)= 4.63e-11, f(   7.390851332165767e-01)=-2.37e-12
f(   7.390851332020247e-01)= 2.20e-11, f(   7.390851332165767e-01)=-2.37e-12
f(   7.390851332093007e-01)= 9.81e-12, f(   7.390851332165767e-01)=-2.37e-12
f(   7.390851332129387e-01)= 3.72e-12, f(   7.390851332165767e-01)=-2.37e-12
f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332165767e-01)=-2.37e-12
f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332156672e-01)=-8.48e-13
f(   7.390851332147577e-01)= 6.74e-13, f(   7.390851332152124e-01)=-8.66e-14
f(   7.390851332149850e-01)= 2.94e-13, f(   7.390851332152124e-01)=-8.66e-14
f(   7.390851332150987e-01)= 1.04e-13, f(   7.390851332152124e-01)=-8.66e-14
f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332152124e-01)=-8.66e-14
f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151840e-01)=-3.91e-14
f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151698e-01)=-1.53e-14
f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151627e-01)=-3.44e-15
f(   7.390851332151627e-01)=-3.44e-15
oyabun%