bisect.m

   1 function root = 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   ## Store these so we don't have to recompute them.
  14   fa = f(a);
  15   fb = f(b);
  16
  17   ## We first check for roots at a,b before bothering to bisect the
  18   ## interval.
  19   if (fa == 0)
  20     root = a;
  21     return;
  22   end
  23
  24   if (fb == 0)
  25     root = b;
  26     return;
  27   end
  28
  29   ## Bisect the interval.
  30   c = (a+b)/2;
  31
  32   ## If we're within the prescribed tolerance, we're done.
  33   if (b-c < epsilon)
  34     root = c;
  35     return;
  36   end
  37
  38   if (has_root(fa,f(c)))
  39     root = bisect(f,a,c,epsilon);
  40   else
  41     root = bisect(f,c,b,epsilon);
  42   end
  43 end