1 {-# LANGUAGE ScopedTypeVariables #-}
2 {-# LANGUAGE TypeFamilies #-}
4 -- | Classical iterative methods to solve the system Ax = b.
6 module Linear.Iteration
9 import Data.List (find)
10 import Data.Maybe (fromJust)
11 import Data.Vector.Fixed (Arity, N1, S)
12 import NumericPrelude hiding ((*))
13 import qualified Algebra.Algebraic as Algebraic
14 import qualified Algebra.Field as Field
15 import qualified Algebra.RealField as RealField
16 import qualified Algebra.ToRational as ToRational
17 import qualified Prelude as P
23 -- | A generalized implementation for Jacobi, Gauss-Seidel, etc. All
24 -- that we really need to know is how to construct the matrix M, so we
25 -- take a function that does it as an argument.
26 classical_iteration :: (Field.C a, Arity m)
27 => (Mat m m a -> Mat m m a)
32 classical_iteration m_function matrix b x_current =
35 big_m = m_function matrix
36 big_n = big_m - matrix
37 rhs = big_n*x_current + b
38 -- TODO: Should be solve below! M might not be lower-triangular.
39 x_next = forward_substitute big_m rhs
42 -- | Perform one iteration of successive over-relaxation.
44 sor_iteration :: forall m a.
47 -> Mat m m a -- ^ Matrix A
48 -> Mat m N1 a -- ^ Vector b
49 -> Mat m N1 a -- ^ Vector x_current
50 -> Mat m N1 a -- ^ Output vector x_next
52 classical_iteration m_function
54 m_function :: Mat m m a -> Mat m m a
56 let diag = (recip omega) *> (diagonal_part matrix)
57 lt = lt_part_strict matrix
62 -- | Compute an infinite list of SOR iterations starting with the
64 sor_iterations :: (Field.C a, Arity m)
70 sor_iterations omega matrix b =
71 iterate (sor_iteration omega matrix b)
74 -- | Perform one iteration of Gauss-Seidel.
75 gauss_seidel_iteration :: (Field.C a, Arity m)
80 gauss_seidel_iteration = sor_iteration one
83 -- | Compute an infinite list of Gauss-Seidel iterations starting with
85 gauss_seidel_iterations :: (Field.C a, Arity m)
90 gauss_seidel_iterations matrix b =
91 iterate (gauss_seidel_iteration matrix b)
94 -- | Perform one Jacobi iteration,
96 -- x1 = M^(-1) * (N*x0 + b)
100 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
101 -- >>> let x0 = vec2d (0, 0::Double)
102 -- >>> let b = vec2d (1, 1::Double)
103 -- >>> jacobi_iteration m b x0
105 -- >>> let x1 = jacobi_iteration m b x0
106 -- >>> jacobi_iteration m b x1
109 jacobi_iteration :: (Field.C a, Arity m)
115 classical_iteration diagonal_part
118 -- | Compute an infinite list of Jacobi iterations starting with the
120 jacobi_iterations :: (Field.C a, Arity m)
125 jacobi_iterations matrix b =
126 iterate (jacobi_iteration matrix b)
129 -- | Solve the system Ax = b using the Jacobi method. This will run
130 -- forever if the iterations do not converge.
134 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
135 -- >>> let x0 = vec2d (0, 0::Double)
136 -- >>> let b = vec2d (1, 1::Double)
137 -- >>> let epsilon = 10**(-6)
138 -- >>> jacobi_method m b x0 epsilon
139 -- ((0.0),(0.4999995231628418))
141 jacobi_method :: (RealField.C a,
142 Algebraic.C a, -- Normed instance
143 ToRational.C a, -- Normed instance
147 Arity n, -- Normed instance
155 classical_method jacobi_iterations
158 -- | Solve the system Ax = b using the Gauss-Seidel method. This will
159 -- run forever if the iterations do not converge.
163 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
164 -- >>> let x0 = vec2d (0, 0::Double)
165 -- >>> let b = vec2d (1, 1::Double)
166 -- >>> let epsilon = 10**(-12)
167 -- >>> gauss_seidel_method m b x0 epsilon
168 -- ((4.547473508864641e-13),(0.49999999999954525))
170 gauss_seidel_method :: (RealField.C a,
171 Algebraic.C a, -- Normed instance
172 ToRational.C a, -- Normed instance
176 Arity n, -- Normed instance
183 gauss_seidel_method =
184 classical_method gauss_seidel_iterations
187 -- | Solve the system Ax = b using the Successive Over-Relaxation
188 -- (SOR) method. This will run forever if the iterations do not
193 -- >>> let m = fromList [[4,2],[2,2]] :: Mat2 Double
194 -- >>> let x0 = vec2d (0, 0::Double)
195 -- >>> let b = vec2d (1, 1::Double)
196 -- >>> let epsilon = 10**(-12)
197 -- >>> sor_method 1.5 m b x0 epsilon
198 -- ((6.567246746413957e-13),(0.4999999999993727))
200 sor_method :: (RealField.C a,
201 Algebraic.C a, -- Normed instance
202 ToRational.C a, -- Normed instance
206 Arity n, -- Normed instance
215 classical_method (sor_iterations omega)
218 -- | General implementation for all classical iteration methods. For
219 -- its first argument, it takes a function which generates the
220 -- sequence of iterates when supplied with the remaining arguments
221 -- (except for the tolerance).
223 classical_method :: forall m n a b.
225 Algebraic.C a, -- Normed instance
226 ToRational.C a, -- Normed instance
230 Arity n, -- Normed instance
232 => (Mat m m a -> Mat m N1 a -> Mat m N1 a -> [Mat m N1 a])
238 classical_method iterations_function matrix b x0 epsilon =
239 -- fromJust is "safe," because the list is infinite. If the
240 -- algorithm doesn't converge, 'find' will search forever and never
242 fst' $ fromJust $ find error_small_enough diff_pairs
244 x_n = iterations_function matrix b x0
246 pairs :: [(Mat m N1 a, Mat m N1 a)]
247 pairs = zip (tail x_n) x_n
249 append_diff :: (Mat m N1 a, Mat m N1 a)
250 -> (Mat m N1 a, Mat m N1 a, b)
251 append_diff (cur,prev) =
254 diff = norm (cur - prev)
256 diff_pairs :: [(Mat m N1 a, Mat m N1 a, b)]
257 diff_pairs = map append_diff pairs
259 fst' :: (c, d, e) -> c
262 error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
263 error_small_enough (_,_,err) = err < epsilon