From: Michael Orlitzky Date: Mon, 6 May 2013 16:29:13 +0000 (-0400) Subject: Add jacobi.m. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=commitdiff_plain;h=b3195a686b04b99c52abe9bcd75125e3747db99f Add jacobi.m. --- diff --git a/iterative/jacobi.m b/iterative/jacobi.m new file mode 100644 index 0000000..3e4f946 --- /dev/null +++ b/iterative/jacobi.m @@ -0,0 +1,91 @@ +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