bisect.m

   1 #!/usr/bin/octave --silent
   2
   3 function root = bisect(f, a, b, epsilon)
   4   ## Find a root of the function `f` on the closed interval [a,b].
   5   ##
   6   ## It is assumed that there are an odd number of roots within [a,b].
   7   ## The intermediate value theorem is used to determine whether or not
   8   ## there is a root on a given interval. So, for example, cos(-2) and
   9   ## cos(2) are both less than zero; the I.V.T. would not conclude that
  10   ## cos(x) has a root on -2<=x<=2.
  11   ##
  12   ## If `f` has more than one root on [a,b], the smallest root will be
  13   ## returned. This is an implementation detail.
  14
  15   ## Store these so we don't have to recompute them.
  16   fa = f(a);
  17   fb = f(b);
  18
  19   ## We first check for roots at a,b before bothering to bisect the
  20   ## interval.
  21   if (fa == 0)
  22     root = a;
  23     return;
  24   end
  25
  26   if (fb == 0)
  27     root = b;
  28     return;
  29   end
  30
  31   ## Bisect the interval.
  32   c = (a+b)/2;
  33
  34   ## If we're within the prescribed tolerance, we're done.
  35   if (b-c < epsilon)
  36     root = c;
  37     return;
  38   end
  39
  40   if (has_root(fa,f(c)))
  41     root = bisect(f,a,c,epsilon);
  42   else
  43     root = bisect(f,c,b,epsilon);
  44   end
  45 end