From: Michael Orlitzky Date: Wed, 17 Oct 2012 16:41:12 +0000 (-0400) Subject: Add an iteration count to the fixed_point function, rename it, and move it to the... X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=c5da1efa77844ae6159dfc781ed886fdffbbf4d1;p=numerical-analysis.git Add an iteration count to the fixed_point function, rename it, and move it to the Roots.Fast module. Implement fixed_point and fixed_point_iteration_count in Roots.Simple in terms of the general function. Add a fixed_point_error_ratios function. --- diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs index 8e49750..d6e4d7c 100644 --- a/src/Roots/Fast.hs +++ b/src/Roots/Fast.hs @@ -6,6 +6,11 @@ module Roots.Fast where +import Data.List (find) + +import Vector + + has_root :: (Fractional a, Ord a, Ord b, Num b) => (a -> b) -- ^ The function @f@ -> a -- ^ The \"left\" endpoint, @a@ @@ -80,3 +85,52 @@ bisect f a b epsilon f_of_a f_of_b Just v -> v c = (a + b) / 2 + + + +-- | Iterate the function @f@ with the initial guess @x0@ in hopes of +-- finding a fixed point. +fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate. + -> a -- ^ The initial value @x0@. + -> [a] -- ^ The resulting sequence of x_{n}. +fixed_point_iterations f x0 = + iterate f x0 + + +-- | Find a fixed point of the function @f@ with the search starting +-- at x0. This will find the first element in the chain f(x0), +-- f(f(x0)),... such that the magnitude of the difference between it +-- and the next element is less than epsilon. +-- +-- We also return the number of iterations required. +-- +fixed_point_with_iterations :: (Vector a, RealFrac b) + => (a -> a) -- ^ The function @f@ to iterate. + -> b -- ^ The tolerance, @epsilon@. + -> a -- ^ The initial value @x0@. + -> (Int, a) -- ^ The (iterations, fixed point) pair +fixed_point_with_iterations f epsilon x0 = + (fst winning_pair) + where + xn = fixed_point_iterations f x0 + xn_plus_one = tail xn + + abs_diff v w = norm (v - w) + + -- The nth entry in this list is the absolute value of x_{n} - + -- x_{n+1}. + differences = zipWith abs_diff xn xn_plus_one + + -- This produces the list [(n, xn)] so that we can determine + -- the number of iterations required. + numbered_xn = zip [0..] xn + + -- A list of pairs, (xn, |x_{n} - x_{n+1}|). + pairs = zip numbered_xn differences + + -- The pair (xn, |x_{n} - x_{n+1}|) with + -- |x_{n} - x_{n+1}| < epsilon. The pattern match on 'Just' is + -- "safe" since the list is infinite. We'll succeed or loop + -- forever. + Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs + diff --git a/src/Roots/Simple.hs b/src/Roots/Simple.hs index 5aed7a1..3237d60 100644 --- a/src/Roots/Simple.hs +++ b/src/Roots/Simple.hs @@ -215,40 +215,51 @@ secant_method f epsilon x0 x1 -fixed_point_iterations :: (a -> a) -- ^ The function @f@ to iterate. - -> a -- ^ The initial value @x0@. - -> [a] -- ^ The resulting sequence of x_{n}. -fixed_point_iterations f x0 = - iterate f x0 - - -- | Find a fixed point of the function @f@ with the search starting --- at x0. This will find the first element in the chain f(x0), --- f(f(x0)),... such that the magnitude of the difference between it --- and the next element is less than epsilon. +-- at x0. We delegate to the version that returns the number of +-- iterations and simply discard the number of iterations. -- -fixed_point :: (Num a, Vector a, RealFrac b) +fixed_point :: (Vector a, RealFrac b) => (a -> a) -- ^ The function @f@ to iterate. -> b -- ^ The tolerance, @epsilon@. -> a -- ^ The initial value @x0@. -> a -- ^ The fixed point. fixed_point f epsilon x0 = - fst winning_pair - where - xn = fixed_point_iterations f x0 - xn_plus_one = tail $ fixed_point_iterations f x0 + snd $ F.fixed_point_with_iterations f epsilon x0 - abs_diff v w = norm (v - w) - -- The nth entry in this list is the absolute value of x_{n} - - -- x_{n+1}. - differences = zipWith abs_diff xn xn_plus_one +-- | Return the number of iterations required to find a fixed point of +-- the function @f@ with the search starting at x0 and tolerance +-- @epsilon@. We delegate to the version that returns the number of +-- iterations and simply discard the fixed point. +fixed_point_iteration_count :: (Vector a, RealFrac b) + => (a -> a) -- ^ The function @f@ to iterate. + -> b -- ^ The tolerance, @epsilon@. + -> a -- ^ The initial value @x0@. + -> Int -- ^ The fixed point. +fixed_point_iteration_count f epsilon x0 = + fst $ F.fixed_point_with_iterations f epsilon x0 - -- A list of pairs, (xn, |x_{n} - x_{n+1}|). - pairs = zip xn differences - -- The pair (xn, |x_{n} - x_{n+1}|) with - -- |x_{n} - x_{n+1}| < epsilon. The pattern match on 'Just' is - -- "safe" since the list is infinite. We'll succeed or loop - -- forever. - Just winning_pair = find (\(_, diff) -> diff < epsilon) pairs +-- | Returns a list of ratios, +-- +-- ||x^{*} - x_{n+1}|| / ||x^{*} - x_{n}||^{p} +-- +-- of fixed point iterations for the function @f@ with initial guess +-- @x0@ and @p@ some positive power. +-- +-- This is used to determine the rate of convergence. +-- +fixed_point_error_ratios :: (Vector a, RealFrac b) + => (a -> a) -- ^ The function @f@ to iterate. + -> a -- ^ The initial value @x0@. + -> a -- ^ The true solution, @x_star@. + -> Integer -- ^ The power @p@. + -> [b] -- ^ The resulting sequence of x_{n}. +fixed_point_error_ratios f x0 x_star p = + zipWith (/) en_plus_one en_exp + where + xn = F.fixed_point_iterations f x0 + en = map (\x -> norm (x_star - x)) xn + en_plus_one = tail en + en_exp = map (^p) en