function x=wJacobi(A,f,omega,x,M,tau); % WJACOBI: Solves Ax=f using weighted Jacobi iterations % USAGE: % x=WJACOBI(A,f,x,M,tau,lambda); % % INPUTS: % A : nxn matrix % f : nxm matrix % [omega]: relaxation parameter. Default: (2/3) % (omega=1 gives the original Jacobi method) % [x] : nxm matrix of starting values. Default: f % [M] : maximum number of iterations. Default: 1000 % [tau] : convergence tolerence. Default: sqrt(eps) % % OUTPUT: % x : nxm matrix of approximate solutions to Ax=f % % Written by: % -- % John L. Weatherwax 2006-12-19 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- if( nargin<3 | isempty(omega) ) omega=2/3; end if( nargin<4 | isempty(x) ) x=f; end if( nargin<5 | isempty(M) ) M=1000; end if( nargin<6 | isempty(tau) ) tau=sqrt(eps); end % obtain the factors D, L and U in the splitting A = D - L - U: D=diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); % compute the Jacobi iteration matrix: PJ = D \ (L+U); % compute the modified right hand side: omegaF = omega * (D \ f); % compute the weighted Jacobi iteration matrix: wPJ = (1-omega)*eye(size(A)) + omega*PJ; divtol=100*max(abs(f-A*x)); for i=1:M x_new = wPJ * x + omegaF; % the new value r=f-A*x_new; % the new residual dx=x_new - x; % convergence check if( all(abs(r)divtol ) disp('WARNING: Iterations diverging'); return end x = x_new; % reassign the new x end if i==M, disp('WARNING: Maximum iterations exceeded'); end