1

私はこのとても素敵なフォーラムをしばらく休んでいます。数値解析コースを受講していて、二分法をプログラムするように求められました。これが私のコードです

/*
 * Bisection.cpp
 *
 *  Created on: 08/10/2012
 *  Author: BRabbit27
 *  École Polytechnique Fédérale de Lausanne - M.Sc. CSE
 */

#include <cmath>
#include <iostream>

using namespace std;

double functionA(double x) {
    return sin(2.0 * x) - 1.0 + x;
}

double functionB(double x) {
    return 2.0 * x / (1.0 + x / 1.5);
}

double bisectionMethod(double (*function)(double), double a, double b, double tolerance) {
    double x;
    double f;
    double error = tolerance + 1;
    int step = 0;
    double fa = (*function)(a);
    double fb = (*function)(b);
    //Check the conditions of a root in the given interval
    if (a < b) {
        if (fa * fb < 0) {
            while (error > tolerance) {
                step++;

                x = (a + b) / 2.0;
                f = (*function)(x);

                if (f == 0) {
                    cout << "Root found in x = " << x;
                    return x;
                } else if (f * fa > 0) {
                    a = x;
                } else if (f * fa < 0) {
                    b = x;
                }
                error = (b - a) / pow(2.0, (double) step + 1);
            }
            cout << "Root found in x = " << x;
            return x;
        } else {
            cout << "There not exist a root in that interval." << endl;
            return -1;
        }
    } else {
        cout << "Mandatory \"a < b\", verify." << endl;
        return -1;
    }
}

int main(int argc, char *argv[]){
    bisectionMethod(functionA, -3.0, 3.0, 10.0e-7);
}

私が抱えている唯一の問題は、x = 0.354492のときにルートが見つかり、実際のルートがx = 1/3にあるため、実際には倍精度または許容誤差のいずれかで問題が発生していることです。このコードを改善してより良い結果を得るにはどうすればよいかわかりません。何か案が？

4

1 に答える 1

2

本当の根はx=1/3ではありません！3.52288です

sin（2.0 * x）-1.0 + x = 0

sin（2.0 * x）= 1-x

1-1/3 = 2/3

sin（2/3）！= 2/3

あなたの寛容の定義は奇妙だと思います。許容範囲は、xあなたが受け入れる意志の範囲でなければなりませんよね？それは簡単です：b - a > tolerance。

double bisectionMethod(double (*function)(double), double a, double b, double tolerance) {
    if (a > b) {
        cout << "Mandatory \"a < b\", verify." << endl;
        return -1;
    }
    double fa = function(a);
    double fb = function(b);
    if (fa * fb > 0) {
        cout << "No change of sign - bijection not possible" << endl;
        return -1;
    }

    do {    
        double x = (a + b) / 2.0;
        double f = (*function)(x);
 
        if (f * fa > 0) {
            a = x;
        } else {
            b = x;
        }
    } while (b - a > tolerance);

    cout << "Root found in x = " << x;
    return x;
}
于 2012-10-08T18:33:06.100 に答える