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 fixed_point_error_ratios,
15 fixed_point_iteration_count,
22 import Data.List (find)
23 import NumericPrelude hiding ( abs )
24 import Algebra.Absolute ( abs )
25 import qualified Algebra.Additive as Additive ( C )
26 import qualified Algebra.Algebraic as Algebraic ( C )
27 import qualified Algebra.Field as Field ( C )
28 import qualified Algebra.RealField as RealField ( C )
29 import qualified Algebra.RealRing as RealRing ( C )
31 import Normed ( Normed(..) )
32 import qualified Roots.Fast as F (
34 fixed_point_iterations,
35 fixed_point_with_iterations,
40 -- | Does the (continuous) function @f@ have a root on the interval
41 -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
42 -- [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
48 -- >>> has_root f (-1) 1 Nothing
51 -- This fails if we don't specify an @epsilon@, because cos(-2) ==
52 -- cos(2) doesn't imply that there's a root on [-2,2].
54 -- >>> has_root cos (-2) 2 Nothing
56 -- >>> has_root cos (-2) 2 (Just 0.001)
59 has_root :: (RealField.C a, RealRing.C b)
60 => (a -> b) -- ^ The function @f@
61 -> a -- ^ The \"left\" endpoint, @a@
62 -> a -- ^ The \"right\" endpoint, @b@
63 -> Maybe a -- ^ The size of the smallest subinterval
64 -- we'll examine, @epsilon@
66 has_root f a b epsilon =
67 F.has_root f a b epsilon Nothing Nothing
70 -- | We are given a function @f@ and an interval [a,b]. The bisection
71 -- method finds a root by splitting [a,b] in half repeatedly.
73 -- If one is found within some prescribed tolerance @epsilon@, it is
74 -- returned. Otherwise, the interval [a,b] is split into two
75 -- subintervals [a,c] and [c,b] of equal length which are then both
76 -- checked via the same process.
78 -- Returns 'Just' the value x for which f(x) == 0 if one is found,
79 -- or Nothing if one of the preconditions is violated.
83 -- >>> let actual = 1.5707963267948966
84 -- >>> let Just root = bisect cos 1 2 0.001
87 -- >>> abs (root - actual) < 0.001
90 -- >>> bisect sin (-1) 1 0.001
93 bisect :: (RealField.C a, RealRing.C b)
94 => (a -> b) -- ^ The function @f@ whose root we seek
95 -> a -- ^ The \"left\" endpoint of the interval, @a@
96 -> a -- ^ The \"right\" endpoint of the interval, @b@
97 -> a -- ^ The tolerance, @epsilon@
99 bisect f a b epsilon =
100 F.bisect f a b epsilon Nothing Nothing
103 -- | We are given a function @f@ and an interval [a,b]. The trisection
104 -- method finds a root by splitting [a,b] into three
105 -- subintervals repeatedly.
107 -- If one is found within some prescribed tolerance @epsilon@, it is
108 -- returned. Otherwise, the interval [a,b] is split into two
109 -- subintervals [a,c] and [c,b] of equal length which are then both
110 -- checked via the same process.
112 -- Returns 'Just' the value x for which f(x) == 0 if one is found,
113 -- or Nothing if one of the preconditions is violated.
117 -- >>> let actual = 1.5707963267948966
118 -- >>> let Just root = trisect cos 1 2 0.001
120 -- 1.5713305898491083
121 -- >>> abs (root - actual) < 0.001
124 -- >>> trisect sin (-1) 1 0.001
127 trisect :: (RealField.C a, RealRing.C b)
128 => (a -> b) -- ^ The function @f@ whose root we seek
129 -> a -- ^ The \"left\" endpoint of the interval, @a@
130 -> a -- ^ The \"right\" endpoint of the interval, @b@
131 -> a -- ^ The tolerance, @epsilon@
133 trisect f a b epsilon =
134 F.trisect f a b epsilon Nothing Nothing
137 -- | Find a fixed point of the function @f@ with the search starting
138 -- at x0. We delegate to the version that returns the number of
139 -- iterations and simply discard the number of iterations.
141 fixed_point :: (Normed a, Additive.C a, Algebraic.C b, RealField.C b)
142 => (a -> a) -- ^ The function @f@ to iterate.
143 -> b -- ^ The tolerance, @epsilon@.
144 -> a -- ^ The initial value @x0@.
145 -> a -- ^ The fixed point.
146 fixed_point f epsilon x0 =
147 snd $ F.fixed_point_with_iterations f epsilon x0
150 -- | Return the number of iterations required to find a fixed point of
151 -- the function @f@ with the search starting at x0 and tolerance
152 -- @epsilon@. We delegate to the version that returns the number of
153 -- iterations and simply discard the fixed point.
154 fixed_point_iteration_count :: (Normed a,
158 => (a -> a) -- ^ The function @f@ to iterate.
159 -> b -- ^ The tolerance, @epsilon@.
160 -> a -- ^ The initial value @x0@.
161 -> Int -- ^ The fixed point.
162 fixed_point_iteration_count f epsilon x0 =
163 fst $ F.fixed_point_with_iterations f epsilon x0
166 -- | Returns a list of ratios,
168 -- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p}
170 -- of fixed point iterations for the function @f@ with initial guess
171 -- @x0@ and @p@ some positive power.
173 -- This is used to determine the rate of convergence.
175 fixed_point_error_ratios :: (Normed a,
179 => (a -> a) -- ^ The function @f@ to iterate.
180 -> a -- ^ The initial value @x0@.
181 -> a -- ^ The true solution, @x_star@.
182 -> Integer -- ^ The power @p@.
183 -> [b] -- ^ The resulting sequence of x_{n}.
184 fixed_point_error_ratios f x0 x_star p =
185 zipWith (/) en_plus_one en_exp
187 xn = F.fixed_point_iterations f x0
188 en = map (\x -> norm (x_star - x)) xn
189 en_plus_one = tail en
194 -- | The sequence x_{n} of values obtained by applying Newton's method
195 -- on the function @f@ and initial guess @x0@.
200 -- >>> let f x = x^6 - x - 1
201 -- >>> let f' x = 6*x^5 - 1
202 -- >>> tail $ take 4 $ newton_iterations f f' 2
203 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362]
205 newton_iterations :: (Field.C a)
206 => (a -> a) -- ^ The function @f@ whose root we seek
207 -> (a -> a) -- ^ The derivative of @f@
208 -> a -- ^ Initial guess, x-naught
210 newton_iterations f f' =
214 xn - ( (f xn) / (f' xn) )
217 -- | Use Newton's method to find a root of @f@ near the initial guess
218 -- @x0@. If your guess is bad, this will recurse forever!
224 -- >>> let f x = x^6 - x - 1
225 -- >>> let f' x = 6*x^5 - 1
226 -- >>> let Just root = newtons_method f f' (1/1000000) 2
228 -- 1.1347241385002211
229 -- >>> abs (f root) < 1/100000
232 -- >>> import Data.Number.BigFloat
233 -- >>> let eps = 1/(10^20) :: BigFloat Prec50
234 -- >>> let Just root = newtons_method f f' eps 2
236 -- 1.13472413840151949260544605450647284028100785303643e0
237 -- >>> abs (f root) < eps
240 newtons_method :: (RealField.C a)
241 => (a -> a) -- ^ The function @f@ whose root we seek
242 -> (a -> a) -- ^ The derivative of @f@
243 -> a -- ^ The tolerance epsilon
244 -> a -- ^ Initial guess, x-naught
246 newtons_method f f' epsilon x0 =
247 find (\x -> abs (f x) < epsilon) x_n
249 x_n = newton_iterations f f' x0
252 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
253 -- to the previous two values. Returns a list containing all
254 -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
258 -- >>> let fibs = iterate2 (+) 0 1
260 -- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
262 iterate2 :: (a -> a -> a) -- ^ The function @f@
263 -> a -- ^ The initial value @x0@
264 -> a -- ^ The second value, @x1@
265 -> [a] -- ^ The result list, [x0, x1, ...]
270 let next = f prev2 prev1 in
274 -- | The sequence x_{n} of values obtained by applying the secant
275 -- method on the function @f@ and initial guesses @x0@, @x1@.
277 -- The recursion more or less implements a two-parameter 'iterate',
278 -- although one list is passed to the next iteration (as opposed to
279 -- one function argument, with iterate). At each step, we peel the
280 -- first two elements off the list and then compute/append elements
281 -- three, four... onto the end of the list.
286 -- >>> let f x = x^6 - x - 1
287 -- >>> take 4 $ secant_iterations f 2 1
288 -- [2.0,1.0,1.0161290322580645,1.190577768676638]
290 secant_iterations :: (Field.C a)
291 => (a -> a) -- ^ The function @f@ whose root we seek
292 -> a -- ^ Initial guess, x-naught
293 -> a -- ^ Second initial guess, x-one
295 secant_iterations f =
299 let x_change = prev1 - prev2
300 y_change = (f prev1) - (f prev2)
302 (prev1 - (f prev1 * (x_change / y_change)))
305 -- | Use the secant method to find a root of @f@ near the initial guesses
306 -- @x0@ and @x1@. If your guesses are bad, this will recurse forever!
311 -- >>> let f x = x^6 - x - 1
312 -- >>> let Just root = secant_method f (1/10^9) 2 1
314 -- 1.1347241384015196
315 -- >>> abs (f root) < (1/10^9)
318 secant_method :: (RealField.C a)
319 => (a -> a) -- ^ The function @f@ whose root we seek
320 -> a -- ^ The tolerance epsilon
321 -> a -- ^ Initial guess, x-naught
322 -> a -- ^ Second initial guess, x-one
324 secant_method f epsilon x0 x1
325 = find (\x -> abs (f x) < epsilon) x_n
327 x_n = secant_iterations f x0 x1