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@.
201 -- >>> let f x = x^6 - x - 1
202 -- >>> let f' x = 6*x^5 - 1
203 -- >>> tail $ take 4 $ newton_iterations f f' 2
204 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362]
206 newton_iterations :: (Field.C a)
207 => (a -> a) -- ^ The function @f@ whose root we seek
208 -> (a -> a) -- ^ The derivative of @f@
209 -> a -- ^ Initial guess, x-naught
211 newton_iterations f f' =
215 xn - ( (f xn) / (f' xn) )
218 -- | Use Newton's method to find a root of @f@ near the initial guess
219 -- @x0@. If your guess is bad, this will recurse forever!
225 -- >>> let f x = x^6 - x - 1
226 -- >>> let f' x = 6*x^5 - 1
227 -- >>> let Just root = newtons_method f f' (1/1000000) 2
229 -- 1.1347241385002211
230 -- >>> abs (f root) < 1/100000
233 -- >>> import Data.Number.BigFloat
234 -- >>> let eps = 1/(10^20) :: BigFloat Prec50
235 -- >>> let Just root = newtons_method f f' eps 2
237 -- 1.13472413840151949260544605450647284028100785303643e0
238 -- >>> abs (f root) < eps
241 newtons_method :: (RealField.C a)
242 => (a -> a) -- ^ The function @f@ whose root we seek
243 -> (a -> a) -- ^ The derivative of @f@
244 -> a -- ^ The tolerance epsilon
245 -> a -- ^ Initial guess, x-naught
247 newtons_method f f' epsilon x0 =
248 find (\x -> abs (f x) < epsilon) x_n
250 x_n = newton_iterations f f' x0
253 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
254 -- to the previous two values. Returns a list containing all
255 -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
259 -- >>> let fibs = iterate2 (+) 0 1
261 -- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
263 iterate2 :: (a -> a -> a) -- ^ The function @f@
264 -> a -- ^ The initial value @x0@
265 -> a -- ^ The second value, @x1@
266 -> [a] -- ^ The result list, [x0, x1, ...]
271 let next = f prev2 prev1 in
275 -- | The sequence x_{n} of values obtained by applying the secant
276 -- method on the function @f@ and initial guesses @x0@, @x1@.
278 -- The recursion more or less implements a two-parameter 'iterate',
279 -- although one list is passed to the next iteration (as opposed to
280 -- one function argument, with iterate). At each step, we peel the
281 -- first two elements off the list and then compute/append elements
282 -- three, four... onto the end of the list.
288 -- >>> let f x = x^6 - x - 1
289 -- >>> take 4 $ secant_iterations f 2 1
290 -- [2.0,1.0,1.0161290322580645,1.190577768676638]
292 secant_iterations :: (Field.C a)
293 => (a -> a) -- ^ The function @f@ whose root we seek
294 -> a -- ^ Initial guess, x-naught
295 -> a -- ^ Second initial guess, x-one
297 secant_iterations f =
301 let x_change = prev1 - prev2
302 y_change = (f prev1) - (f prev2)
304 (prev1 - (f prev1 * (x_change / y_change)))
307 -- | Use the secant method to find a root of @f@ near the initial guesses
308 -- @x0@ and @x1@. If your guesses are bad, this will recurse forever!
314 -- >>> let f x = x^6 - x - 1
315 -- >>> let Just root = secant_method f (1/10^9) 2 1
317 -- 1.1347241384015196
318 -- >>> abs (f root) < (1/10^9)
321 secant_method :: (RealField.C a)
322 => (a -> a) -- ^ The function @f@ whose root we seek
323 -> a -- ^ The tolerance epsilon
324 -> a -- ^ Initial guess, x-naught
325 -> a -- ^ Second initial guess, x-one
327 secant_method f epsilon x0 x1
328 = find (\x -> abs (f x) < epsilon) x_n
330 x_n = secant_iterations f x0 x1