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 checks 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 -- >>> bisect cos 1 2 0.001
74 -- >>> bisect sin (-1) 1 0.001
77 bisect :: (RealField.C a, RealRing.C b)
78 => (a -> b) -- ^ The function @f@ whose root we seek
79 -> a -- ^ The \"left\" endpoint of the interval, @a@
80 -> a -- ^ The \"right\" endpoint of the interval, @b@
81 -> a -- ^ The tolerance, @epsilon@
83 bisect f a b epsilon =
84 F.bisect f a b epsilon Nothing Nothing
87 -- | Find a fixed point of the function @f@ with the search starting
88 -- at x0. We delegate to the version that returns the number of
89 -- iterations and simply discard the number of iterations.
91 fixed_point :: (Normed a, Additive.C a, Algebraic.C b, RealField.C b)
92 => (a -> a) -- ^ The function @f@ to iterate.
93 -> b -- ^ The tolerance, @epsilon@.
94 -> a -- ^ The initial value @x0@.
95 -> a -- ^ The fixed point.
96 fixed_point f epsilon x0 =
97 snd $ F.fixed_point_with_iterations f epsilon x0
100 -- | Return the number of iterations required to find a fixed point of
101 -- the function @f@ with the search starting at x0 and tolerance
102 -- @epsilon@. We delegate to the version that returns the number of
103 -- iterations and simply discard the fixed point.
104 fixed_point_iteration_count :: (Normed a,
108 => (a -> a) -- ^ The function @f@ to iterate.
109 -> b -- ^ The tolerance, @epsilon@.
110 -> a -- ^ The initial value @x0@.
111 -> Int -- ^ The fixed point.
112 fixed_point_iteration_count f epsilon x0 =
113 fst $ F.fixed_point_with_iterations f epsilon x0
116 -- | Returns a list of ratios,
118 -- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p}
120 -- of fixed point iterations for the function @f@ with initial guess
121 -- @x0@ and @p@ some positive power.
123 -- This is used to determine the rate of convergence.
125 fixed_point_error_ratios :: (Normed a,
129 => (a -> a) -- ^ The function @f@ to iterate.
130 -> a -- ^ The initial value @x0@.
131 -> a -- ^ The true solution, @x_star@.
132 -> Integer -- ^ The power @p@.
133 -> [b] -- ^ The resulting sequence of x_{n}.
134 fixed_point_error_ratios f x0 x_star p =
135 zipWith (/) en_plus_one en_exp
137 xn = F.fixed_point_iterations f x0
138 en = map (\x -> norm (x_star - x)) xn
139 en_plus_one = tail en
144 -- | The sequence x_{n} of values obtained by applying Newton's method
145 -- on the function @f@ and initial guess @x0@.
150 -- >>> let f x = x^6 - x - 1
151 -- >>> let f' x = 6*x^5 - 1
152 -- >>> tail $ take 4 $ newton_iterations f f' 2
153 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362]
155 newton_iterations :: (Field.C a)
156 => (a -> a) -- ^ The function @f@ whose root we seek
157 -> (a -> a) -- ^ The derivative of @f@
158 -> a -- ^ Initial guess, x-naught
160 newton_iterations f f' x0 =
164 xn - ( (f xn) / (f' xn) )
167 -- | Use Newton's method to find a root of @f@ near the initial guess
168 -- @x0@. If your guess is bad, this will recurse forever!
174 -- >>> let f x = x^6 - x - 1
175 -- >>> let f' x = 6*x^5 - 1
176 -- >>> let Just root = newtons_method f f' (1/1000000) 2
178 -- 1.1347241385002211
179 -- >>> abs (f root) < 1/100000
182 -- >>> import Data.Number.BigFloat
183 -- >>> let eps = 1/(10^20) :: BigFloat Prec50
184 -- >>> let Just root = newtons_method f f' eps 2
186 -- 1.13472413840151949260544605450647284028100785303643e0
187 -- >>> abs (f root) < eps
190 newtons_method :: (RealField.C a)
191 => (a -> a) -- ^ The function @f@ whose root we seek
192 -> (a -> a) -- ^ The derivative of @f@
193 -> a -- ^ The tolerance epsilon
194 -> a -- ^ Initial guess, x-naught
196 newtons_method f f' epsilon x0 =
197 find (\x -> abs (f x) < epsilon) x_n
199 x_n = newton_iterations f f' x0
202 -- | Takes a function @f@ of two arguments and repeatedly applies @f@
203 -- to the previous two values. Returns a list containing all
204 -- generated values, f(x0, x1), f(x1, x2), f(x2, x3)...
208 -- >>> let fibs = iterate2 (+) 0 1
210 -- [0,1,1,2,3,5,8,13,21,34,55,89,144,233,377]
212 iterate2 :: (a -> a -> a) -- ^ The function @f@
213 -> a -- ^ The initial value @x0@
214 -> a -- ^ The second value, @x1@
215 -> [a] -- ^ The result list, [x0, x1, ...]
220 let next = f prev2 prev1 in
224 -- | The sequence x_{n} of values obtained by applying the secant
225 -- method on the function @f@ and initial guesses @x0@, @x1@.
227 -- The recursion more or less implements a two-parameter 'iterate',
228 -- although one list is passed to the next iteration (as opposed to
229 -- one function argument, with iterate). At each step, we peel the
230 -- first two elements off the list and then compute/append elements
231 -- three, four... onto the end of the list.
236 -- >>> let f x = x^6 - x - 1
237 -- >>> take 4 $ secant_iterations f 2 1
238 -- [2.0,1.0,1.0161290322580645,1.190577768676638]
240 secant_iterations :: (Field.C a)
241 => (a -> a) -- ^ The function @f@ whose root we seek
242 -> a -- ^ Initial guess, x-naught
243 -> a -- ^ Second initial guess, x-one
245 secant_iterations f x0 x1 =
249 let x_change = prev1 - prev2
250 y_change = (f prev1) - (f prev2)
252 (prev1 - (f prev1 * (x_change / y_change)))
255 -- | Use the secant method to find a root of @f@ near the initial guesses
256 -- @x0@ and @x1@. If your guesses are bad, this will recurse forever!
261 -- >>> let f x = x^6 - x - 1
262 -- >>> let Just root = secant_method f (1/10^9) 2 1
264 -- 1.1347241384015196
265 -- >>> abs (f root) < (1/10^9)
268 secant_method :: (RealField.C a)
269 => (a -> a) -- ^ The function @f@ whose root we seek
270 -> a -- ^ The tolerance epsilon
271 -> a -- ^ Initial guess, x-naught
272 -> a -- ^ Second initial guess, x-one
274 secant_method f epsilon x0 x1
275 = find (\x -> abs (f x) < epsilon) x_n
277 x_n = secant_iterations f x0 x1