Update numeric-prelude and fixed-vector.

[numerical-analysis.git] / src / Roots / Simple.hs

diff --git a/src/Roots/Simple.hs b/src/Roots/Simple.hs
index 0a1debff0a0b1db653efabdf2be00ce9a4b91ed4..a6aa09e497ba841d2c3a8e0be2749aab27c72662 100644 (file)
--- a/src/Roots/Simple.hs
+++ b/src/Roots/Simple.hs
@@ -18,8 +18,9 @@ import Normed
 import qualified Roots.Fast as F
 
 import NumericPrelude hiding (abs)
 import qualified Roots.Fast as F
 
 import NumericPrelude hiding (abs)

-import qualified Algebra.Absolute as Absolute

 import Algebra.Absolute (abs)
 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
 import qualified Algebra.Field as Field
 import qualified Algebra.RealField as RealField
 import qualified Algebra.RealRing as RealRing


@@ -55,7 +56,7 @@ has_root f a b epsilon =
 
 
 -- | We are given a function @f@ and an interval [a,b]. The bisection
 
 
 -- | 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
 --
 --   If one is found within some prescribed tolerance @epsilon@, it is
 --   returned. Otherwise, the interval [a,b] is split into two


@@ -67,8 +68,12 @@ has_root f a b epsilon =
 --
 --   Examples:
 --
 --
 --   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 sin (-1) 1 0.001
 --   Just 0.0


@@ -83,11 +88,45 @@ bisect f a b epsilon =
   F.bisect f a b epsilon Nothing Nothing
 
 
   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.
 --
 -- | 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, RealField.C 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@.
             => (a -> a) -- ^ The function @f@ to iterate.
             -> b       -- ^ The tolerance, @epsilon@.
             -> a       -- ^ The initial value @x0@.


@@ -100,7 +139,10 @@ fixed_point f epsilon x0 =
 --   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.
 --   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 :: (Normed a, RealField.C b)
+fixed_point_iteration_count :: (Normed a,
+                                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@.
                             => (a -> a) -- ^ The function @f@ to iterate.
                             -> b       -- ^ The tolerance, @epsilon@.
                             -> a       -- ^ The initial value @x0@.


@@ -118,7 +160,10 @@ fixed_point_iteration_count f epsilon x0 =
 --
 --   This is used to determine the rate of convergence.
 --
 --
 --   This is used to determine the rate of convergence.
 --

-fixed_point_error_ratios :: (Normed a, RealField.C b)
+fixed_point_error_ratios :: (Normed a,
+                             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@.
                    => (a -> a) -- ^ The function @f@ to iterate.
                    -> a       -- ^ The initial value @x0@.
                    -> a       -- ^ The true solution, @x_star@.