Bisection method lost of precision

Question

I've been absent from this very nice forum for a while. I'm taking a Numerical Analysis course and I was asked to program the bisection method, here's my code

/*
 * 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);
}

The only problem I have is that the root is found when x = 0.354492 and the real root is in x=1/3 so actually either I have something bad with double precision or with my tolerance. I don't know how can I improve this code to have a better result. Any idea?

Answer 1

The real root is not x = 1/3! It's 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

---

Your definition of tolerance strikes me as odd. Tolerance should be the range in x you're will to accept, right? Well that's easy: 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;
}