gaussseidel.m

% Gauss-Seidel迭代法
% 注释的地方需要修改
clear

% 以下A,b即原方程组化为矩阵形式Ax=b的A,b
A=[5,2,1;-1,4,2;2,-3,10];     % 系数矩阵A，逗号分隔行，分号分隔列
b=[-12;20;3];                 % 向量b
x0=[0;0;0];                   % 迭代初始向量
p=inf;                        % 范数，inf表示无穷范数
ep=1e-4;                      % 精度要求，1e-n表示10的-n次幂
maxiter=500;                  % 最大迭代次数
% 若迭代到最大迭代次数后，误差的p范数仍然达不到精度要求，则判断迭代发散

n=length(b);
x=zeros(n,1);
iter=1;
exitflag=1;
xs(iter,:)=x0;
err(iter,1)=nan;
U=triu(A,1);
L=tril(A);
while exitflag
    x=flipud(triu_solve(rot90(L,2),flipud(b(:)-U*x0(:))));
    xs(iter+1,:)=x;
    err(iter+1,1)=norm(x-x0(:));
    if err(iter+1,1)<=ep
        break
    end
    x0=x;
    iter=iter+1;
    if iter>maxiter
        exitflag=0;
    end
end

if exitflag             % 若迭代发散，不会输出结果
vpa(x,10)               % 线性方程组的解有效数字位数
iter                    % 迭代次数，不需要可删除
vpa(xs,10)              % 完整的迭代序列的有效数字位数，不需要可删除
vpa(err(end),10)        % 最终误差的有效数字位数，不需要可删除
end

function x=triu_solve(varargin)
if nargin==1
    U=varargin{1}(:,1:end-1);
    y=varargin{1}(:,end);
elseif nargin==2
    U=varargin{1};
    y=varargin{2};
end
n=length(y); 
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
    x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);
end
end