解线性方程组——(Gauss-Seidel)高斯-赛德尔迭代法 | 北太天元

一、Gauss-Seidel迭代法

$n=3$时
A = ( a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33 ) , b = ( b 1 b 2 b 3 ) , A=

,\quad b=, A= ​a11​a21​a31​​a12​a22​a32​​a13​a23​a33​​ ​,b= ​b1​b2​b3​​ ​,

先看一下Jacobi公式的特点

Jacobi公式为
x 1 ( k + 1 ) = b 1 − a 12 x 2 ( k ) − a 13 x 3 ( k ) a 11 x 2 ( k + 1 ) = b 2 − a 21 x 1 ( k ) − a 23 x 3 ( k ) a 22 x 3 ( k + 1 ) = b 3 − a 31 x 1 ( k ) − a 32 x 2 ( k ) a 33

x1(k+1)​=a11​b1​−a12​x2(k)​−a13​x3(k)​​x2(k+1)​=a22​b2​−a21​x1(k)​−a23​x3(k)​​x3(k+1)​=a33​b3​−a31​x1(k)​−a32​x2(k)​​​

可以看出Jacobi方法每次迭代使用的是上一步的结果，每行全部独立计算完后才进入下一轮迭代。

而实际上计算$x_2$时, 本次迭代产生的$x_1$已经更新, 使用$x_3$时， $x_1,x_2$已经更新。由此想法（尽可能使用最新产生的结果），得到Gauss-Seidel方法

Gauss-Seidel公式为
x 1 ( k + 1 ) = b 1 − a 12 x 2 ( k ) − a 13 x 3 ( k ) a 11 x 2 ( k + 1 ) = b 2 − a 21 x 1 ( k + 1 ) − a 23 x 3 ( k ) a 22 x 3 ( k + 1 ) = b 3 − a 31 x 1 ( k + 1 ) − a 32 x 2 ( k + 1 ) a 33

x1(k+1)​x2(k+1)​x3(k+1)​​=a11​b1​−a12​x2(k)​−a13​x3(k)​​=a22​b2​−a21​x1(k+1)​−a23​x3(k)​​=a33​b3​−a31​x1(k+1)​−a32​x2(k+1)​​​

或等价的，将 $A$ 分解为 $A=D-L-U$，其中
$$D=diag(a_{11},a_{22},a_{33}),$$ L = − [ 0 0 0 a 21 0 0 a 31 a 32 0 ] , U = − [ 0 a 12 a 13 0 0 a 23 0 0 0 ] . L=-

,\quad U=-. L=− ​0a21​a31​​00a32​​000​ ​,U=− ​000​a12​00​a13​a23​0​ ​. ( D − L − U ) x = b D x = b + ( L + U ) x D x ( k + 1 ) = b + L x ( k + 1 ) + U x ( k ) ( D − L ) x ( k + 1 ) = b + U x ( k ) x ( k + 1 ) = ( D − L ) − 1 ( b + U x ( k ) )

(D−L−U)x=bDx=b+(L+U)xDx(k+1)=b+Lx(k+1)+Ux(k)(D−L)x(k+1)=b+Ux(k)x(k+1)=(D−L)−1(b+Ux(k))" role="presentation">(D−L−U)xDxDx(k+1)(D−L)x(k+1)x(k+1)=b=b+(L+U)x=b+Lx(k+1)+Ux(k)=b+Ux(k)=(D−L)−1(b+Ux(k))$$\begin{array}{rl}(D-L-U)x& =b\\ Dx& =b+(L+U)x\\ D{x}^{(k+1)}& =b+L{x}^{(k+1)}+U{x}^{(k)}\\ (D-L)x^{(k+1)}& =b+U{x}^{(k)}\\ x^{(k+1)}& =(D-L{)}^{-1}(b+U{x}^{(k)})\end{array}$$

(D−L−U)xDxDx(k+1)(D−L)x(k+1)x(k+1)​=b=b+(L+U)x=b+Lx(k+1)+Ux(k)=b+Ux(k)=(D−L)−1(b+Ux(k))​

其中 $x^{(k+1)}=D^{-1}(b+L{x}^{(k+1)}+U{x}^{(k)})$
写成矩阵的形式
[ x 1 ( k + 1 ) x 2 ( k + 1 ) x 3 ( k + 1 ) ] = [ 1 a 11 1 a 22 1 a 33 ] ( [ b 1 b 2 b 3 ] − [ a 21 a 31 a 32   ] [ x 1 ( k + 1 ) x 2 ( k + 1 ) x 3 ( k + 1 ) ] − [   a 12 a 13 a 23 ] [ x 1 ( k ) x 2 ( k ) x 3 ( k ) ] )

= \left( - -\right) ​x1(k+1)​x2(k+1)​x3(k+1)​​ ​= ​a11​1​​a22​1​​a33​1​​ ​ ​ ​b1​b2​b3​​ ​− ​a21​a31​​a32​​ ​ ​ ​x1(k+1)​x2(k+1)​x3(k+1)​​ ​− ​ ​a12​​a13​a23​​ ​ ​x1(k)​x2(k)​x3(k)​​ ​ ​
下面是其一般形式下的算法

