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

[octave.git] / iterative / successive_over_relaxation.m

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