Add jacobi.m.
[octave.git] / 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 norm of the
36 % residual divided by the norm of ``b``) is less
37 % 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 norms
49 % at each iteration. If not requested, they will
50 % not be computed to save space.
51 %
52 save_residual_norms = false;
53 if (nargout > 2)
54 save_residual_norms = true;
55 residual_norms = [];
56 end
57
58 if (nargin < 4)
59 tolerance = 1e-10;
60 end
61
62 if (nargin < 5)
63 max_iterations = intmax();
64 end
65
66 % This creates a diagonal matrix from the diagonal entries of A.
67 M = diag(diag(A));
68 N = M - A;
69
70 % Avoid recomputing this each iteration.
71 b_norm = norm(b);
72
73 k = 0;
74 xk = x0;
75 rk = b - A*xk;
76
77 while (norm(rk) > tolerance*b_norm && k < max_iterations)
78 x_next = M \ (N*xk + b);
79 r_next = b - A*x_next;
80
81 k = k + 1;
82 xk = x_next;
83 rk = r_next;
84 if (save_residual_norms)
85 residual_norms = [ residual_norms; norm(rk) ];
86 end
87 end
88
89 iterations = k;
90 x = xk;
91 end