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
59 max_iterations = intmax();
62 M_of_A = @(A) diag(diag(A));
64 [x, iterations, residual_norms] = ...
65 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
68 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);