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 Algebra.Absolute (abs)
22 import qualified Algebra.Additive as Additive
23 import qualified Algebra.Algebraic as Algebraic
24 import qualified Algebra.Field as Field
25 import qualified Algebra.RealField as RealField
26 import qualified Algebra.RealRing as RealRing
28 -- | Does the (continuous) function @f@ have a root on the interval
29 -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in
30 -- [a,b] by the intermediate value theorem. Likewise when f(a) >= 0
36 -- >>> has_root f (-1) 1 Nothing
39 -- This fails if we don't specify an @epsilon@, because cos(-2) ==
40 -- cos(2) doesn't imply that there's a root on [-2,2].
42 -- >>> has_root cos (-2) 2 Nothing
44 -- >>> has_root cos (-2) 2 (Just 0.001)
47 has_root :: (RealField.C a, RealRing.C b)
48 => (a -> b) -- ^ The function @f@
49 -> a -- ^ The \"left\" endpoint, @a@
50 -> a -- ^ The \"right\" endpoint, @b@
51 -> Maybe a -- ^ The size of the smallest subinterval
52 -- we'll examine, @epsilon@
54 has_root f a b epsilon =
55 F.has_root f a b epsilon Nothing Nothing
58 -- | We are given a function @f@ and an interval [a,b]. The bisection
59 -- method finds a root by splitting [a,b] in half repeatedly.
61 -- If one is found within some prescribed tolerance @epsilon@, it is
62 -- returned. Otherwise, the interval [a,b] is split into two
63 -- subintervals [a,c] and [c,b] of equal length which are then both
64 -- checked via the same process.
66 -- Returns 'Just' the value x for which f(x) == 0 if one is found,
67 -- or Nothing if one of the preconditions is violated.
71 -- >>> let actual = 1.5707963267948966
72 -- >>> let Just root = bisect cos 1 2 0.001
75 -- >>> abs (root - actual) < 0.001
78 -- >>> bisect sin (-1) 1 0.001
81 bisect :: (RealField.C a, RealRing.C b)
82 => (a -> b) -- ^ The function @f@ whose root we seek
83 -> a -- ^ The \"left\" endpoint of the interval, @a@
84 -> a -- ^ The \"right\" endpoint of the interval, @b@
85 -> a -- ^ The tolerance, @epsilon@
87 bisect f a b epsilon =
88 F.bisect f a b epsilon Nothing Nothing
91 -- | We are given a function @f@ and an interval [a,b]. The trisection
92 -- method finds a root by splitting [a,b] into three
93 -- subintervals repeatedly.
95 -- If one is found within some prescribed tolerance @epsilon@, it is
96 -- returned. Otherwise, the interval [a,b] is split into two
97 -- subintervals [a,c] and [c,b] of equal length which are then both
98 -- checked via the same process.
100 -- Returns 'Just' the value x for which f(x) == 0 if one is found,
101 -- or Nothing if one of the preconditions is violated.
105 -- >>> let actual = 1.5707963267948966
106 -- >>> let Just root = trisect cos 1 2 0.001
108 -- 1.5713305898491083
109 -- >>> abs (root - actual) < 0.001
112 -- >>> trisect sin (-1) 1 0.001
115 trisect :: (RealField.C a, RealRing.C b)
116 => (a -> b) -- ^ The function @f@ whose root we seek
117 -> a -- ^ The \"left\" endpoint of the interval, @a@
118 -> a -- ^ The \"right\" endpoint of the interval, @b@
119 -> a -- ^ The tolerance, @epsilon@
121 trisect f a b epsilon =
122 F.trisect f a b epsilon Nothing Nothing
125 -- | Find a fixed point of the function @f@ with the search starting
126 -- at x0. We delegate to the version that returns the number of
127 -- iterations and simply discard the number of iterations.
129 fixed_point :: (Normed a, Additive.C a, Algebraic.C b, RealField.C b)
130 => (a -> a) -- ^ The function @f@ to iterate.
131 -> b -- ^ The tolerance, @epsilon@.
132 -> a -- ^ The initial value @x0@.
133 -> a -- ^ The fixed point.
134 fixed_point f epsilon x0 =
135 snd $ F.fixed_point_with_iterations f epsilon x0
138 -- | Return the number of iterations required to find a fixed point of
139 -- the function @f@ with the search starting at x0 and tolerance
140 -- @epsilon@. We delegate to the version that returns the number of
141 -- iterations and simply discard the fixed point.
142 fixed_point_iteration_count :: (Normed a,
146 => (a -> a) -- ^ The function @f@ to iterate.
147 -> b -- ^ The tolerance, @epsilon@.
148 -> a -- ^ The initial value @x0@.
149 -> Int -- ^ The fixed point.
150 fixed_point_iteration_count f epsilon x0 =
151 fst $ F.fixed_point_with_iterations f epsilon x0
154 -- | Returns a list of ratios,
156 -- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p}
158 -- of fixed point iterations for the function @f@ with initial guess
159 -- @x0@ and @p@ some positive power.
161 -- This is used to determine the rate of convergence.
163 fixed_point_error_ratios :: (Normed a,
167 => (a -> a) -- ^ The function @f@ to iterate.
168 -> a -- ^ The initial value @x0@.
169 -> a -- ^ The true solution, @x_star@.
170 -> Integer -- ^ The power @p@.
171 -> [b] -- ^ The resulting sequence of x_{n}.
172 fixed_point_error_ratios f x0 x_star p =
173 zipWith (/) en_plus_one en_exp
175 xn = F.fixed_point_iterations f x0
176 en = map (\x -> norm (x_star - x)) xn
177 en_plus_one = tail en
182 -- | The sequence x_{n} of values obtained by applying Newton's method
183 -- on the function @f@ and initial guess @x0@.
188 -- >>> let f x = x^6 - x - 1
189 -- >>> let f' x = 6*x^5 - 1
190 -- >>> tail $ take 4 $ newton_iterations f f' 2
191 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362]
193 newton_iterations :: (Field.C a)
194 => (a -> a) -- ^ The function @f@ whose root we seek
195 -> (a -> a) -- ^ The derivative of @f@
196 -> a -- ^ Initial guess, x-naught
198 newton_iterations f f' x0 =
202 xn - ( (f xn) / (f' xn) )
205 -- | Use Newton's method to find a root of @f@ near the initial guess
206 -- @x0@. If your guess is bad, this will recurse forever!
212 -- >>> let f x = x^6 - x - 1
213 -- >>> let f' x = 6*x^5 - 1
214 -- >>> let Just root = newtons_method f f' (1/1000000) 2
216 -- 1.1347241385002211
217 -- >>> abs (f root) < 1/100000
220 -- >>> import Data.Number.BigFloat
221 -- >>> let eps = 1/(10^20) :: BigFloat Prec50
222 -- >>> let Just root = newtons_method f f' eps 2
224 -- 1.13472413840151949260544605450647284028100785303643e0
225 -- >>> abs (f root) < eps
228 newtons_method :: (RealField.C a)
229 => (a -> a) -- ^ The function @f@ whose root we seek
230 -> (a -> a) -- ^ The derivative of @f@
231 -> a -- ^ The tolerance epsilon
232 -> a -- ^ Initial guess, x-naught
234 newtons_method f f' epsilon x0 =
235 find (\x -> abs (f x) < epsilon) x_n
237 x_n = newton_iterations f f' x0
240 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
241 -- to the previous two values. Returns a list containing all
242 -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
246 -- >>> let fibs = iterate2 (+) 0 1
248 -- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
250 iterate2 :: (a -> a -> a) -- ^ The function @f@
251 -> a -- ^ The initial value @x0@
252 -> a -- ^ The second value, @x1@
253 -> [a] -- ^ The result list, [x0, x1, ...]
258 let next = f prev2 prev1 in
262 -- | The sequence x_{n} of values obtained by applying the secant
263 -- method on the function @f@ and initial guesses @x0@, @x1@.
265 -- The recursion more or less implements a two-parameter 'iterate',
266 -- although one list is passed to the next iteration (as opposed to
267 -- one function argument, with iterate). At each step, we peel the
268 -- first two elements off the list and then compute/append elements
269 -- three, four... onto the end of the list.
274 -- >>> let f x = x^6 - x - 1
275 -- >>> take 4 $ secant_iterations f 2 1
276 -- [2.0,1.0,1.0161290322580645,1.190577768676638]
278 secant_iterations :: (Field.C a)
279 => (a -> a) -- ^ The function @f@ whose root we seek
280 -> a -- ^ Initial guess, x-naught
281 -> a -- ^ Second initial guess, x-one
283 secant_iterations f x0 x1 =
287 let x_change = prev1 - prev2
288 y_change = (f prev1) - (f prev2)
290 (prev1 - (f prev1 * (x_change / y_change)))
293 -- | Use the secant method to find a root of @f@ near the initial guesses
294 -- @x0@ and @x1@. If your guesses are bad, this will recurse forever!
299 -- >>> let f x = x^6 - x - 1
300 -- >>> let Just root = secant_method f (1/10^9) 2 1
302 -- 1.1347241384015196
303 -- >>> abs (f root) < (1/10^9)
306 secant_method :: (RealField.C a)
307 => (a -> a) -- ^ The function @f@ whose root we seek
308 -> a -- ^ The tolerance epsilon
309 -> a -- ^ Initial guess, x-naught
310 -> a -- ^ Second initial guess, x-one
312 secant_method f epsilon x0 x1
313 = find (\x -> abs (f x) < epsilon) x_n
315 x_n = secant_iterations f x0 x1