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 %
54 if (nargin < 4)
55 tolerance = 1e-10;
56 end
58 if (nargin < 5)
59 max_iterations = intmax();
60 end
62 M_of_A = @(A) diag(diag(A));
63 if (nargout > 2)
64 [x, iterations, residual_norms] = ...
65 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
66 else
67 [x, iterations] = ...
68 classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
69 end
70 end