线性方程的Seidel迭代法

function X=Seidel(A,b,x0,e)

%  线性方程的Seidel迭代法

%  A为系数矩阵n*n;b为右值n*1;x0为初值n*1,默认为全0;e为精度,默认为10^(-6)

%  BJ为Jacobi迭代矩阵,B1为BJ下三角部分,B2为BJ上三角部分

%  其基本迭代公式为  X = inv(E-B1)*B2*X + inv(E-B1)*g

%  BS为Seidel迭代矩阵，另一基本迭代公式为 X= BS*X + gBS

%  编程过程中要用到的几个迭代收敛的充要条件：

%  1,BS的谱半径小于1

%  2,BJ的第一范式或无穷范式小于1

%  3,系数矩阵A是严格对角占优

%  以上三个条件只要满足其中一个即可使此迭代法对任意的初始向量x0都收敛

%

%  作者：野渡无人

%  最后修改日期：2008.4.9

%

% A=[10 -2 -1;-2 10 -1;-1 -2 5]

% b=[3 15 10]'

%

% >> X=Seidel(A,b)

%  0        0.000000    0.000000    0.000000

%  1        0.300000    1.560000    2.684000

%  2        0.880400    1.944480    2.953872

%  3        0.984283    1.992244    2.993754

%  4        0.997824    1.998940    2.999141

%  5        0.999702    1.999855    2.999882

%  6        0.999959    1.999980    2.999984

%  7        0.999994    1.999997    2.999998

%  8        0.999999    2.000000    3.000000

%  9        1.000000    2.000000    3.000000

%

% X =

%

%    1.0000

%    2.0000

%    3.0000

    [n,m]=size(A);

 

if n~=m

    error('请输入方阵');

end

 

 

if nargin==3

    e=10^(-6);  % 如果只输入3个参数，则默认e为10^(-6)

end

 

if nargin==2

    e=10^(-6);

    x0=zeros(n,1); % 如果只输入2个参数，则默认e为10^(-6),x0为全0列向量

end

 

% 初始化BJ,g,X ,flag,E

BJ=zeros(n);

g=zeros(n,1);

X=zeros(n,1);

flag=0;  % 用于检验是否满足迭代收敛的几个条件

E=eye(n);

 

sn=zeros(n,1);   % 用于判别系数A是否为严格对角占优矩阵

for i=1:n

    tmp=0;

    for j=1:n

        tmp=tmp+abs(A(i,j));

    end

    sn(i,1)=abs(2*A(i,i))-tmp;

end

if min(sn)>0

    %fprintf('A为严格对角占优/n');

    flag=flag+1;   % 满足收敛条件3:系数矩阵A是严格对角占优，即可进行迭代

end

 

%  求迭代矩阵BJ的子过程

for i=1:n

    for j=1:n

        if j==i

            continue

        end

        BJ(i,j)=-A(i,j)/A(i,i);

    end

   g(i)=b(i)/A(i,i);

end

 

if norm(BJ,1)<1

    flag=flag+1;  % 满足收敛条件2:BJ的第一范式小于1，即可以进行迭代

end

 

if norm(BJ,inf)<1

    flag=flag+1;  % 满足收敛条件2:BJ的无穷范式小于1，即可以进行迭代

end

 

B1=tril(BJ);   % 取得BJ的下三角部分

B2=triu(BJ);   % 取得BJ的上三角部分

BS=inv(E-B1)*B2;

gBS=inv(E-B1)*g;

 

lamda=eig(BS);   % lamda为BS的特征值，以向量形式存放

if max(abs(lamda))<1

    flag=flag+1; % 满足收敛条件1:BS的谱半径小于1，即可以进行迭代

end

 

if flag==0

    error('所输入的参数不满足迭代收敛3个条件中的任一个，迭代不能进行！');

end

 

% 格式化输出迭代初值

l=0;  %初始化迭代次数l

fprintf('%3d    ',l);

fprintf('    %f',x0);

fprintf('/n');

 

% 格式化输出第一次迭代值

l=l+1; 

X=BS*x0+gBS;

 

fprintf('%3d    ',l);

for i=1:n

    fprintf('    %f',X(i));

end

fprintf('/n');

 

while max(abs(X-x0))>e 

    x0=X;

    X=BS*x0+gBS;

    %max(abs(X-x0));

    l=l+1;

    fprintf('%3d    ',l);

    for i=1:n

        fprintf('    %f',X(i));

    end

    fprintf('/n');

end