Finding the root of a function by Bisection Method

Question

The program has to look for a root in an interval [a,b]. The root should be declared with a certain accuracy eps. I.e it should look for a part-interval in [a,b], which has the length of eps.

The bisection algorithm should be:

1. Save the interval boundaries
2. Look if [a,b] has a root. (original given interval)
3. look if a-b < eps. If yes, part-interval found.
4. If no, divide [a,b] in half and continue with point 2. etc.

(We can assume that there is already a root in the given original interval [a,b])

I've tested it with following functions:

log x + 2

x*x*x

-3.*x*x + x + 2.

Are there any improvements possible or did I make a big mistake?

Bisection.cpp

#include "Bisection.hpp"
#include <cmath>

using namespace std;


void bisection(std::vector<double>& interval, double eps) {
  double a;
  double b;
  double c;

  a = interval[0];
  b = interval[1];

  if (b - a < eps) // interval already 'found'
  {
    return;
  }

  do
  {
    c = (a + b) / 2.0;

    if (f(c) == 0)
    {
      a = c;
      b = c; 
    }
    else 
    {
      if (f(a) * f(c) > 0)
      {
        a = c;
      }
      else
      {
        b = c;
      }                  
    }


  } while (f(c) != 0 || b - a > eps);

  /*found part-interval with the length of < eps*/
  interval[0] = a;
  interval[1] = b;
}

Answer 1

Headers and namespaces

• The <cmath> header is referenced, but never used; it can be removed.

• If we're to use std::vector, we require

  #include <vector>
  

• This is a danger sign:

  using namespace std;
  

  Bringing all names in from a namespace is problematic; namespace std particularly so. See Why is “using namespace std” considered bad practice?. Happily, we can just delete this line.

Function interface

• As Zeta's answer says, a vector is overkill to represent a numeric range.

• By preference, functions arguments should be inputs, and the output should be as a return value, something like:

  double find_root(double min, double max, double epsilon);
  

  It's also better to pass the function as an argument, rather than via the (undeclared) global f:

  double find_root(std::function<double(double)> f,
                   double min, double max, double epsilon);
  

  (You could pass a function pointer instead, but std::function is more general, and ought to have no overhead when wrapping a function pointer).

If the caller actually needs the lower and upper bounds of the improved range, you could return both values by using a suitable compound type. Choose from std::pair, std::tuple, std::array - or define your own "interval" type (which you could also use as the input range).

Implementation

• Instead of a do-while preceded by the same test, we could flip it into while-do:

  while (min + epsilon < max) { ... }
  

• We can avoid any risk of overflow in computing the midpoint if we divide first, then add:

   auto mid = 0.5 * min + 0.5 * max;
  

  If we want to move towards linear interpolation, we'd want

  auto mid = (min * f_max - max * f_min) / (f_max - f_min)
  

  and then move the furthest bound to the mid.

• We're calling f() repeatedly - this may be inefficient if it is a weighty computation. We probably want to keep a record of f(min), f(max) and f(mid), and only update as we narrow the range.

• Instead of multiplying two values to determine whether they have the same sign, it's more efficient to obtain each sign as a boolean, and compare those:

  if (f(min) * f(mid) > 0) // expensive
  if ((f(min) < 0) == (f(mid) < 0)) // cheaper
  

• There's little point ever comparing equality of a floating-point value; simply omit those tests and allow the function to converge.

The main() function

• No main() was supplied, but it was hinted at.

---

Improved code

#include <functional>

double find_root(std::function<double(double)> f,
                 double min, double max, double epsilon)
{
    auto f_min = f(min);
    while (min + epsilon < max) {
        auto const mid = 0.5 * min + 0.5 * max;
        auto const f_mid = f(mid);

        if ((f_min < 0) == (f_mid < 0)) {
            min = mid;
            f_min = f_mid;
        } else {
            max = mid;
        }
    }

    return min;
}

#include <cmath>
#include <iostream>
int main()
{
    std::cout << "log x + 2: "
              << find_root([](double x) { return std::log(x) + 2; }, 0.001, 1.0, 1e-7)
              << std::endl;
    std::cout << "x*x*x: "
              << find_root([](double x) { return x*x*x; }, -0.1, 1.0, 1e-11)
              << std::endl;
    std::cout << "-3.*x*x + x + 2.: "
              << find_root([](double x) { return -3.*x*x + x + 2.; }, -5, 5, 1e-7)
              << std::endl;
}

