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)
-> 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
-- >>> 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,
-> 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
projects / numerical-analysis.git / commitdiff