src/Roots.hs

   1 -- | The Roots module contains root-finding algorithms. That is,
   2 --   procedures to (numerically) find solutions to the equation,
   3 --
   4 --   > f(x) = 0
   5 --
   6 --   where f is assumed to be continuous on the interval of interest.
   7 --
   8
   9 module Roots
  10 where
  11
  12
  13 -- | Does the (continuous) function @f@ have a root on the interval
  14 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
  15 --   [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
  16 --   and f(b) <= 0.
  17 --
  18 --   Examples:
  19 --
  20 --   >>> let f x = x**3
  21 --   >>> has_root f (-1) 1 Nothing
  22 --   True
  23 --
  24 --   This fails if we don't specify an @epsilon@, because cos(-2) ==
  25 --   cos(2) doesn't imply that there's a root on [-2,2].
  26 --
  27 --   >>> has_root cos (-2) 2 Nothing
  28 --   False
  29 --   >>> has_root cos (-2) 2 (Just 0.001)
  30 --   True
  31 --
  32 has_root :: (Fractional a, Ord a, Ord b, Num b)
  33          => (a -> b) -- ^ The function @f@
  34          -> a       -- ^ The \"left\" endpoint, @a@
  35          -> a       -- ^ The \"right\" endpoint, @b@
  36          -> Maybe a -- ^ The size of the smallest subinterval
  37                    --   we'll examine, @epsilon@
  38          -> Bool
  39 has_root f a b epsilon =
  40   if not ((signum (f a)) * (signum (f b)) == 1) then
  41     -- We don't care about epsilon here, there's definitely a root!
  42     True
  43   else
  44     if (b - a) <= epsilon' then
  45       -- Give up, return false.
  46       False
  47     else
  48       -- If either [a,c] or [c,b] have roots, we do too.
  49       (has_root f a c (Just epsilon')) || (has_root f c b (Just epsilon'))
  50   where
  51     -- If the size of the smallest subinterval is not specified,
  52     -- assume we just want to check once on all of [a,b].
  53     epsilon' = case epsilon of
  54                  Nothing -> (b-a)
  55                  Just eps -> eps
  56     c = (a + b)/2
  57
  58
  59
  60 -- | We are given a function @f@ and an interval [a,b]. The bisection
  61 --   method checks finds a root by splitting [a,b] in half repeatedly.
  62 --
  63 --   If one is found within some prescribed tolerance @epsilon@, it is
  64 --   returned. Otherwise, the interval [a,b] is split into two
  65 --   subintervals [a,c] and [c,b] of equal length which are then both
  66 --   checked via the same process.
  67 --
  68 --   Returns 'Just' the value x for which f(x) == 0 if one is found,
  69 --   or Nothing if one of the preconditions is violated.
  70 --
  71 --   Examples:
  72 --
  73 --   >>> bisect cos 1 2 0.001
  74 --   Just 1.5712890625
  75 --
  76 --   >>> bisect sin (-1) 1 0.001
  77 --   Just 0.0
  78 --
  79 bisect :: (Fractional a, Ord a, Num b, Ord b)
  80        => (a -> b) -- ^ The function @f@ whose root we seek
  81        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  82        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  83        -> a       -- ^ The tolerance, @epsilon@
  84        -> Maybe a
  85 bisect f a b epsilon
  86   -- We pass @epsilon@ to the 'has_root' function because if we want a
  87   -- result within epsilon of the true root, we need to know that
  88   -- there *is* a root within an interval of length epsilon.
  89   | not (has_root f a b (Just epsilon)) = Nothing
  90   | f a == 0 = Just a
  91   | f b == 0 = Just b
  92   | (b - c) < epsilon = Just c
  93   | otherwise =
  94       if (has_root f a c (Just epsilon)) then bisect f a c epsilon
  95       else bisect f c b epsilon
  96   where
  97     c = (a + b) / 2