二、算法

$♡$ Gauss-Seidel迭代法

主要思路

输入 : $A,b,x^{(0)}$
输出 : $x$
x ( 0 ) =  initial vector  x ( k + 1 ) = ( D − L ) − 1 ( b + U x ( k ) )

x(0)x(k+1)​= initial vector =(D−L)−1(b+Ux(k))​

添加一些限制

• 容许误差 e_tol,
• 最大迭代步 $N$.

当 残差 $\ge N$ 时，都会停止迭代，输出结果

实现步骤

• 步骤 $1$: $k=0$, $x=x^{(0)}$；

• 步骤 $2$: 计算残差 $r=\Vert b-Ax\Vert$，

  • 如果残差 $r$ >e_tol 且 $k<N$，转步骤 $3$;
  • 否则, 转步骤 $5$；

• 步骤 $3$: 更新解向量$$x=(D-L{)}^{-1}(b+U{x}^{(0)})$$

• 步骤 $4$：$x0=x$, $k=k+1$，转步骤 $2$;

• 步骤 $5$: 输出 $x$.

三、北太天元源程序

Gauss-Seidel

function [x,k,r] = myGS(A,b,x0,e_tol,N)
% Gauss-Seidel迭代法解线性方程组
% Input: A, b(列向量), x0(初始值)
%        e_tol: error tolerant  
%        N: 限制迭代次数小于 N 次
% Output: x , k(迭代次数),r:残差
%   Version:            1.0
%   last modified:      01/29/2024
    n = length(b); k = 0; 
    x=zeros(n,N); % 记录每一次迭代的结果，方便后续作误差分析
    x(:,1)=x0; 
    L = -tril(A,-1); U = -triu(A,1); D = diag(diag(A));
    r = norm(b - A*x,2);
    while r > e_tol && k < N
        x(:,k+2) = inv(D-L)*(b+U*x(:,k+1)); % 不同之处
        r = norm(b - A*x(:,k+2),2); % 残差
        k = k+1;
    end
    x = x(:,2:k+1); % x取迭代时的结果
    
    if k>N
        fprintf('迭代超出最大迭代次数');
    else
        fprintf('迭代次数=%i\n',k);
    end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

保存为 myGS.m 文件

四、数值算例

（此例子与 Jaboci迭代法 文章中的例子相同，可以对比着来看）

$Ax=b$,求 $x$
A = ( 10 − 1 2 0 − 1 11 − 1 3 2 − 1 10 − 1 0 3 − 1 8 ) b = ( 6 25 − 11 15 ) A =

\quad b = A= ​10−120​−111−13​2−110−1​03−18​ ​b= ​625−1115​ ​

用Gauss列主元消去法,得

x = 
   1.000000000000000
   2.000000000000000
  -1.000000000000000
   1.000000000000000
1
2
3
4
5

下面我们用Gauss-Seidel 迭代法进行求解

%% Gauss-Seidel test
% time :      4/24/2024

%% example 1
clc;clear all,format long;
N = 100; e_tol = 1e-8;
A=[10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8];
b=[6; 25; -11; 15];
x0=[0; 0; 0; 0];
    [x11,k1] = myGS(A,b,x0,e_tol,N)
    [x12,k2] = myJacobi(A,b,x0,e_tol,N)
% 作图查看误差变化
    x_exact=[1;2;-1;1]; %真解
    n = length(b);
    error=zeros(n,k1);% 每个分量的误差
        error = abs(x_exact - x11) 
    res =zeros(1,k1); % 残差
    res1 = res;
    res2 = res;
    for i=1:1:k1
        res1(i) = norm(b-A*x11(:,i),2);
        
    end 
    for i=1:1:k2
        res2(i) = norm(b-A*x12(:,i),2);
    end
    % 数值解
        figure(1);
        plot(1:k1,x11(1,:),'-*r',1:k1,x11(2,:),'-og', 1:k1,x11(3,:),'-+b',1:k1,x11(4,:),'-dk');
        legend('x_1','x_2','x_3','x_4');
        title('G-S下每个数值解的变化')
    % 残差变化
        figure(2);
        plot(1:k1,res1,'-*r');
        hold on
        plot(1:k2,res2,'-*b');
        legend('G-S','Jacobi');
        title('残差变化')
    % 误差
        figure(3);
        plot(1:k1,error(1,:),'-*r',1:k1,error(2,:),'-og', 1:k1,error(3,:),'-+b',1:k1,error(4,:),'-dk');
        legend('x_1','x_2','x_3','x_4');
        title('G-S下每个数值解的误差变化')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

运行后得到

和Jacobi相比，其 达到相同精度 所需要得迭代步数更少，如下图

Gauss-Seidel 需进行 10次迭代

而 Jacobi 需进行 26 次迭代