function [x,n]=SOR(A,b,x0,w,M,eps)
%采用超松弛迭代法求線性方程組Ax=b的解
%線性方程組的系數矩陣：A
%線性方程組中的常數向量：b
%迭代初始向量：x0
%松弛因子：w
%解的精度控制：eps
%迭代步數控制：M
%線性方程組的解：x
%求出所需精度的解實際的迭代步數：n
if nargin==4
    eps= 1.0e-6;
    M  = 10000;
elseif nargin == 5
    eps  = 1.0e-6;
end      

if(w<=0 || w>=2)
    error;
    return;
end

D=diag(diag(A));    %求A的對角矩陣
L=-tril(A,-1);      %求A的下三角陣
U=-triu(A,1);       %求A的上三角陣
B=inv(D-L*w)*((1-w)*D+w*U);
f=w*inv((D-L*w))*b;
x=x0;
n=0;                  %迭代次數
tol=1;                %前后兩次迭代結果誤差
while tol>=eps
    x = B*x0+f;       %迭代公式
    n = n+1;
    tol = norm(x-x0);
    x0 = x;
    if(n>=M)
        disp('Warning: 迭代次數太多，可能不收斂！');
        return;
    end
end