iterative/jacobi.m

   1 function [x, iterations, residual_norms] = ...
   2          jacobi(A, b, x0, tolerance, max_iterations)
   3   %
   4   % Solve the system,
   5   %
   6   %   Ax = b
   7   %
   8   % iteratively using Jacobi iterations. That is, we let,
   9   %
  10   %   A = M - N = D - (L + U)
  11   %
  12   % where D is a diagonal matrix consisting of the diagonal entries of
  13   % A (the rest zeros), and N = (L + U) are the remaining upper- and
  14   % lower-triangular parts of A. Now,
  15   %
  16   %   Ax = (M - N)x = Mx - Nx = b
  17   %
  18   % has solution,
  19   %
  20   %   x = M^(-1)*(Nx + b)
  21   %
  22   % Thus, our iterations are of the form,
  23   %
  24   %   x_{k+1} = M^(-1)*(Nx_{k} + b)
  25   %
  26   % INPUT:
  27   %
  28   %   ``A`` -- The n-by-n coefficient matrix of the system.
  29   %
  30   %   ``b`` -- An n-by-1 vector; the right-hand side of the system.
  31   %
  32   %   ``x0`` -- An n-by-1 vector; an initial guess to the solution.
  33   %
  34   %   ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance.
  35   %                    we stop when the relative error (the infinity norm
  36   %                     of the residual divided by the infinity norm of
  37   %                     ``b``) is less than ``tolerance``.
  38   %
  39   %   ``max_iterations`` -- (optional; default: intmax) the maximum
  40   %                         number of iterations we will perform.
  41   %
  42   % OUTPUT:
  43   %
  44   %   ``x`` -- An n-by-1 vector; the approximate solution to the system.
  45   %
  46   %   ``iterations`` -- The number of iterations taken.
  47   %
  48   %   ``residual_norms`` -- An n-by-iterations vector of the residual
  49   %                         (infinity)norms at each iteration. If not
  50   %                         requested, they will not be computed to
  51   %                         save space.
  52   %
  53   save_residual_norms = false;
  54   if (nargout > 2)
  55      save_residual_norms = true;
  56      residual_norms = [];
  57   end
  58
  59   if (nargin < 4)
  60     tolerance = 1e-10;
  61   end
  62
  63   if (nargin < 5)
  64     max_iterations = intmax();
  65   end
  66
  67   % This creates a diagonal matrix from the diagonal entries of A.
  68   M = diag(diag(A));
  69   N = M - A;
  70
  71   % Avoid recomputing this at the beginning of each iteration.
  72   b_norm = norm(b, 'inf');
  73   relative_tolerance = tolerance*b_norm;
  74
  75   k = 0;
  76   xk = x0;
  77   rk = b - A*xk;
  78   rk_norm = norm(rk, 'inf');
  79
  80   while (rk_norm > relative_tolerance && k < max_iterations)
  81     % This is almost certainly slower than updating the individual
  82     % components of xk, but it's "fast enough."
  83     x_next = M \ (N*xk + b);
  84     r_next = b - A*x_next;
  85
  86     k = k + 1;
  87     xk = x_next;
  88     rk = r_next;
  89
  90     % We store the current residual norm so that, if we're saving
  91     % them, we don't need to recompute it at the beginning of the next
  92     % iteration.
  93     rk_norm = norm(rk, 'inf');
  94     if (save_residual_norms)
  95       residual_norms = [ residual_norms; rk_norm ];
  96     end
  97   end
  98
  99    iterations = k;
 100    x = xk;
 101 end