1 function [x, iterations, residual_norms] = ...
2 jacobi(A, b, x0, tolerance, max_iterations)
8 % iteratively using Jacobi iterations. That is, we let,
10 % A = M - N = D - (L + U)
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,
16 % Ax = (M - N)x = Mx - Nx = b
22 % Thus, our iterations are of the form,
24 % x_{k+1} = M^(-1)*(Nx_{k} + b)
28 % ``A`` -- The n-by-n coefficient matrix of the system.
30 % ``b`` -- An n-by-1 vector; the right-hand side of the system.
32 % ``x0`` -- An n-by-1 vector; an initial guess to the solution.
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``.
39 % ``max_iterations`` -- (optional; default: intmax) the maximum
40 % number of iterations we will perform.
44 % ``x`` -- An n-by-1 vector; the approximate solution to the system.
46 % ``iterations`` -- The number of iterations taken.
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
53 save_residual_norms = false;
55 save_residual_norms = true;
64 max_iterations = intmax();
67 % This creates a diagonal matrix from the diagonal entries of A.
71 % Avoid recomputing this at the beginning of each iteration.
72 b_norm = norm(b, 'inf');
73 relative_tolerance = tolerance*b_norm;
78 rk_norm = norm(rk, 'inf');
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;
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
93 rk_norm = norm(rk, 'inf');
94 if (save_residual_norms)
95 residual_norms = [ residual_norms; rk_norm ];