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 Data.List (find)
  13
  14 import qualified Roots.Fast as F
  15
  16 -- | Does the (continuous) function @f@ have a root on the interval
  17 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
  18 --   [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
  19 --   and f(b) <= 0.
  20 --
  21 --   Examples:
  22 --
  23 --   >>> let f x = x**3
  24 --   >>> has_root f (-1) 1 Nothing
  25 --   True
  26 --
  27 --   This fails if we don't specify an @epsilon@, because cos(-2) ==
  28 --   cos(2) doesn't imply that there's a root on [-2,2].
  29 --
  30 --   >>> has_root cos (-2) 2 Nothing
  31 --   False
  32 --   >>> has_root cos (-2) 2 (Just 0.001)
  33 --   True
  34 --
  35 has_root :: (Fractional a, Ord a, Ord b, Num b)
  36          => (a -> b) -- ^ The function @f@
  37          -> a       -- ^ The \"left\" endpoint, @a@
  38          -> a       -- ^ The \"right\" endpoint, @b@
  39          -> Maybe a -- ^ The size of the smallest subinterval
  40                    --   we'll examine, @epsilon@
  41          -> Bool
  42 has_root f a b epsilon =
  43   F.has_root f a b epsilon Nothing Nothing
  44
  45
  46
  47
  48 -- | We are given a function @f@ and an interval [a,b]. The bisection
  49 --   method checks finds a root by splitting [a,b] in half repeatedly.
  50 --
  51 --   If one is found within some prescribed tolerance @epsilon@, it is
  52 --   returned. Otherwise, the interval [a,b] is split into two
  53 --   subintervals [a,c] and [c,b] of equal length which are then both
  54 --   checked via the same process.
  55 --
  56 --   Returns 'Just' the value x for which f(x) == 0 if one is found,
  57 --   or Nothing if one of the preconditions is violated.
  58 --
  59 --   Examples:
  60 --
  61 --   >>> bisect cos 1 2 0.001
  62 --   Just 1.5712890625
  63 --
  64 --   >>> bisect sin (-1) 1 0.001
  65 --   Just 0.0
  66 --
  67 bisect :: (Fractional a, Ord a, Num b, Ord b)
  68        => (a -> b) -- ^ The function @f@ whose root we seek
  69        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  70        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  71        -> a       -- ^ The tolerance, @epsilon@
  72        -> Maybe a
  73 bisect f a b epsilon =
  74   F.bisect f a b epsilon Nothing Nothing
  75
  76
  77
  78 -- | The sequence x_{n} of values obtained by applying Newton's method
  79 --   on the function @f@ and initial guess @x0@.
  80 newton_iterations :: (Fractional a, Ord a)
  81                     => (a -> a) -- ^ The function @f@ whose root we seek
  82                     -> (a -> a) -- ^ The derivative of @f@
  83                     -> a       -- ^ Initial guess, x-naught
  84                     -> [a]
  85 newton_iterations f f' x0 =
  86   iterate next x0
  87   where
  88   next xn =
  89     xn - ( (f xn) / (f' xn) )
  90
  91
  92
  93 newtons_method :: (Fractional a, Ord a)
  94                  => (a -> a) -- ^ The function @f@ whose root we seek
  95                  -> (a -> a) -- ^ The derivative of @f@
  96                  -> a       -- ^ The tolerance epsilon
  97                  -> a       -- ^ Initial guess, x-naught
  98                  -> Maybe a
  99 newtons_method f f' epsilon x0
 100   = find (\x -> abs (f x) < epsilon) x_n
 101   where
 102     x_n = newton_iterations f f' x0