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