数值分析算法 MATLAB 实践 线性方程组迭代法

Gauss-Seidel迭代法

%% 求线性方程组的Gauss-Seidel迭代法，调用格式为
%  [x, k] = guaseidel(A,b,x0,eps,it_max)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-5,
%  it_max 为最大迭代次数，默认为100
%  x 为线性方程组的解，k迭代次数
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x3,k3_cnt] = GuaSeidelFunc(A,b,x0,eps,it_max);
disp('迭代次数:k3_cnt=');
disp(k3_cnt)
disp(['方程组的解：x3 = ']);
disp(x3)


%% 求线性方程组的Gauss-Seidel迭代法，调用格式为
%  [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，ep为精度要求，默认为1e-6,
%  N 为最大迭代次数，默认为500
%  x 为线性方程组的解，k迭代次数
x0=[0,0,0]';%[x1;x2;x3]列向量
it_max = 1000;eps=1e-6;
[x4,k4_cnt]=GuaSeideFunmethod(A,b,x0,eps,it_max);
disp('迭代次数:k4_cnt=');
disp(k4_cnt)
disp(['方程组的解：x4 = ']);
disp(x4)

%% 求线性方程组的GuessSeidel迭代法，调用格式为
%  function [x,k] = GuessSeidel(A,b,eps,it_max)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，eps为精度要求，默认为1e-6,
%  it_max 为最大迭代次数，默认为500
%  x 为线性方程组的解，k迭代次数
it_max = 1000;eps=1e-6;
[x5,k5_cnt] = GuessSeidel(A,b,eps,it_max);
disp('迭代次数:k5_cnt=');
disp(k5_cnt)
disp(['方程组的解：x5 = ']);
disp(x5)

function [x,k] = GuessSeidel(A,b,eps,it_max)
%%GuessSeidel:高斯-赛德尔方法求解线性方程组
%高斯-赛德尔方法是一种迭代法，首先猜测各个xi的初始值（一个简单的方法是设各个xi为0）
%将这些初始值带入到第一个方程解出x1,然后更新x1，将xi带入第二个方程x2，更新x2
%依次迭代，直至数值解非常接近真实值为止
%判断条件：对任意的i，有ea(i)=abs((x(i)-xold(i))/x(i))<=es
%即：max(ea)<=es
%%输入
%A=系数矩阵
%b=右侧矩阵
%es=终止准则（default = 0.00001%）
%maxit=最大迭代次数（default = 500）
%输出：
%x=解向量
%%代码实现
%思路：解向量可以简单的表示为x=d-C*x
%其中di=b_i/a_ii,C的对角线元素为0。
    if nargin<2
        error('至少输入系数矩阵和右侧矩阵')
    end
    if nargin<4||isempty(it_max)
        it_max=500;
    end
    if nargin<3||isempty(eps)
        eps=1e-6;
    end
    [m,n]=size(A);
    if m~=n
         error('系数矩阵必须为方阵')
    end
%求解C
    k = 0;
    C = A;
    for i = 1:n
        C(i,i)=0;
        x(i) = 0;%顺便求初始x
    end
    x=x';
    for i = 1:n
        C(i,1:n) = C(i,1:n)/A(i,i);
    end
    %求解d
    for i=1:n
        d(i)=b(i)/A(i,i);
    end
    %开始迭代
        iter=0;
    while(1)
        xprev=x;%记录上次的x
        for i=1:n
            x(i)=d(i) - C(i,:)*x;%求解并更新xi
            if x(i)~=0
                ea(i)=abs((x(i)-xprev(i))/x(i));
            end
        end
        iter = iter+1;
        if max(ea)=it_max
            break
        end
        k=k+1;
    end
end

function [x,k] = GuessSeidel(A,b,eps,it_max)
%%GuessSeidel:高斯-赛德尔方法求解线性方程组
%高斯-赛德尔方法是一种迭代法，首先猜测各个xi的初始值（一个简单的方法是设各个xi为0）
%将这些初始值带入到第一个方程解出x1,然后更新x1，将xi带入第二个方程x2，更新x2
%依次迭代，直至数值解非常接近真实值为止
%判断条件：对任意的i，有ea(i)=abs((x(i)-xold(i))/x(i))<=es
%即：max(ea)<=es
%%输入
%A=系数矩阵
%b=右侧矩阵
%es=终止准则（default = 0.00001%）
%maxit=最大迭代次数（default = 500）
%输出：
%x=解向量
%%代码实现
%思路：解向量可以简单的表示为x=d-C*x
%其中di=b_i/a_ii,C的对角线元素为0。
    if nargin<2
        error('至少输入系数矩阵和右侧矩阵')
    end
    if nargin<4||isempty(it_max)
        it_max=500;
    end
    if nargin<3||isempty(eps)
        eps=1e-6;
    end
    [m,n]=size(A);
    if m~=n
         error('系数矩阵必须为方阵')
    end
%求解C
    k = 0;
    C = A;
    for i = 1:n
        C(i,i)=0;
        x(i) = 0;%顺便求初始x
    end
    x=x';
    for i = 1:n
        C(i,1:n) = C(i,1:n)/A(i,i);
    end
    %求解d
    for i=1:n
        d(i)=b(i)/A(i,i);
    end
    %开始迭代
        iter=0;
    while(1)
        xprev=x;%记录上次的x
        for i=1:n
            x(i)=d(i) - C(i,:)*x;%求解并更新xi
            if x(i)~=0
                ea(i)=abs((x(i)-xprev(i))/x(i));
            end
        end
        iter = iter+1;
        if max(ea)=it_max
            break
        end
        k=k+1;
    end
end

function [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
% 求线性方程组的Gauss-Seidel迭代法，调用格式为
%  [x,k]=GuaSeideFunmethod(A,b,x0,ep,N)
%  其中, A 为线性方程组的系数矩阵，b 为常数项，ep为精度要求，默认为1e-6,
%  N 为最大迭代次数，默认为500
%  x 为线性方程组的解，k迭代次数
    n=length(b); 
    if nargin<5 
        N=500; 
    end  
    if nargin<4 
       ep=1e-6;  
    end
    if nargin<3
        x0=zeros(n,1);
        k=0; 
    end
    x=zeros(n,1);
    k=0;
    while k<N
        for i=1:n
            if i==1 
          x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);        %开始迭代变量
            elseif i==n  
           x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);  %最后迭代变量
            else                                         %其它迭代变量                    
 x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);
            end
        end
        if  norm(x-x0,inf)<ep
            break;
        end
         x0=x;   
        %disp('x=');
        %disp(x); % 此两行代码可输出中间结果
    k=k+1;
    end
    if k==N
        warning('已到达迭代次数上限！');
    end
 end