sor.m

辛普森方法

function x=SOR_iterative(A,b,w,x0,kmax)
% This function finds a soltuion to Ax = b by SOR

if nargin < 5
	tol = 1e-6; kmax = 100;             %called by SOR(A,b,w,x0)
elseif kmax < 1
	tol = max(kmax,1e-16); kmax = 1000; %called by SOR(A,b,w,x0,tol)
else 
	tol = 1e-6;                         %called by SOR(A,b,w,x0,kmax)
end

if nargin < 4, x0 = zeros(size(b)); end       %called by SOR(A,b,w)
if nargin < 3, w = 1; end                     %called by SOR(A,b)

n = size(A,1);
x = zeros(size(b));

k=1;
% iteration
while (k<=kmax)
	x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);
	for i=2:n
		x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
    end
	if max(abs(x-x0)) <= tol, break; end
	k=k+1;
	x0=x;
end