Answer 2

A std::vector<T> vec is an arbitrary collection of Ts. vec might contain none, one or many elements. However, you want to model an interval, which consists of exactly two elements. A std::vector<T> therefore isn't the correct container for an interval.

You could use std::pair<double, double> instead, or just write your own variant:

struct interval_type {
    double lower_bound;
    double upper_bound;
};

Next, try to initialize the variables instead of just declaring them. While we're at it, try to come up with better names, e.g. instead of

  double a;
  double b;
  double c;

  a = interval[0];
  b = interval[1];

write

  double lower = interval.lower_bound; // if you use the struct from above
  double upper = interval.upper_bound;

c's scope is too large, by the way. It's just necessary inside the loop, so we can limit its life time, but we need to break in this case:

  do
  {
    const double midpoint = (lower + upper) / 2.0;

    if (f(midpoint) == 0)
    {
      a = midpoint;
      b = midpoint; 
      break;
    }
    else 
    {
      if (f(lower) * f(midpoint) > 0)
      {
        lower = midpoint;
      }
      else
      {
        upper = midpoint;
      }                  
    }
  } while (upper - lower > eps);

Other than that, the test f(value) == 0 is usually too strict, since floating point arithmetic can be a finicky beast. std::abs(f(value)) <= some_eps works well enough most of the times.

Other than that, I suggest you to allow other functions than a fixed f. Something along

template <typename F>
interval_type bisection(F && f, interval_type interval, double eps = 1e-8){

}

comes to mind. As you can see, there is no reason to change the interval in place. Just return a new one.

Answer 3

Here are some additional test cases on the interval [0.0, 1.0] with epsilon = 0.1:

1. f(x) = (x - 0.2) * (x - 0.2)
2. f(x) = 1e-200 * (x - 0.8) + 1e-202
3. f(x) = x

Function #1 has a root within the interval, but the OP function bisection() fails because it does not check for a precondition of the bisection algorithm, namely that the function being evaluated is not non-zero and opposite signs at the interval endpoints. If the preconditions are not met, either return a separate error code / error condition, or return NAN for the returned containing end points.

Function #2 has a root at 0.8, but the OP function bisection() erroneously reports it as being in the interval [0.0, 0.0625] because the expression f(a) * f(b) is equal to 0 due to floating point underflow.

Function #3 has a root at 0.0, the end point. Note that using a condition of (f(a) < 0.0) == (f(b) < 0.0) fails this test case if checking the endpoints for zero function values is not done before entering the loop.

Unless you know the bisection() function will be used only in rather protected circumstances, it is a good idea to check its input parameters early in the function. Here is one example that passes the function f as a parameter, checks parameters for validity before continuing, avoids some other overflow exposures, avoids redundant calls to evaluate f (in case they are expensive) and returns an interval of length zero if an exact zero is found:

#include <cmath>

void bisection(double (*f)(double), double& lower, double& upper, double epsilon) {
  double f_lower;
  double f_upper;


  if (!(std::isfinite(lower) && std::isfinite(upper) && epsilon > 0.0 && lower <= upper))
  {
    lower = upper = NAN;
    return;
  }
  f_lower = f(lower);
  f_upper = f(upper);
  if (!(f_lower <= 0.0 && f_upper >= 0.0 || f_lower >= 0.0 && f_upper <= 0.0))
  {
    lower = upper = NAN;
    return;
  }

  if (f_lower == 0) {
    upper = lower;
  }
  if (f_upper == 0) {
    lower = upper;
  }

  while (lower + epsilon < upper)
  {
    double mid = lower * 0.5 + upper * 0.5;
    double f_mid = f(mid);

    if (f_mid == 0)
    {
      lower = mid;
      upper = mid; 
    }
    else 
    {
      if (f_lower > 0.0 && f_mid > 0.0 || f_lower < 0.0 && f_mid < 0.0)
      {
        lower = mid;
        f_lower = f_mid;
      }
      else
      {
        upper = mid;
      }                  
    }
  }
}