function [u,k] = nma_SD( A,f,tol,max_iterations) %function solves \$Au=f\$ using the method of steepest descent. %------------------------------------------------------ %file name: nma_SD.m % This function solves Au=f using the method of % steepest descent. % % INPUT: % A: System matrix, must be PSD matrix % f: RHS, column vector % tol: tolerance on the residual ||r||. When the residual (f-Au) becomes % smaller than this value, then the iterative process will stop and % u is returned. % max_iterations: if this number is reached before tolerance is reached % then process is stopped and error is thrown. % % Output: % u: solution to Au=f % k: number of iterations used to reach the solution % % see nma_generate_A_f_poisson_2D_dirichlet_zero_square.m for an % example on how to generate A,f for the 2D poisson pde on unit square % % EXAMPLE: % This example shows using this function to solve the 2D Poisson PDE % problem on unit square % % clear all % myforce = @(X,Y) -exp( -(X-0.25).^2 - (Y-0.6).^2 ); % n = 31; % for grid 33x33 % [A,f]=nma_generate_A_f_poisson_2D_dirichlet_zero_square(n,myforce); % A=full(A); % tol=10^-2; max_iterations=10000; % [u,k]=nma_SD( -A,f,tol,max_iterations) % u=reshape(-u,n,n); % mesh(u) % compare with exact solution % uu=-A\f; % size(uu) % uu=reshape(uu,n,n); % figure % mesh(uu) % % by Nasser M. Abbasi % Matlab 2010 a % Free to use at your own risk if nargin ~= 4 error('nma_steepest_descent_solver:: number of arguments must be 2'); end if not(nma_is_spd_matrix(A)) error('nma_steepest_descent_solver:: A is not SPD'); end [nRowF,nColF] = size(f(:)); if nColF ~= 1 error('nma_steepest_descent_solver:: f must be a vector'); end [n,~] = size(A); if n ~= nRowF error('nma_steepest_descent_solver:: f size not compatible with A size'); end u = zeros(n,1); r = f; k = 0; while (norm(r,2)>tol) && k<=max_iterations fprintf('k=%d, norm=%f\n',k,norm(r,2)); w = A*r; alpha = (r'*r)/(r'*w); u = u + alpha*r; r = r - alpha*w; k = k + 1; end if k>max_iterations error('nma_steepest_descent_solver:: no convergence before reaching max_iterations'); end end