From: Michael Orlitzky Date: Sun, 21 Jul 2013 21:21:12 +0000 (-0400) Subject: Add Gauss-Seidel and SOR iterative methods. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=eb2fb001d7b7a7769844b6e5c3623e7501a7ac12 Add Gauss-Seidel and SOR iterative methods. --- diff --git a/src/Linear/Iteration.hs b/src/Linear/Iteration.hs index 31c7c73..89f2f0c 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