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

diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs
index 31c7c73..7daf3b6 100644
--- a/src/Linear/Iteration.hs
+++ b/src/Linear/Iteration.hs
@@ -20,6 +20,77 @@ import Linear.Matrix
 import Linear.System
 import 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)
+
+
 -- | Perform one Jacobi iteration,
 --
 --     x1 = M^(-1) * (N*x0 + b)
@@ -40,13 +111,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
@@ -72,8 +138,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 +151,97 @@ 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:
+--
+--   >>> 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:
+--
+--   >>> 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 +261,28 @@ 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:
+--
+--   >>> 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)