From: Michael Orlitzky Date: Thu, 13 Sep 2012 18:00:51 +0000 (-0400) Subject: Add two initial octave functions, has_root and bisect. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=78254b47cdaf395eac1e31845a2cc9edb55a4aca;p=octave.git Add two initial octave functions, has_root and bisect. --- 78254b47cdaf395eac1e31845a2cc9edb55a4aca 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