X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FIteration.hs;h=06d30a802aa55717b6a9ec45b42c69a8b8f8b473;hb=ae914d13235a4582077a5cb2b1edd630d9c6ad62;hp=31c7c738edbfb0500d4689af4875ba5fb08ddc48;hpb=2e2df5a19f1ef4fbc1a5194d0970d3886c8d5fd5;p=numerical-analysis.git

diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs
index 31c7c73..06d30a8 100644
--- a/src/Linear/Iteration.hs
+++ b/src/Linear/Iteration.hs
@@ -1,24 +1,113 @@
+{-# LANGUAGE NoImplicitPrelude #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
 
 -- | Classical iterative methods to solve the system Ax = b.
 
-module Linear.Iteration
+module Linear.Iteration (
+  gauss_seidel_iteration,
+  gauss_seidel_iterations,
+  gauss_seidel_method,
+  jacobi_iteration,
+  jacobi_iterations,
+  jacobi_method,
+  rayleigh_quotient,
+  sor_iteration,
+  sor_iterations,
+  sor_method )
 where
 
-import Data.List (find)
-import Data.Maybe (fromJust)
-import Data.Vector.Fixed (Arity, N1, S)
-import NumericPrelude hiding ((*))
-import qualified Algebra.Algebraic as Algebraic
-import qualified Algebra.Field as Field
-import qualified Algebra.RealField as RealField
-import qualified Algebra.ToRational as ToRational
-import qualified Prelude as P
+import Data.List ( find )
+import Data.Maybe ( fromJust )
+import Data.Vector.Fixed ( Arity, N1, S )
+import NumericPrelude hiding ( (*) )
+import qualified Algebra.Algebraic as Algebraic ( C )
+import qualified Algebra.Field as Field ( C )
+import qualified Algebra.RealField as RealField ( C )
+import qualified Algebra.ToRational as ToRational ( C )
+
+import Linear.Matrix (
+  Mat(..),
+  (!!!),
+  (*),
+  diagonal_part,
+  dot,
+  lt_part_strict,
+  transpose )
+import Linear.System ( forward_substitute )
+import Normed ( Normed(..) )
+
+
+-- | A generalized implementation for Jacobi, Gauss-Seidel, etc. All
+--   that we really need to know is how to construct the matrix M, so we
+--   take a function that does it as an argument.
+classical_iteration :: (Field.C a, Arity m)
+                 => (Mat m m  a -> Mat m m  a)
+                 -> Mat m m  a
+                 -> Mat m N1 a
+                 -> Mat m N1 a
+                 -> Mat m N1 a
+classical_iteration m_function matrix b x_current =
+  x_next
+  where
+    big_m = m_function matrix
+    big_n = big_m - matrix
+    rhs = big_n*x_current + b
+    -- TODO: Should be solve below! M might not be lower-triangular.
+    x_next = forward_substitute big_m rhs
+
+
+-- | Perform one iteration of successive over-relaxation.
+--
+sor_iteration :: forall m a.
+                 (Field.C a, Arity m)
+              => a -- ^ Omega
+              -> Mat m m  a -- ^ Matrix A
+              -> Mat m N1 a -- ^ Vector b
+              -> Mat m N1 a -- ^ Vector x_current
+              -> Mat m N1 a -- ^ Output vector x_next
+sor_iteration omega =
+  classical_iteration m_function
+  where
+    m_function :: Mat m m a -> Mat m m a
+    m_function matrix =
+      let diag = (recip omega) *> (diagonal_part matrix)
+          lt   = lt_part_strict matrix
+      in
+        diag + lt
+
+
+-- | Compute an infinite list of SOR iterations starting with the
+--   vector x0.
+sor_iterations :: (Field.C a, Arity m)
+               => a
+               -> Mat m m  a
+               -> Mat m N1 a
+               -> Mat m N1 a
+               -> [Mat m N1 a]
+sor_iterations omega matrix b =
+  iterate (sor_iteration omega matrix b)
+
+
+-- | Perform one iteration of Gauss-Seidel.
+gauss_seidel_iteration :: (Field.C a, Arity m)
+                       => Mat m m  a
+                       -> Mat m N1 a
+                       -> Mat m N1 a
+                       -> Mat m N1 a
+gauss_seidel_iteration = sor_iteration one
+
+
+-- | Compute an infinite list of Gauss-Seidel iterations starting with
+--   the vector x0.
+gauss_seidel_iterations :: (Field.C a, Arity m)
+                        => Mat m m  a
+                        -> Mat m N1 a
+                        -> Mat m N1 a
+                        -> [Mat m N1 a]
+gauss_seidel_iterations matrix b =
+  iterate (gauss_seidel_iteration matrix b)
 
-import Linear.Matrix
-import Linear.System
-import Normed
 
 -- | Perform one Jacobi iteration,
 --
@@ -26,6 +115,8 @@ import Normed
 --
 --   Examples:
 --
+--   >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 --   >>> let x0 = vec2d (0, 0::Double)
 --   >>> let b  = vec2d (1, 1::Double)
@@ -40,13 +131,8 @@ jacobi_iteration :: (Field.C a, Arity m)
                  -> Mat m N1 a
                  -> Mat m N1 a
                  -> Mat m N1 a
-jacobi_iteration matrix b x_current =
-  x_next
-  where
-    big_m = diagonal matrix
-    big_n = big_m - matrix
-    rhs = big_n*x_current + b
-    x_next = forward_substitute big_m rhs
+jacobi_iteration =
+  classical_iteration diagonal_part
 
 
 -- | Compute an infinite list of Jacobi iterations starting with the
@@ -65,6 +151,8 @@ jacobi_iterations matrix b =
 --
 --   Examples:
 --
+--   >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
 --   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
 --   >>> let x0 = vec2d (0, 0::Double)
 --   >>> let b  = vec2d (1, 1::Double)
@@ -72,8 +160,7 @@ jacobi_iterations matrix b =
 --   >>> jacobi_method m b x0 epsilon
 --   ((0.0),(0.4999995231628418))
 --
-jacobi_method :: forall m n a b.
-                 (RealField.C a,
+jacobi_method :: (RealField.C a,
                   Algebraic.C a, -- Normed instance
                   ToRational.C a, -- Normed instance
                   Algebraic.C b,
@@ -86,13 +173,101 @@ jacobi_method :: forall m n a b.
               -> Mat m N1 a
               -> b
               -> Mat m N1 a
-jacobi_method matrix b x0 epsilon =
+jacobi_method =
+  classical_method jacobi_iterations
+
+
+-- | Solve the system Ax = b using the Gauss-Seidel method. This will
+--   run forever if the iterations do not converge.
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
+--   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
+--   >>> let x0 = vec2d (0, 0::Double)
+--   >>> let b  = vec2d (1, 1::Double)
+--   >>> let epsilon = 10**(-12)
+--   >>> gauss_seidel_method m b x0 epsilon
+--   ((4.547473508864641e-13),(0.49999999999954525))
+--
+gauss_seidel_method :: (RealField.C a,
+                        Algebraic.C a, -- Normed instance
+                        ToRational.C a, -- Normed instance
+                        Algebraic.C b,
+                        RealField.C b,
+                        Arity m,
+                        Arity n, -- Normed instance
+                        m ~ S n)
+                    => Mat m m  a
+                    -> Mat m N1 a
+                    -> Mat m N1 a
+                    -> b
+                    -> Mat m N1 a
+gauss_seidel_method =
+  classical_method gauss_seidel_iterations
+
+
+-- | Solve the system Ax = b using the Successive Over-Relaxation
+--   (SOR) method. This will run forever if the iterations do not
+--   converge.
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
+--   >>> let m  = fromList [[4,2],[2,2]] :: Mat2 Double
+--   >>> let x0 = vec2d (0, 0::Double)
+--   >>> let b  = vec2d (1, 1::Double)
+--   >>> let epsilon = 10**(-12)
+--   >>> sor_method 1.5 m b x0 epsilon
+--   ((6.567246746413957e-13),(0.4999999999993727))
+--
+sor_method :: (RealField.C a,
+               Algebraic.C a, -- Normed instance
+               ToRational.C a, -- Normed instance
+               Algebraic.C b,
+               RealField.C b,
+               Arity m,
+               Arity n, -- Normed instance
+               m ~ S n)
+           => a
+           -> Mat m m  a
+           -> Mat m N1 a
+           -> Mat m N1 a
+           -> b
+           -> Mat m N1 a
+sor_method omega =
+  classical_method (sor_iterations omega)
+
+
+-- | General implementation for all classical iteration methods. For
+--   its first argument, it takes a function which generates the
+--   sequence of iterates when supplied with the remaining arguments
+--   (except for the tolerance).
+--
+classical_method :: forall m n a b.
+                    (RealField.C a,
+                     Algebraic.C a, -- Normed instance
+                     ToRational.C a, -- Normed instance
+                     Algebraic.C b,
+                     RealField.C b,
+                     Arity m,
+                     Arity n, -- Normed instance
+                     m ~ S n)
+                 => (Mat m m  a -> Mat m N1 a -> Mat m N1 a -> [Mat m N1 a])
+                 -> Mat m m  a
+                 -> Mat m N1 a
+                 -> Mat m N1 a
+                 -> b
+                 -> Mat m N1 a
+classical_method iterations_function matrix b x0 epsilon =
   -- fromJust is "safe," because the list is infinite. If the
   -- algorithm doesn't converge, 'find' will search forever and never
   -- return Nothing.
   fst' $ fromJust $ find error_small_enough diff_pairs
   where
-    x_n = jacobi_iterations matrix b x0
+    x_n = iterations_function matrix b x0
 
     pairs :: [(Mat m N1 a, Mat m N1 a)]
     pairs = zip (tail x_n) x_n
@@ -112,3 +287,30 @@ jacobi_method matrix b x0 epsilon =
 
     error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
     error_small_enough (_,_,err) = err < epsilon
+
+
+
+-- | Compute the Rayleigh quotient of  @matrix@ and @vector@.
+--
+--   Examples:
+--
+--   >>> import Linear.Matrix ( Mat2, fromList, vec2d )
+--
+--   >>> let m = fromList [[3,1],[1,2]] :: Mat2 Rational
+--   >>> let v = vec2d (1, 1::Rational)
+--   >>> rayleigh_quotient m v
+--   7 % 2
+--
+rayleigh_quotient :: (RealField.C a,
+                      Arity m,
+                      Arity n,
+                      m ~ S n)
+                  => (Mat m m a)
+                  -> (Mat m N1 a)
+                  -> a
+rayleigh_quotient matrix vector =
+  (vector `dot` (matrix * vector)) / (norm_squared vector)
+  where
+    -- We don't use the norm function here to avoid the algebraic
+    -- requirement on our field.
+    norm_squared v = ((transpose v) * v) !!! (0,0)