From 78254b47cdaf395eac1e31845a2cc9edb55a4aca Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 13 Sep 2012 14:00:51 -0400
Subject: [PATCH] Add two initial octave functions, has_root and bisect.

---
 bisect.m   | 45 +++++++++++++++++++++++++++++++++++++++++++++
 has_root.m | 24 ++++++++++++++++++++++++
 2 files changed, 69 insertions(+)
 create mode 100755 bisect.m
 create mode 100755 has_root.m

diff --git a/bisect.m b/bisect.m
new file mode 100755
index 0000000..56e099e
--- /dev/null
+++ b/bisect.m
@@ -0,0 +1,45 @@
+#!/usr/bin/octave --silent
+
+function root = bisect(f, a, b, epsilon)
+  ## Find a root of the function `f` on the closed interval [a,b].
+  ##
+  ## It is assumed that there are an odd number of roots within [a,b].
+  ## The intermediate value theorem is used to determine whether or not
+  ## there is a root on a given interval. So, for example, cos(-2) and
+  ## cos(2) are both less than zero; the I.V.T. would not conclude that
+  ## cos(x) has a root on -2<=x<=2.
+  ##
+  ## If `f` has more than one root on [a,b], the smallest root will be
+  ## returned. This is an implementation detail.
+
+  ## Store these so we don't have to recompute them.
+  fa = f(a);
+  fb = f(b);
+
+  ## We first check for roots at a,b before bothering to bisect the
+  ## interval.  
+  if (fa == 0)
+    root = a;
+    return;
+  end
+
+  if (fb == 0)
+    root = b;
+    return;
+  end
+
+  ## Bisect the interval.
+  c = (a+b)/2;
+  
+  ## If we're within the prescribed tolerance, we're done.
+  if (b-c < epsilon)
+    root = c;
+    return;
+  end
+
+  if (has_root(fa,f(c)))
+    root = bisect(f,a,c,epsilon);
+  else
+    root = bisect(f,c,b,epsilon);
+  end
+end
diff --git a/has_root.m b/has_root.m
new file mode 100755
index 0000000..5c835d8
--- /dev/null
+++ b/has_root.m
@@ -0,0 +1,24 @@
+#!/usr/bin/octave --silent
+
+function has_root = has_root(fa, fb)
+  ## Use the intermediate value theorem to determine whether or not some
+  ## function has an odd number of roots on an interval. If the function
+  ## in question has an even number of roots, the result will be
+  ## incorrect.
+  ##
+  ## Call the function whose roots we're concerned with 'f'. The two
+  ## parameters `fa` and `fb` should correspond to f(a) and f(b).
+  ##
+
+  ## If either f(a) or f(b) is zero, the product of their signs will be
+  ## zero and either a or b is a root. If the product of their signs is
+  ## negative, then f(a) and f(b) are non-zero and have opposite sign,
+  ## so there must be a root on (a,b). The only case we don't want is
+  ## when f(a) and f(b) have the same sign; in this case, the product of
+  ## their signs would be one.
+  if (sign(fa) * sign(fb) != 1)
+    has_root = true;
+  else
+    has_root = false;
+  end
+end
-- 
2.43.2