X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FRoots%2FSimple.hs;h=a6aa09e497ba841d2c3a8e0be2749aab27c72662;hb=04f56a8882bb0c574b603f8c3fed9481ea934f7f;hp=44d3d62112d2c4f339fa16b65c143e65ed5bb83a;hpb=fe73028041fe3becce6ce1ff268181d55d54a011;p=numerical-analysis.git diff --git a/src/Roots/Simple.hs b/src/Roots/Simple.hs index 44d3d62..a6aa09e 100644 --- a/src/Roots/Simple.hs +++ b/src/Roots/Simple.hs @@ -18,9 +18,12 @@ import Normed import qualified Roots.Fast as F import NumericPrelude hiding (abs) -import Algebra.Absolute -import Algebra.Field -import Algebra.Ring +import Algebra.Absolute (abs) +import qualified Algebra.Additive as Additive +import qualified Algebra.Algebraic as Algebraic +import qualified Algebra.Field as Field +import qualified Algebra.RealField as RealField +import qualified Algebra.RealRing as RealRing -- | Does the (continuous) function @f@ have a root on the interval -- [a,b]? If f(a) <] 0 and f(b) ]> 0, we know that there's a root in @@ -41,11 +44,7 @@ import Algebra.Ring -- >>> has_root cos (-2) 2 (Just 0.001) -- True -- -has_root :: (Algebra.Field.C a, - Ord a, - Algebra.Ring.C b, - Algebra.Absolute.C b, - Ord b) +has_root :: (RealField.C a, RealRing.C b) => (a -> b) -- ^ The function @f@ -> a -- ^ The \"left\" endpoint, @a@ -> a -- ^ The \"right\" endpoint, @b@ @@ -57,7 +56,7 @@ has_root f a b epsilon = -- | We are given a function @f@ and an interval [a,b]. The bisection --- method checks finds a root by splitting [a,b] in half repeatedly. +-- method finds a root by splitting [a,b] in half repeatedly. -- -- If one is found within some prescribed tolerance @epsilon@, it is -- returned. Otherwise, the interval [a,b] is split into two @@ -69,17 +68,17 @@ has_root f a b epsilon = -- -- Examples: -- --- >>> bisect cos 1 2 0.001 --- Just 1.5712890625 +-- >>> let actual = 1.5707963267948966 +-- >>> let Just root = bisect cos 1 2 0.001 +-- >>> root +-- 1.5712890625 +-- >>> abs (root - actual) < 0.001 +-- True -- -- >>> bisect sin (-1) 1 0.001 -- Just 0.0 -- -bisect :: (Algebra.Field.C a, - Ord a, - Algebra.Ring.C b, - Algebra.Absolute.C b, - Ord b) +bisect :: (RealField.C a, RealRing.C b) => (a -> b) -- ^ The function @f@ whose root we seek -> a -- ^ The \"left\" endpoint of the interval, @a@ -> a -- ^ The \"right\" endpoint of the interval, @b@ @@ -89,14 +88,45 @@ bisect f a b epsilon = F.bisect f a b epsilon Nothing Nothing +-- | We are given a function @f@ and an interval [a,b]. The trisection +-- method finds a root by splitting [a,b] into three +-- subintervals repeatedly. +-- +-- If one is found within some prescribed tolerance @epsilon@, it is +-- returned. Otherwise, the interval [a,b] is split into two +-- subintervals [a,c] and [c,b] of equal length which are then both +-- checked via the same process. +-- +-- Returns 'Just' the value x for which f(x) == 0 if one is found, +-- or Nothing if one of the preconditions is violated. +-- +-- Examples: +-- +-- >>> let actual = 1.5707963267948966 +-- >>> let Just root = trisect cos 1 2 0.001 +-- >>> root +-- 1.5713305898491083 +-- >>> abs (root - actual) < 0.001 +-- True +-- +-- >>> trisect sin (-1) 1 0.001 +-- Just 0.0 +-- +trisect :: (RealField.C a, RealRing.C b) + => (a -> b) -- ^ The function @f@ whose root we seek + -> a -- ^ The \"left\" endpoint of the interval, @a@ + -> a -- ^ The \"right\" endpoint of the interval, @b@ + -> a -- ^ The tolerance, @epsilon@ + -> Maybe a +trisect f a b epsilon = + F.trisect f a b epsilon Nothing Nothing + + -- | Find a fixed point of the function @f@ with the search starting -- at x0. We delegate to the version that returns the number of -- iterations and simply discard the number of iterations. -- -fixed_point :: (Normed a, - Algebra.Field.C b, - Algebra.Absolute.C b, - Ord b) +fixed_point :: (Normed a, Additive.C a, Algebraic.C b, RealField.C b) => (a -> a) -- ^ The function @f@ to iterate. -> b -- ^ The tolerance, @epsilon@. -> a -- ^ The initial value @x0@. @@ -110,9 +140,9 @@ fixed_point f epsilon x0 = -- @epsilon@. We delegate to the version that returns the number of -- iterations and simply discard the fixed point. fixed_point_iteration_count :: (Normed a, - Algebra.Field.C b, - Algebra.Absolute.C b, - Ord b) + Additive.C a, + RealField.C b, + Algebraic.C b) => (a -> a) -- ^ The function @f@ to iterate. -> b -- ^ The tolerance, @epsilon@. -> a -- ^ The initial value @x0@. @@ -131,9 +161,9 @@ fixed_point_iteration_count f epsilon x0 = -- This is used to determine the rate of convergence. -- fixed_point_error_ratios :: (Normed a, - Algebra.Field.C b, - Algebra.Absolute.C b, - Ord b) + Additive.C a, + RealField.C b, + Algebraic.C b) => (a -> a) -- ^ The function @f@ to iterate. -> a -- ^ The initial value @x0@. -> a -- ^ The true solution, @x_star@. @@ -160,7 +190,7 @@ fixed_point_error_ratios f x0 x_star p = -- >>> tail $ take 4 $ newton_iterations f f' 2 -- [1.6806282722513088,1.4307389882390624,1.2549709561094362] -- -newton_iterations :: (Algebra.Field.C a) +newton_iterations :: (Field.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> (a -> a) -- ^ The derivative of @f@ -> a -- ^ Initial guess, x-naught @@ -195,7 +225,7 @@ newton_iterations f f' x0 = -- >>> abs (f root) < eps -- True -- -newtons_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a) +newtons_method :: (RealField.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> (a -> a) -- ^ The derivative of @f@ -> a -- ^ The tolerance epsilon @@ -245,7 +275,7 @@ iterate2 f x0 x1 = -- >>> take 4 $ secant_iterations f 2 1 -- [2.0,1.0,1.0161290322580645,1.190577768676638] -- -secant_iterations :: (Algebra.Field.C a) +secant_iterations :: (Field.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> a -- ^ Initial guess, x-naught -> a -- ^ Second initial guess, x-one @@ -273,7 +303,7 @@ secant_iterations f x0 x1 = -- >>> abs (f root) < (1/10^9) -- True -- -secant_method :: (Algebra.Field.C a, Algebra.Absolute.C a, Ord a) +secant_method :: (RealField.C a) => (a -> a) -- ^ The function @f@ whose root we seek -> a -- ^ The tolerance epsilon -> a -- ^ Initial guess, x-naught