src/Roots/Simple.hs

   1 {-# LANGUAGE RebindableSyntax #-}
   2
   3 -- | The Roots.Simple module contains root-finding algorithms. That
   4 --   is, procedures to (numerically) find solutions to the equation,
   5 --
   6 --   > f(x) = 0
   7 --
   8 --   where f is assumed to be continuous on the interval of interest.
   9 --
  10
  11 module Roots.Simple
  12 where
  13
  14 import Data.List (find)
  15
  16 import Normed
  17
  18 import qualified Roots.Fast as F
  19
  20 import NumericPrelude hiding (abs)
  21 import qualified Algebra.Absolute as Absolute
  22 import Algebra.Absolute (abs)
  23 import qualified Algebra.Field as Field
  24 import qualified Algebra.RealField as RealField
  25 import qualified Algebra.RealRing as RealRing
  26
  27 -- | Does the (continuous) function @f@ have a root on the interval
  28 --   [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
  29 --   [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
  30 --   and f(b) <= 0.
  31 --
  32 --   Examples:
  33 --
  34 --   >>> let f x = x**3
  35 --   >>> has_root f (-1) 1 Nothing
  36 --   True
  37 --
  38 --   This fails if we don't specify an @epsilon@, because cos(-2) ==
  39 --   cos(2) doesn't imply that there's a root on [-2,2].
  40 --
  41 --   >>> has_root cos (-2) 2 Nothing
  42 --   False
  43 --   >>> has_root cos (-2) 2 (Just 0.001)
  44 --   True
  45 --
  46 has_root :: (RealField.C a, RealRing.C b)
  47          => (a -> b) -- ^ The function @f@
  48          -> a       -- ^ The \"left\" endpoint, @a@
  49          -> a       -- ^ The \"right\" endpoint, @b@
  50          -> Maybe a -- ^ The size of the smallest subinterval
  51                    --   we'll examine, @epsilon@
  52          -> Bool
  53 has_root f a b epsilon =
  54   F.has_root f a b epsilon Nothing Nothing
  55
  56
  57 -- | We are given a function @f@ and an interval [a,b]. The bisection
  58 --   method checks finds a root by splitting [a,b] in half repeatedly.
  59 --
  60 --   If one is found within some prescribed tolerance @epsilon@, it is
  61 --   returned. Otherwise, the interval [a,b] is split into two
  62 --   subintervals [a,c] and [c,b] of equal length which are then both
  63 --   checked via the same process.
  64 --
  65 --   Returns 'Just' the value x for which f(x) == 0 if one is found,
  66 --   or Nothing if one of the preconditions is violated.
  67 --
  68 --   Examples:
  69 --
  70 --   >>> bisect cos 1 2 0.001
  71 --   Just 1.5712890625
  72 --
  73 --   >>> bisect sin (-1) 1 0.001
  74 --   Just 0.0
  75 --
  76 bisect :: (RealField.C a, RealRing.C b)
  77        => (a -> b) -- ^ The function @f@ whose root we seek
  78        -> a       -- ^ The \"left\" endpoint of the interval, @a@
  79        -> a       -- ^ The \"right\" endpoint of the interval, @b@
  80        -> a       -- ^ The tolerance, @epsilon@
  81        -> Maybe a
  82 bisect f a b epsilon =
  83   F.bisect f a b epsilon Nothing Nothing
  84
  85
  86 -- | Find a fixed point of the function @f@ with the search starting
  87 --   at x0. We delegate to the version that returns the number of
  88 --   iterations and simply discard the number of iterations.
  89 --
  90 fixed_point :: (Normed a, RealField.C b)
  91             => (a -> a) -- ^ The function @f@ to iterate.
  92             -> b       -- ^ The tolerance, @epsilon@.
  93             -> a       -- ^ The initial value @x0@.
  94             -> a       -- ^ The fixed point.
  95 fixed_point f epsilon x0 =
  96   snd $ F.fixed_point_with_iterations f epsilon x0
  97
  98
  99 -- | Return the number of iterations required to find a fixed point of
 100 --   the function @f@ with the search starting at x0 and tolerance
 101 --   @epsilon@. We delegate to the version that returns the number of
 102 --   iterations and simply discard the fixed point.
 103 fixed_point_iteration_count :: (Normed a, RealField.C b)
 104                             => (a -> a) -- ^ The function @f@ to iterate.
 105                             -> b       -- ^ The tolerance, @epsilon@.
 106                             -> a       -- ^ The initial value @x0@.
 107                             -> Int       -- ^ The fixed point.
 108 fixed_point_iteration_count f epsilon x0 =
 109   fst $ F.fixed_point_with_iterations f epsilon x0
 110
 111
 112 -- | Returns a list of ratios,
 113 --
 114 --     ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p}
 115 --
 116 --   of fixed point iterations for the function @f@ with initial guess
 117 --   @x0@ and @p@ some positive power.
 118 --
 119 --   This is used to determine the rate of convergence.
 120 --
 121 fixed_point_error_ratios :: (Normed a, RealField.C b)
 122                    => (a -> a) -- ^ The function @f@ to iterate.
 123                    -> a       -- ^ The initial value @x0@.
 124                    -> a       -- ^ The true solution, @x_star@.
 125                    -> Integer -- ^ The power @p@.
 126                    -> [b]     -- ^ The resulting sequence of x_{n}.
 127 fixed_point_error_ratios f x0 x_star p =
 128   zipWith (/) en_plus_one en_exp
 129   where
 130     xn = F.fixed_point_iterations f x0
 131     en = map (\x -> norm (x_star - x)) xn
 132     en_plus_one = tail en
 133     en_exp = map (^p) en
 134
 135
 136
 137 -- | The sequence x_{n} of values obtained by applying Newton's method
 138 --   on the function @f@ and initial guess @x0@.
 139 --
 140 --   Examples:
 141 --
 142 --   Atkinson, p. 60.
 143 --   >>> let f x = x^6 - x - 1
 144 --   >>> let f' x = 6*x^5 - 1
 145 --   >>> tail $ take 4 $ newton_iterations f f' 2
 146 --   [1.6806282722513088,1.4307389882390624,1.2549709561094362]
 147 --
 148 newton_iterations :: (Field.C a)
 149                     => (a -> a) -- ^ The function @f@ whose root we seek
 150                     -> (a -> a) -- ^ The derivative of @f@
 151                     -> a       -- ^ Initial guess, x-naught
 152                     -> [a]
 153 newton_iterations f f' x0 =
 154   iterate next x0
 155   where
 156   next xn =
 157     xn - ( (f xn) / (f' xn) )
 158
 159
 160 -- | Use Newton's method to find a root of @f@ near the initial guess
 161 --   @x0@. If your guess is bad, this will recurse forever!
 162 --
 163 --   Examples:
 164 --
 165 --   Atkinson, p. 60.
 166 --
 167 --   >>> let f x = x^6 - x - 1
 168 --   >>> let f' x = 6*x^5 - 1
 169 --   >>> let Just root = newtons_method f f' (1/1000000) 2
 170 --   >>> root
 171 --   1.1347241385002211
 172 --   >>> abs (f root) < 1/100000
 173 --   True
 174 --
 175 --   >>> import Data.Number.BigFloat
 176 --   >>> let eps = 1/(10^20) :: BigFloat Prec50
 177 --   >>> let Just root = newtons_method f f' eps 2
 178 --   >>> root
 179 --   1.13472413840151949260544605450647284028100785303643e0
 180 --   >>> abs (f root) < eps
 181 --   True
 182 --
 183 newtons_method :: (RealField.C a)
 184                  => (a -> a) -- ^ The function @f@ whose root we seek
 185                  -> (a -> a) -- ^ The derivative of @f@
 186                  -> a       -- ^ The tolerance epsilon
 187                  -> a       -- ^ Initial guess, x-naught
 188                  -> Maybe a
 189 newtons_method f f' epsilon x0 =
 190   find (\x -> abs (f x) < epsilon) x_n
 191   where
 192     x_n = newton_iterations f f' x0
 193
 194
 195 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
 196 --   to the previous two values. Returns a list containing all
 197 --   generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
 198 --
 199 --   Examples:
 200 --
 201 --   >>> let fibs = iterate2 (+) 0 1
 202 --   >>> take 15 fibs
 203 --   [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
 204 --
 205 iterate2 :: (a -> a -> a) -- ^ The function @f@
 206          -> a           -- ^ The initial value @x0@
 207          -> a           -- ^ The second value, @x1@
 208          -> [a]         -- ^ The result list, [x0, x1, ...]
 209 iterate2 f x0 x1 =
 210   x0 : x1 : (go x0 x1)
 211   where
 212     go prev2 prev1 =
 213       let next = f prev2 prev1 in
 214         next : go prev1 next
 215
 216
 217 -- | The sequence x_{n} of values obtained by applying the secant
 218 --   method on the function @f@ and initial guesses @x0@, @x1@.
 219 --
 220 --   The recursion more or less implements a two-parameter 'iterate',
 221 --   although one list is passed to the next iteration (as opposed to
 222 --   one function argument, with iterate). At each step, we peel the
 223 --   first two elements off the list and then compute/append elements
 224 --   three, four... onto the end of the list.
 225 --
 226 --   Examples:
 227 --
 228 --   Atkinson, p. 67.
 229 --   >>> let f x = x^6 - x - 1
 230 --   >>> take 4 $ secant_iterations f 2 1
 231 --   [2.0,1.0,1.0161290322580645,1.190577768676638]
 232 --
 233 secant_iterations :: (Field.C a)
 234                     => (a -> a) -- ^ The function @f@ whose root we seek
 235                     -> a       -- ^ Initial guess, x-naught
 236                     -> a       -- ^ Second initial guess, x-one
 237                     -> [a]
 238 secant_iterations f x0 x1 =
 239   iterate2 g x0 x1
 240   where
 241   g prev2 prev1 =
 242     let x_change = prev1 - prev2
 243         y_change = (f prev1) - (f prev2)
 244     in
 245       (prev1 - (f prev1 * (x_change / y_change)))
 246
 247
 248 -- | Use the secant method to find a root of @f@ near the initial guesses
 249 --   @x0@ and @x1@. If your guesses are bad, this will recurse forever!
 250 --
 251 --   Examples:
 252 --
 253 --   Atkinson, p. 67.
 254 --   >>> let f x = x^6 - x - 1
 255 --   >>> let Just root = secant_method f (1/10^9) 2 1
 256 --   >>> root
 257 --   1.1347241384015196
 258 --   >>> abs (f root) < (1/10^9)
 259 --   True
 260 --
 261 secant_method :: (RealField.C a)
 262                  => (a -> a) -- ^ The function @f@ whose root we seek
 263                  -> a       -- ^ The tolerance epsilon
 264                  -> a       -- ^ Initial guess, x-naught
 265                  -> a       -- ^ Second initial guess, x-one
 266                  -> Maybe a
 267 secant_method f epsilon x0 x1
 268   = find (\x -> abs (f x) < epsilon) x_n
 269   where
 270     x_n = secant_iterations f x0 x1