1 function [x, iterations, residual_norms] = ...
2 successive_over_relaxation(omega, A, b, x0, ...
3 tolerance, max_iterations)
4 %
5 % Solve the system,
6 %
7 % Ax = b
8 %
9 % iteratively using SOR. That is, we let,
10 %
11 % A = M - N = ((1/omega)*D + L) - U
12 %
13 % where D is a diagonal matrix consisting of the diagonal entries of
14 % A (the rest zeros), and L,U are the upper- and lower-triangular
15 % parts of A respectively. Now,
16 %
17 % Ax = (M - N)x = Mx - Nx = b
18 %
19 % has solution,
20 %
21 % x = M^(-1)*(Nx + b)
22 %
23 % Thus, our iterations are of the form,
24 %
25 % x_{k+1} = M^(-1)*(Nx_{k} + b)
26 %
27 % INPUT:
28 %
29 % ``omega`` -- The relaxation factor.
30 %
31 % ``A`` -- The n-by-n coefficient matrix of the system.
32 %
33 % ``b`` -- An n-by-1 vector; the right-hand side of the system.
34 %
35 % ``x0`` -- An n-by-1 vector; an initial guess to the solution.
36 %
37 % ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance.
38 % we stop when the relative error (the infinity norm
39 % of the residual divided by the infinity norm of
40 % ``b``) is less than ``tolerance``.
41 %
42 % ``max_iterations`` -- (optional; default: intmax) the maximum
43 % number of iterations we will perform.
44 %
45 % OUTPUT:
46 %
47 % ``x`` -- An n-by-1 vector; the approximate solution to the system.
48 %
49 % ``iterations`` -- The number of iterations taken.
50 %
51 % ``residual_norms`` -- An n-by-iterations vector of the residual
52 % (infinity)norms at each iteration. If not
53 % requested, they will not be computed to
54 % save space.
55 %
57 if (nargin < 5)
58 tolerance = 1e-10;
59 end
61 if (nargin < 6)
62 max_iterations = intmax();
63 end
65 M_of_A = @(A) (1/omega)*diag(diag(A)) + tril(A,-1);
67 if (nargout > 2)
68 [x, iterations, residual_norms] = ...
69 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
70 else
71 [x, iterations] = ...
72 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
73 end
74 end