Add classical_iteration() sans comments, refactor jacobi() to use it.

Add gauss_seidel() and successive_over_relaxation() also using classical_iteration().

diff --git a/iterative/classical_iteration.m b/iterative/classical_iteration.m
new file mode 100644 (file)
index 0000000..4e43ff1
--- /dev/null
+++ b/iterative/classical_iteration.m
@@ -0,0 +1,44 @@
+function [x, iterations, residual_norms] = ...
+        classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations)
+
+  save_residual_norms = false;
+  if (nargout > 2)
+     save_residual_norms = true;
+     residual_norms = [];
+  end
+
+  % This creates a diagonal matrix from the diagonal entries of A.
+  M = M_of_A(A);
+  N = M - A;
+
+  % Avoid recomputing this at the beginning of each iteration.
+  b_norm = norm(b, 'inf');
+  relative_tolerance = tolerance*b_norm;
+
+  k = 0;
+  xk = x0;
+  rk = b - A*xk;
+  rk_norm = norm(rk, 'inf');
+
+  while (rk_norm > relative_tolerance && k < max_iterations)
+    % This is almost certainly slower than updating the individual
+    % components of xk, but it's "fast enough."
+    x_next = M \ (N*xk + b);
+    r_next = b - A*x_next;
+
+    k = k + 1;
+    xk = x_next;
+    rk = r_next;
+
+    % We store the current residual norm so that, if we're saving
+    % them, we don't need to recompute it at the beginning of the next
+    % iteration.
+    rk_norm = norm(rk, 'inf');
+    if (save_residual_norms)
+      residual_norms = [ residual_norms; rk_norm ];
+    end
+  end
+
+   iterations = k;
+   x = xk;
+end

diff --git a/iterative/gauss_seidel.m b/iterative/gauss_seidel.m
new file mode 100644 (file)
index 0000000..a844142
--- /dev/null
+++ b/iterative/gauss_seidel.m
@@ -0,0 +1,21 @@
+function [x, iterations, residual_norms] = ...
+        gauss_seidel(A, b, x0, tolerance, max_iterations)
+
+  if (nargin < 4)
+    tolerance = 1e-10;
+  end
+
+  if (nargin < 5)
+    max_iterations = intmax();
+  end
+
+  omega = 1;
+
+  if (nargout > 2)
+    [x, iterations, residual_norms] = ...
+      successive_over_relaxation(A, b, omega, x0, tolerance, max_iterations);
+  else
+    [x, iterations] = ...
+      successive_over_relaxation(A, b, omega, x0, tolerance, max_iterations);
+  end
+end

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
index 347861d5aee459bb10b74c6457ad215c8e32b567..ffee379a20ffddcdfb36b1db09a6ebd8ace83742 100644 (file)
--- a/iterative/jacobi.m
+++ b/iterative/jacobi.m
@@ -50,11 +50,6 @@ function [x, iterations, residual_norms] = ...
   %                         requested, they will not be computed to
   %                         save space.
   %
-  save_residual_norms = false;
-  if (nargout > 2)
-     save_residual_norms = true;
-     residual_norms = [];
-  end
 
   if (nargin < 4)
     tolerance = 1e-10;
@@ -64,38 +59,12 @@ function [x, iterations, residual_norms] = ...
     max_iterations = intmax();
   end
 
-  % This creates a diagonal matrix from the diagonal entries of A.
-  M = diag(diag(A));
-  N = M - A;
-
-  % Avoid recomputing this at the beginning of each iteration.
-  b_norm = norm(b, 'inf');
-  relative_tolerance = tolerance*b_norm;
-
-  k = 0;
-  xk = x0;
-  rk = b - A*xk;
-  rk_norm = norm(rk, 'inf');
-
-  while (rk_norm > relative_tolerance && k < max_iterations)
-    % This is almost certainly slower than updating the individual
-    % components of xk, but it's "fast enough."
-    x_next = M \ (N*xk + b);
-    r_next = b - A*x_next;
-
-    k = k + 1;
-    xk = x_next;
-    rk = r_next;
-
-    % We store the current residual norm so that, if we're saving
-    % them, we don't need to recompute it at the beginning of the next
-    % iteration.
-    rk_norm = norm(rk, 'inf');
-    if (save_residual_norms)
-      residual_norms = [ residual_norms; rk_norm ];
-    end
+  M_of_A = @(A) diag(diag(A));
+  if (nargout > 2)
+    [x, iterations, residual_norms] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
+  else
+    [x, iterations] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
   end
-
-   iterations = k;
-   x = xk;
 end

diff --git a/iterative/successive_over_relaxation.m b/iterative/successive_over_relaxation.m
new file mode 100644 (file)
index 0000000..2bb96f8
--- /dev/null
+++ b/iterative/successive_over_relaxation.m
@@ -0,0 +1,74 @@
+function [x, iterations, residual_norms] = ...
+        successive_over_relaxation(A, b, omega, x0, ...
+                                   tolerance, max_iterations)
+  %
+  % Solve the system,
+  %
+  %   Ax = b
+  %
+  % iteratively using SOR. That is, we let,
+  %
+  %   A = M - N = ((1/omega)*D + L) - U
+  %
+  % where D is a diagonal matrix consisting of the diagonal entries of
+  % A (the rest zeros), and L,U are the upper- and lower-triangular
+  % parts of A respectively. Now,
+  %
+  %   Ax = (M - N)x = Mx - Nx = b
+  %
+  % has solution,
+  %
+  %   x = M^(-1)*(Nx + b)
+  %
+  % Thus, our iterations are of the form,
+  %
+  %   x_{k+1} = M^(-1)*(Nx_{k} + b)
+  %
+  % INPUT:
+  %
+  %   ``A`` -- The n-by-n coefficient matrix of the system.
+  %
+  %   ``b`` -- An n-by-1 vector; the right-hand side of the system.
+  %
+  %   ``x0`` -- An n-by-1 vector; an initial guess to the solution.
+  %
+  %   ``omega`` -- The relaxation factor.
+  %
+  %   ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance.
+  %                    we stop when the relative error (the infinity norm
+  %                     of the residual divided by the infinity norm of
+  %                     ``b``) is less than ``tolerance``.
+  %
+  %   ``max_iterations`` -- (optional; default: intmax) the maximum
+  %                         number of iterations we will perform.
+  %
+  % OUTPUT:
+  %
+  %   ``x`` -- An n-by-1 vector; the approximate solution to the system.
+  %
+  %   ``iterations`` -- The number of iterations taken.
+  %
+  %   ``residual_norms`` -- An n-by-iterations vector of the residual
+  %                         (infinity)norms at each iteration. If not
+  %                         requested, they will not be computed to
+  %                         save space.
+  %
+
+  if (nargin < 5)
+    tolerance = 1e-10;
+  end
+
+  if (nargin < 6)
+    max_iterations = intmax();
+  end
+
+  M_of_A = @(A) (1/omega)*diag(diag(A)) + tril(A,-1);
+
+  if (nargout > 2)
+    [x, iterations, residual_norms] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
+  else
+    [x, iterations] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
+  end
+end