]> gitweb.michael.orlitzky.com - numerical-analysis.git/blob - src/Linear/Iteration.hs
89f2f0c0eab7788a088642afd81e98525478e610
[numerical-analysis.git] / src / Linear / Iteration.hs
1 {-# LANGUAGE ScopedTypeVariables #-}
2 {-# LANGUAGE TypeFamilies #-}
3
4 -- | Classical iterative methods to solve the system Ax = b.
5
6 module Linear.Iteration
7 where
8
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
18
19 import Linear.Matrix
20 import Linear.System
21 import Normed
22
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)
28 -> Mat m m a
29 -> Mat m N1 a
30 -> Mat m N1 a
31 -> Mat m N1 a
32 classical_iteration m_function matrix b x_current =
33 x_next
34 where
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
40
41
42 -- | Perform one iteration of successive over-relaxation.
43 --
44 sor_iteration :: forall m a.
45 (Field.C a, Arity m)
46 => a -- ^ Omega
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
51 sor_iteration omega =
52 classical_iteration m_function
53 where
54 m_function :: Mat m m a -> Mat m m a
55 m_function matrix =
56 let diag = (recip omega) *> (diagonal_part matrix)
57 lt = lt_part_strict matrix
58 in
59 diag + lt
60
61
62 -- | Compute an infinite list of SOR iterations starting with the
63 -- vector x0.
64 sor_iterations :: (Field.C a, Arity m)
65 => a
66 -> Mat m m a
67 -> Mat m N1 a
68 -> Mat m N1 a
69 -> [Mat m N1 a]
70 sor_iterations omega matrix b =
71 iterate (sor_iteration omega matrix b)
72
73
74 -- | Perform one iteration of Gauss-Seidel.
75 gauss_seidel_iteration :: (Field.C a, Arity m)
76 => Mat m m a
77 -> Mat m N1 a
78 -> Mat m N1 a
79 -> Mat m N1 a
80 gauss_seidel_iteration = sor_iteration one
81
82
83 -- | Compute an infinite list of Gauss-Seidel iterations starting with
84 -- the vector x0.
85 gauss_seidel_iterations :: (Field.C a, Arity m)
86 => Mat m m a
87 -> Mat m N1 a
88 -> Mat m N1 a
89 -> [Mat m N1 a]
90 gauss_seidel_iterations matrix b =
91 iterate (gauss_seidel_iteration matrix b)
92
93
94 -- | Perform one Jacobi iteration,
95 --
96 -- x1 = M^(-1) * (N*x0 + b)
97 --
98 -- Examples:
99 --
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
104 -- ((0.25),(0.5))
105 -- >>> let x1 = jacobi_iteration m b x0
106 -- >>> jacobi_iteration m b x1
107 -- ((0.0),(0.25))
108 --
109 jacobi_iteration :: (Field.C a, Arity m)
110 => Mat m m a
111 -> Mat m N1 a
112 -> Mat m N1 a
113 -> Mat m N1 a
114 jacobi_iteration =
115 classical_iteration diagonal_part
116
117
118 -- | Compute an infinite list of Jacobi iterations starting with the
119 -- vector x0.
120 jacobi_iterations :: (Field.C a, Arity m)
121 => Mat m m a
122 -> Mat m N1 a
123 -> Mat m N1 a
124 -> [Mat m N1 a]
125 jacobi_iterations matrix b =
126 iterate (jacobi_iteration matrix b)
127
128
129 -- | Solve the system Ax = b using the Jacobi method. This will run
130 -- forever if the iterations do not converge.
131 --
132 -- Examples:
133 --
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))
140 --
141 jacobi_method :: (RealField.C a,
142 Algebraic.C a, -- Normed instance
143 ToRational.C a, -- Normed instance
144 Algebraic.C b,
145 RealField.C b,
146 Arity m,
147 Arity n, -- Normed instance
148 m ~ S n)
149 => Mat m m a
150 -> Mat m N1 a
151 -> Mat m N1 a
152 -> b
153 -> Mat m N1 a
154 jacobi_method =
155 classical_method jacobi_iterations
156
157
158 -- | Solve the system Ax = b using the Gauss-Seidel method. This will
159 -- run forever if the iterations do not converge.
160 --
161 -- Examples:
162 --
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))
169 --
170 gauss_seidel_method :: (RealField.C a,
171 Algebraic.C a, -- Normed instance
172 ToRational.C a, -- Normed instance
173 Algebraic.C b,
174 RealField.C b,
175 Arity m,
176 Arity n, -- Normed instance
177 m ~ S n)
178 => Mat m m a
179 -> Mat m N1 a
180 -> Mat m N1 a
181 -> b
182 -> Mat m N1 a
183 gauss_seidel_method =
184 classical_method gauss_seidel_iterations
185
186
187 -- | Solve the system Ax = b using the Successive Over-Relaxation
188 -- (SOR) method. This will run forever if the iterations do not
189 -- converge.
190 --
191 -- Examples:
192 --
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))
199 --
200 sor_method :: (RealField.C a,
201 Algebraic.C a, -- Normed instance
202 ToRational.C a, -- Normed instance
203 Algebraic.C b,
204 RealField.C b,
205 Arity m,
206 Arity n, -- Normed instance
207 m ~ S n)
208 => a
209 -> Mat m m a
210 -> Mat m N1 a
211 -> Mat m N1 a
212 -> b
213 -> Mat m N1 a
214 sor_method omega =
215 classical_method (sor_iterations omega)
216
217
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).
222 --
223 classical_method :: forall m n a b.
224 (RealField.C a,
225 Algebraic.C a, -- Normed instance
226 ToRational.C a, -- Normed instance
227 Algebraic.C b,
228 RealField.C b,
229 Arity m,
230 Arity n, -- Normed instance
231 m ~ S n)
232 => (Mat m m a -> Mat m N1 a -> Mat m N1 a -> [Mat m N1 a])
233 -> Mat m m a
234 -> Mat m N1 a
235 -> Mat m N1 a
236 -> b
237 -> 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
241 -- return Nothing.
242 fst' $ fromJust $ find error_small_enough diff_pairs
243 where
244 x_n = iterations_function matrix b x0
245
246 pairs :: [(Mat m N1 a, Mat m N1 a)]
247 pairs = zip (tail x_n) x_n
248
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) =
252 (cur,prev,diff)
253 where
254 diff = norm (cur - prev)
255
256 diff_pairs :: [(Mat m N1 a, Mat m N1 a, b)]
257 diff_pairs = map append_diff pairs
258
259 fst' :: (c, d, e) -> c
260 fst' (x,_,_) = x
261
262 error_small_enough :: (Mat m N1 a, Mat m N1 a, b)-> Bool
263 error_small_enough (_,_,err) = err < epsilon