function [x, iterations, residual_norms] = ... jacobi(A, b, x0, tolerance, max_iterations) % % Solve the system, % % Ax = b % % iteratively using Jacobi iterations. That is, we let, % % A = M - N = D - (L + U) % % where D is a diagonal matrix consisting of the diagonal entries of % A (the rest zeros), and N = (L + U) are the remaining upper- and % lower-triangular parts of A. 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. % % ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance. % we stop when the relative error (the norm of the % residual divided by the 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 norms % at each iteration. If not 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; end if (nargin < 5) 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 each iteration. b_norm = norm(b); k = 0; xk = x0; rk = b - A*xk; while (norm(rk) > tolerance*b_norm && k < max_iterations) x_next = M \ (N*xk + b); r_next = b - A*x_next; k = k + 1; xk = x_next; rk = r_next; if (save_residual_norms) residual_norms = [ residual_norms; norm(rk) ]; end end iterations = k; x = xk; end