Add the rayleigh_quotient function.

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

diff --git a/src/Roots/Fast.hs b/src/Roots/Fast.hs
index 0deb1fd6237a5909ea7d15a252f09b524093d1bf..8b69786379218b12448b114d4183ca9c3ef64c5d 100644 (file)
--- a/src/Roots/Fast.hs
+++ b/src/Roots/Fast.hs
@@ -9,6 +9,7 @@ module Roots.Fast
 where
 
 import Data.List (find)
 where
 
 import Data.List (find)

+import Data.Maybe (fromMaybe)

 
 import Normed
 
 
 import Normed
 


@@ -16,7 +17,6 @@ import NumericPrelude hiding (abs)
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Algebraic as Algebraic
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Algebraic as Algebraic

-import qualified Algebra.Field as Field

 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.RealField as RealField
 
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.RealField as RealField
 


@@ -31,33 +31,20 @@ has_root :: (RealField.C a,
          -> Maybe b -- ^ Precoumpted f(a)
          -> Maybe b -- ^ Precoumpted f(b)
          -> Bool
          -> Maybe b -- ^ Precoumpted f(a)
          -> Maybe b -- ^ Precoumpted f(b)
          -> Bool

-has_root f a b epsilon f_of_a f_of_b =
-  if not ((signum (f_of_a')) * (signum (f_of_b')) == 1) then
-    -- We don't care about epsilon here, there's definitely a root!
-    True
-  else
-    if (b - a) <= epsilon' then
-      -- Give up, return false.
-      False
-    else
-      -- If either [a,c] or [c,b] have roots, we do too.
+has_root f a b epsilon f_of_a f_of_b
+  | (signum (f_of_a')) * (signum (f_of_b')) /= 1 = True
+  | (b - a) <= epsilon' = False
+  | otherwise =

       (has_root f a c (Just epsilon') (Just f_of_a') Nothing) ||
         (has_root f c b (Just epsilon') Nothing (Just f_of_b'))
   where
     -- If the size of the smallest subinterval is not specified,
     -- assume we just want to check once on all of [a,b].
       (has_root f a c (Just epsilon') (Just f_of_a') Nothing) ||
         (has_root f c b (Just epsilon') Nothing (Just f_of_b'))
   where
     -- If the size of the smallest subinterval is not specified,
     -- assume we just want to check once on all of [a,b].

-    epsilon' = case epsilon of
-                 Nothing -> (b-a)
-                 Just eps -> eps
+    epsilon' = fromMaybe (b-a) epsilon

 
     -- Compute f(a) and f(b) only if needed.
 
     -- Compute f(a) and f(b) only if needed.

-    f_of_a'  = case f_of_a of
-                 Nothing -> f a
-                 Just v -> v
-
-    f_of_b'  = case f_of_b of
-                 Nothing -> f b
-                 Just v -> v
+    f_of_a'  = fromMaybe (f a) f_of_a
+    f_of_b'  = fromMaybe (f b) f_of_b

 
     c = (a + b)/2
 
 
     c = (a + b)/2
 


@@ -88,26 +75,64 @@ bisect f a b epsilon f_of_a f_of_b
         else bisect f c b epsilon (Just f_of_c') (Just f_of_b')
   where
     -- Compute f(a) and f(b) only if needed.
         else bisect f c b epsilon (Just f_of_c') (Just f_of_b')
   where
     -- Compute f(a) and f(b) only if needed.

-    f_of_a'  = case f_of_a of
-                 Nothing -> f a
-                 Just v -> v
-
-    f_of_b'  = case f_of_b of
-                 Nothing -> f b
-                 Just v -> v
+    f_of_a'  = fromMaybe (f a) f_of_a
+    f_of_b'  = fromMaybe (f b) f_of_b

 
     c = (a + b) / 2
 
 
 
 
     c = (a + b) / 2
 
 
 

+trisect :: (RealField.C a,
+            RealRing.C b,
+            Absolute.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 b -- ^ Precomputed f(a)
+        -> Maybe b -- ^ Precomputed f(b)
+        -> Maybe a
+trisect f a b epsilon f_of_a f_of_b
+  -- We pass @epsilon@ to the 'has_root' function because if we want a
+  -- result within epsilon of the true root, we need to know that
+  -- there *is* a root within an interval of length epsilon.
+  | not (has_root f a b (Just epsilon) (Just f_of_a') (Just f_of_b')) = Nothing
+  | f_of_a' == 0 = Just a
+  | f_of_b' == 0 = Just b
+  | otherwise =
+      -- Use a 'prime' just for consistency.
+    let (a', b', fa', fb')
+          | has_root f d b (Just epsilon) (Just f_of_d') (Just f_of_b') =
+              (d, b, f_of_d', f_of_b')
+          | has_root f c d (Just epsilon) (Just f_of_c') (Just f_of_d') =
+              (c, d, f_of_c', f_of_d')
+          | otherwise =
+              (a, c, f_of_a', f_of_c')
+    in
+      if (b-a) < 2*epsilon
+      then Just ((b+a)/2)
+      else trisect f a' b' epsilon (Just fa') (Just fb')
+  where
+    -- Compute f(a) and f(b) only if needed.
+    f_of_a'  = fromMaybe (f a) f_of_a
+    f_of_b'  = fromMaybe (f b) f_of_b
+
+    c = (2*a + b) / 3
+
+    d = (a + 2*b) / 3
+
+    f_of_c' = f c
+    f_of_d' = f d
+
+

 
 -- | 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}.
 
 -- | 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
+fixed_point_iterations =
+  iterate

 
 
 -- | Find a fixed point of the function @f@ with the search starting
 
 
 -- | Find a fixed point of the function @f@ with the search starting


@@ -118,7 +143,7 @@ fixed_point_iterations f x0 =
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,
 --   We also return the number of iterations required.
 --
 fixed_point_with_iterations :: (Normed a,

-                                Algebraic.C a,
+                                Additive.C a,

                                 RealField.C b,
                                 Algebraic.C b)
                             => (a -> a)  -- ^ The function @f@ to iterate.
                                 RealField.C b,
                                 Algebraic.C b)
                             => (a -> a)  -- ^ The function @f@ to iterate.