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 -- | Use Newton's method to find a root of @f@ near the initial guess
  94 --   @x0@. If your guess is bad, this will recurse forever!
  95 newtons_method :: (Fractional a, Ord a)
  96                  => (a -> a) -- ^ The function @f@ whose root we seek
  97                  -> (a -> a) -- ^ The derivative of @f@
  98                  -> a       -- ^ The tolerance epsilon
  99                  -> a       -- ^ Initial guess, x-naught
 100                  -> Maybe a
 101 newtons_method f f' epsilon x0
 102   = find (\x -> abs (f x) < epsilon) x_n
 103   where
 104     x_n = newton_iterations f f' x0
 105
 106
 107
 108 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
 109 --   to the previous two values. Returns a list containing all
 110 --   generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
 111 --
 112 --   Examples:
 113 --
 114 --   >>> let fibs = iterate2 (+) 0 1
 115 --   >>> take 15 fibs
 116 --   [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
 117 --
 118 iterate2 :: (a -> a -> a) -- ^ The function @f@
 119          -> a           -- ^ The initial value @x0@
 120          -> a           -- ^ The second value, @x1@
 121          -> [a]         -- ^ The result list, [x0, x1, ...]
 122 iterate2 f x0 x1 =
 123   x0 : x1 : (go x0 x1)
 124   where
 125     go prev2 prev1 =
 126       let next = f prev2 prev1 in
 127         next : go prev1 next
 128
 129 -- | The sequence x_{n} of values obtained by applying the secant
 130 --   method on the function @f@ and initial guesses @x0@, @x1@.
 131 --
 132 --   The recursion more or less implements a two-parameter 'iterate',
 133 --   although one list is passed to the next iteration (as opposed to
 134 --   one function argument, with iterate). At each step, we peel the
 135 --   first two elements off the list and then compute/append elements
 136 --   three, four... onto the end of the list.
 137 --
 138 --   Examples:
 139 --
 140 --   Atkinson, p. 67.
 141 --   >>> let f x = x^6 - x - 1
 142 --   >>> take 4 $ secant_iterations f 2 1
 143 --   [2.0,1.0,1.0161290322580645,1.190577768676638]
 144 --
 145 secant_iterations :: (Fractional a, Ord a)
 146                     => (a -> a) -- ^ The function @f@ whose root we seek
 147                     -> a       -- ^ Initial guess, x-naught
 148                     -> a       -- ^ Second initial guess, x-one
 149                     -> [a]
 150 secant_iterations f x0 x1 =
 151   iterate2 g x0 x1
 152   where
 153   g prev2 prev1 =
 154     let x_change = prev1 - prev2
 155         y_change = (f prev1) - (f prev2)
 156     in
 157       (prev1 - (f prev1 * (x_change / y_change)))
 158
 159
 160 -- | Use the secant method to find a root of @f@ near the initial guesses
 161 --   @x0@ and @x1@. If your guesses are bad, this will recurse forever!
 162 --
 163 --   Examples:
 164 --
 165 --   Atkinson, p. 67.
 166 --   >>> let f x = x^6 - x - 1
 167 --   >>> let Just root = secant_method f (1/10^9) 2 1
 168 --   >>> root
 169 --   1.1347241384015196
 170 --   >>> abs (f root) < (1/10^9)
 171 --   True
 172 --
 173 secant_method :: (Fractional a, Ord a)
 174                  => (a -> a) -- ^ The function @f@ whose root we seek
 175                  -> a       -- ^ The tolerance epsilon
 176                  -> a       -- ^ Initial guess, x-naught
 177                  -> a       -- ^ Second initial guess, x-one
 178                  -> Maybe a
 179 secant_method f epsilon x0 x1
 180   = find (\x -> abs (f x) < epsilon) x_n
 181   where
 182     x_n = secant_iterations f x0 x1