src/Roots/Simple.hs

   1 -- | The Roots.Simple module contains root-finding algorithms. That
   2 --   is, 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.Simple
  10 where
  11
  12 import qualified Roots.Fast as F
  13
  14
  15 -- | Does the (continuous) function @f@ have a root on the interval
  16 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
  17 --   [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
  18 --   and f(b) <= 0.
  19 --
  20 --   Examples:
  21 --
  22 --   >>> let f x = x**3
  23 --   >>> has_root f (-1) 1 Nothing
  24 --   True
  25 --
  26 --   This fails if we don't specify an @epsilon@, because cos(-2) ==
  27 --   cos(2) doesn't imply that there's a root on [-2,2].
  28 --
  29 --   >>> has_root cos (-2) 2 Nothing
  30 --   False
  31 --   >>> has_root cos (-2) 2 (Just 0.001)
  32 --   True
  33 --
  34 has_root :: (Fractional a, Ord a, Ord b, Num b)
  35          => (a -> b) -- ^ The function @f@
  36          -> a       -- ^ The \"left\" endpoint, @a@
  37          -> a       -- ^ The \"right\" endpoint, @b@
  38          -> Maybe a -- ^ The size of the smallest subinterval
  39                    --   we'll examine, @epsilon@
  40          -> Bool
  41 has_root f a b epsilon =
  42   F.has_root f a b epsilon Nothing Nothing
  43
  44
  45
  46
  47 -- | We are given a function @f@ and an interval [a,b]. The bisection
  48 --   method checks finds a root by splitting [a,b] in half repeatedly.
  49 --
  50 --   If one is found within some prescribed tolerance @epsilon@, it is
  51 --   returned. Otherwise, the interval [a,b] is split into two
  52 --   subintervals [a,c] and [c,b] of equal length which are then both
  53 --   checked via the same process.
  54 --
  55 --   Returns 'Just' the value x for which f(x) == 0 if one is found,
  56 --   or Nothing if one of the preconditions is violated.
  57 --
  58 --   Examples:
  59 --
  60 --   >>> bisect cos 1 2 0.001
  61 --   Just 1.5712890625
  62 --
  63 --   >>> bisect sin (-1) 1 0.001
  64 --   Just 0.0
  65 --
  66 bisect :: (Fractional a, Ord a, Num b, Ord b)
  67        => (a -> b) -- ^ The function @f@ whose root we seek
  68        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  69        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  70        -> a       -- ^ The tolerance, @epsilon@
  71        -> Maybe a
  72 bisect f a b epsilon =
  73   F.bisect f a b epsilon Nothing Nothing