From 40a26bb2dfcf53296d9996b88008a605ad59573c Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 6 May 2013 21:28:32 -0400 Subject: [PATCH] Add classical_iteration() sans comments, refactor jacobi() to use it. Add gauss_seidel() and successive_over_relaxation() also using classical_iteration(). --- iterative/classical_iteration.m | 44 +++++++++++++++ iterative/gauss_seidel.m | 21 ++++++++ iterative/jacobi.m | 45 +++------------- iterative/successive_over_relaxation.m | 74 ++++++++++++++++++++++++++ 4 files changed, 146 insertions(+), 38 deletions(-) create mode 100644 iterative/classical_iteration.m create mode 100644 iterative/gauss_seidel.m create mode 100644 iterative/successive_over_relaxation.m diff --git a/iterative/classical_iteration.m b/iterative/classical_iteration.m new file mode 100644 index 0000000..4e43ff1 --- /dev/null +++ b/iterative/classical_iteration.m @@ -0,0 +1,44 @@ +function [x, iterations, residual_norms] = ... + classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations) + + save_residual_norms = false; + if (nargout > 2) + save_residual_norms = true; + residual_norms = []; + end + + % This creates a diagonal matrix from the diagonal entries of A. + M = M_of_A(A); + N = M - A; + + % Avoid recomputing this at the beginning of each iteration. + b_norm = norm(b, 'inf'); + relative_tolerance = tolerance*b_norm; + + k = 0; + xk = x0; + rk = b - A*xk; + rk_norm = norm(rk, 'inf'); + + while (rk_norm > relative_tolerance && k < max_iterations) + % This is almost certainly slower than updating the individual + % components of xk, but it's "fast enough." + x_next = M \ (N*xk + b); + r_next = b - A*x_next; + + k = k + 1; + xk = x_next; + rk = r_next; + + % We store the current residual norm so that, if we're saving + % them, we don't need to recompute it at the beginning of the next + % iteration. + rk_norm = norm(rk, 'inf'); + if (save_residual_norms) + residual_norms = [ residual_norms; rk_norm ]; + end + end + + iterations = k; + x = xk; +end diff --git a/iterative/gauss_seidel.m b/iterative/gauss_seidel.m new file mode 100644 index 0000000..a844142 --- /dev/null +++ b/iterative/gauss_seidel.m @@ -0,0 +1,21 @@ +function [x, iterations, residual_norms] = ... + gauss_seidel(A, b, x0, tolerance, max_iterations) + + if (nargin < 4) + tolerance = 1e-10; + end + + if (nargin < 5) + max_iterations = intmax(); + end + + omega = 1; + + if (nargout > 2) + [x, iterations, residual_norms] = ... + successive_over_relaxation(A, b, omega, x0, tolerance, max_iterations); + else + [x, iterations] = ... + successive_over_relaxation(A, b, omega, x0, tolerance, max_iterations); + end +end diff --git a/iterative/jacobi.m b/iterative/jacobi.m index 347861d..ffee379 100644 --- a/iterative/jacobi.m +++ b/iterative/jacobi.m @@ -50,11 +50,6 @@ function [x, iterations, residual_norms] = ... % 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; @@ -64,38 +59,12 @@ function [x, iterations, residual_norms] = ... 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 at the beginning of each iteration. - b_norm = norm(b, 'inf'); - relative_tolerance = tolerance*b_norm; - - k = 0; - xk = x0; - rk = b - A*xk; - rk_norm = norm(rk, 'inf'); - - while (rk_norm > relative_tolerance && k < max_iterations) - % This is almost certainly slower than updating the individual - % components of xk, but it's "fast enough." - x_next = M \ (N*xk + b); - r_next = b - A*x_next; - - k = k + 1; - xk = x_next; - rk = r_next; - - % We store the current residual norm so that, if we're saving - % them, we don't need to recompute it at the beginning of the next - % iteration. - rk_norm = norm(rk, 'inf'); - if (save_residual_norms) - residual_norms = [ residual_norms; rk_norm ]; - end + M_of_A = @(A) diag(diag(A)); + if (nargout > 2) + [x, iterations, residual_norms] = ... + classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations); + else + [x, iterations] = ... + classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations); end - - iterations = k; - x = xk; end diff --git a/iterative/successive_over_relaxation.m b/iterative/successive_over_relaxation.m new file mode 100644 index 0000000..2bb96f8 --- /dev/null +++ b/iterative/successive_over_relaxation.m @@ -0,0 +1,74 @@ +function [x, iterations, residual_norms] = ... + successive_over_relaxation(A, b, omega, x0, ... + tolerance, max_iterations) + % + % Solve the system, + % + % Ax = b + % + % iteratively using SOR. That is, we let, + % + % A = M - N = ((1/omega)*D + L) - U + % + % where D is a diagonal matrix consisting of the diagonal entries of + % A (the rest zeros), and L,U are the upper- and lower-triangular + % parts of A respectively. 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. + % + % ``omega`` -- The relaxation factor. + % + % ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance. + % we stop when the relative error (the infinity norm + % of the residual divided by the infinity 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 + % (infinity)norms at each iteration. If not + % requested, they will not be computed to + % save space. + % + + if (nargin < 5) + tolerance = 1e-10; + end + + if (nargin < 6) + max_iterations = intmax(); + end + + M_of_A = @(A) (1/omega)*diag(diag(A)) + tril(A,-1); + + if (nargout > 2) + [x, iterations, residual_norms] = ... + classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations); + else + [x, iterations] = ... + classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations); + end +end -- 2.44.2