1 {-# LANGUAGE RebindableSyntax #-}
3 -- | The Roots.Simple module contains root-finding algorithms. That
4 -- is, procedures to (numerically) find solutions to the equation,
8 -- where f is assumed to be continuous on the interval of interest.
14 import Data.List (find)
18 import qualified Roots.Fast as F
20 import NumericPrelude hiding (abs)
21 import qualified Algebra.Absolute as Absolute
22 import Algebra.Absolute (abs)
23 import qualified Algebra.Additive as Additive
24 import qualified Algebra.Algebraic as Algebraic
25 import qualified Algebra.Field as Field
26 import qualified Algebra.RealField as RealField
27 import qualified Algebra.RealRing as RealRing
29 -- | Does the (continuous) function @f@ have a root on the interval
30 -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
31 -- [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
37 -- >>> has_root f (-1) 1 Nothing
40 -- This fails if we don't specify an @epsilon@, because cos(-2) ==
41 -- cos(2) doesn't imply that there's a root on [-2,2].
43 -- >>> has_root cos (-2) 2 Nothing
45 -- >>> has_root cos (-2) 2 (Just 0.001)
48 has_root :: (RealField.C a, RealRing.C b)
49 => (a -> b) -- ^ The function @f@
50 -> a -- ^ The \"left\" endpoint, @a@
51 -> a -- ^ The \"right\" endpoint, @b@
52 -> Maybe a -- ^ The size of the smallest subinterval
53 -- we'll examine, @epsilon@
55 has_root f a b epsilon =
56 F.has_root f a b epsilon Nothing Nothing
59 -- | We are given a function @f@ and an interval [a,b]. The bisection
60 -- method checks finds a root by splitting [a,b] in half repeatedly.
62 -- If one is found within some prescribed tolerance @epsilon@, it is
63 -- returned. Otherwise, the interval [a,b] is split into two
64 -- subintervals [a,c] and [c,b] of equal length which are then both
65 -- checked via the same process.
67 -- Returns 'Just' the value x for which f(x) == 0 if one is found,
68 -- or Nothing if one of the preconditions is violated.
72 -- >>> bisect cos 1 2 0.001
75 -- >>> bisect sin (-1) 1 0.001
78 bisect :: (RealField.C a, RealRing.C b)
79 => (a -> b) -- ^ The function @f@ whose root we seek
80 -> a -- ^ The \"left\" endpoint of the interval, @a@
81 -> a -- ^ The \"right\" endpoint of the interval, @b@
82 -> a -- ^ The tolerance, @epsilon@
84 bisect f a b epsilon =
85 F.bisect f a b epsilon Nothing Nothing
88 -- | Find a fixed point of the function @f@ with the search starting
89 -- at x0. We delegate to the version that returns the number of
90 -- iterations and simply discard the number of iterations.
92 fixed_point :: (Normed a, Algebraic.C a, Algebraic.C b, RealField.C b)
93 => (a -> a) -- ^ The function @f@ to iterate.
94 -> b -- ^ The tolerance, @epsilon@.
95 -> a -- ^ The initial value @x0@.
96 -> a -- ^ The fixed point.
97 fixed_point f epsilon x0 =
98 snd $ F.fixed_point_with_iterations f epsilon x0
101 -- | Return the number of iterations required to find a fixed point of
102 -- the function @f@ with the search starting at x0 and tolerance
103 -- @epsilon@. We delegate to the version that returns the number of
104 -- iterations and simply discard the fixed point.
105 fixed_point_iteration_count :: (Normed a,
109 => (a -> a) -- ^ The function @f@ to iterate.
110 -> b -- ^ The tolerance, @epsilon@.
111 -> a -- ^ The initial value @x0@.
112 -> Int -- ^ The fixed point.
113 fixed_point_iteration_count f epsilon x0 =
114 fst $ F.fixed_point_with_iterations f epsilon x0
117 -- | Returns a list of ratios,
119 -- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p}
121 -- of fixed point iterations for the function @f@ with initial guess
122 -- @x0@ and @p@ some positive power.
124 -- This is used to determine the rate of convergence.
126 fixed_point_error_ratios :: (Normed a,
130 => (a -> a) -- ^ The function @f@ to iterate.
131 -> a -- ^ The initial value @x0@.
132 -> a -- ^ The true solution, @x_star@.
133 -> Integer -- ^ The power @p@.
134 -> [b] -- ^ The resulting sequence of x_{n}.
135 fixed_point_error_ratios f x0 x_star p =
136 zipWith (/) en_plus_one en_exp
138 xn = F.fixed_point_iterations f x0
139 en = map (\x -> norm (x_star - x)) xn
140 en_plus_one = tail en
145 -- | The sequence x_{n} of values obtained by applying Newton's method
146 -- on the function @f@ and initial guess @x0@.
151 -- >>> let f x = x^6 - x - 1
152 -- >>> let f' x = 6*x^5 - 1
153 -- >>> tail $ take 4 $ newton_iterations f f' 2
154 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362]
156 newton_iterations :: (Field.C a)
157 => (a -> a) -- ^ The function @f@ whose root we seek
158 -> (a -> a) -- ^ The derivative of @f@
159 -> a -- ^ Initial guess, x-naught
161 newton_iterations f f' x0 =
165 xn - ( (f xn) / (f' xn) )
168 -- | Use Newton's method to find a root of @f@ near the initial guess
169 -- @x0@. If your guess is bad, this will recurse forever!
175 -- >>> let f x = x^6 - x - 1
176 -- >>> let f' x = 6*x^5 - 1
177 -- >>> let Just root = newtons_method f f' (1/1000000) 2
179 -- 1.1347241385002211
180 -- >>> abs (f root) < 1/100000
183 -- >>> import Data.Number.BigFloat
184 -- >>> let eps = 1/(10^20) :: BigFloat Prec50
185 -- >>> let Just root = newtons_method f f' eps 2
187 -- 1.13472413840151949260544605450647284028100785303643e0
188 -- >>> abs (f root) < eps
191 newtons_method :: (RealField.C a)
192 => (a -> a) -- ^ The function @f@ whose root we seek
193 -> (a -> a) -- ^ The derivative of @f@
194 -> a -- ^ The tolerance epsilon
195 -> a -- ^ Initial guess, x-naught
197 newtons_method f f' epsilon x0 =
198 find (\x -> abs (f x) < epsilon) x_n
200 x_n = newton_iterations f f' x0
203 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
204 -- to the previous two values. Returns a list containing all
205 -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
209 -- >>> let fibs = iterate2 (+) 0 1
211 -- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
213 iterate2 :: (a -> a -> a) -- ^ The function @f@
214 -> a -- ^ The initial value @x0@
215 -> a -- ^ The second value, @x1@
216 -> [a] -- ^ The result list, [x0, x1, ...]
221 let next = f prev2 prev1 in
225 -- | The sequence x_{n} of values obtained by applying the secant
226 -- method on the function @f@ and initial guesses @x0@, @x1@.
228 -- The recursion more or less implements a two-parameter 'iterate',
229 -- although one list is passed to the next iteration (as opposed to
230 -- one function argument, with iterate). At each step, we peel the
231 -- first two elements off the list and then compute/append elements
232 -- three, four... onto the end of the list.
237 -- >>> let f x = x^6 - x - 1
238 -- >>> take 4 $ secant_iterations f 2 1
239 -- [2.0,1.0,1.0161290322580645,1.190577768676638]
241 secant_iterations :: (Field.C a)
242 => (a -> a) -- ^ The function @f@ whose root we seek
243 -> a -- ^ Initial guess, x-naught
244 -> a -- ^ Second initial guess, x-one
246 secant_iterations f x0 x1 =
250 let x_change = prev1 - prev2
251 y_change = (f prev1) - (f prev2)
253 (prev1 - (f prev1 * (x_change / y_change)))
256 -- | Use the secant method to find a root of @f@ near the initial guesses
257 -- @x0@ and @x1@. If your guesses are bad, this will recurse forever!
262 -- >>> let f x = x^6 - x - 1
263 -- >>> let Just root = secant_method f (1/10^9) 2 1
265 -- 1.1347241384015196
266 -- >>> abs (f root) < (1/10^9)
269 secant_method :: (RealField.C a)
270 => (a -> a) -- ^ The function @f@ whose root we seek
271 -> a -- ^ The tolerance epsilon
272 -> a -- ^ Initial guess, x-naught
273 -> a -- ^ Second initial guess, x-one
275 secant_method f epsilon x0 x1
276 = find (\x -> abs (f x) < epsilon) x_n
278 x_n = secant_iterations f x0 x1