bisect.m

   1 function [root,iterations] = bisect(f, a, b, epsilon)
   2   ## Find a root of the function `f` on the closed interval [a,b].
   3   ##
   4   ## It is assumed that there are an odd number of roots within [a,b].
   5   ## The intermediate value theorem is used to determine whether or not
   6   ## there is a root on a given interval. So, for example, cos(-2) and
   7   ## cos(2) are both less than zero; the I.V.T. would not conclude that
   8   ## cos(x) has a root on -2<=x<=2.
   9   ##
  10   ## If `f` has more than one root on [a,b], the smallest root will be
  11   ## returned. This is an implementation detail.
  12   ##
  13   ##
  14   ## INPUTS:
  15   ##
  16   ##   * ``f`` - The function whose root we seek.
  17   ##
  18   ##   * ``a`` - The "left" endpoint of the interval in which we search.
  19   ##
  20   ##   * ``b`` - The "right" endpoint of the interval in which we search.
  21   ##
  22   ##   * ``epsilon`` - How close we should be to the actual root before we
  23   ##     halt the search and return the current approximation.
  24   ##
  25   ## OUTPUTS:
  26   ##
  27   ##   * ``root`` - The root that we found.
  28   ##
  29   ##   * ``iterations`` - The number of bisections that we performed
  30   ##     during the search.
  31   ##
  32
  33   ## Store these so we don't have to recompute them.
  34   fa = f(a);
  35   fb = f(b);
  36
  37   ## We first check for roots at a,b before bothering to bisect the
  38   ## interval.
  39   if (fa == 0)
  40     root = a;
  41     iterations = 0;
  42     return;
  43   end
  44
  45   if (fb == 0)
  46     root = b;
  47     iterations = 0;
  48     return;
  49   end
  50
  51   ## Bisect the interval.
  52   c = (a+b)/2;
  53   iterations = 1;
  54
  55   ## If we're within the prescribed tolerance, we're done.
  56   if (b-c < epsilon)
  57     root = c;
  58     return;
  59   end
  60
  61   if (has_root(fa,f(c)))
  62     [root,subiterations] = bisect(f,a,c,epsilon);
  63     iterations = iterations + subiterations;
  64   else
  65     [root,subiterations] = bisect(f,c,b,epsilon);
  66     iterations = iterations + subiterations;
  67   end
  68 end