iterative/successive_over_relaxation.m

   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   %
  56
  57   if (nargin < 5)
  58     tolerance = 1e-10;
  59   end
  60
  61   if (nargin < 6)
  62     max_iterations = intmax();
  63   end
  64
  65   M_of_A = @(A) (1/omega)*diag(diag(A)) + tril(A,-1);
  66
  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