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

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

一、Gauss-Seidel迭代法

$n = 3$ 时
$A=\begin{pmatrix} a_{11} & a_{12} &a_{13}\\ a_{21} & a_{22} &a_{23}\\ a_{31} & a_{32} &a_{33}\\ \end{pmatrix} ,\quad b=\begin{pmatrix} b_1\\b_2\\ b_3 \end{pmatrix},$

先看一下Jacobi公式的特点

Jacobi公式为
$\begin{gathered} x_1^{(k+1)} =\frac{b_{1}-a_{12}x_{2}^{(k)}-a_{13}x_{3}^{(k)}}{a_{11}} \\ x_2^{(k+1)} =\frac{b_{2}-a_{21}x_{1}^{(k)}-a_{23}x_{3}^{(k)}}{a_{22}} \\ x_3^{(k+1)} =\frac{b_{3}-a_{31}x_{1}^{(k)}-a_{32}x_{2}^{(k)}}{a_{33}} \end{gathered}$

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

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

Gauss-Seidel公式为
$\begin{aligned} x_1^{(k+1)} &=\frac{b_{1}-a_{12}x_{2}^{(k)}-a_{13}x_{3}^{(k)}}{a_{11}} \\ x_2^{(k+1)} &=\frac{b_{2}-a_{21}x_{1}^{(k+1)}-a_{23}x_{3}^{(k)}}{a_{22}} \\ x_3^{(k+1)} &=\frac{b_{3}-a_{31}x_{1}^{(k+1)}-a_{32}x_{2}^{(k+1)}}{a_{33}} \end{aligned}$

或等价的，将 $A$ 分解为 $A = D - L - U$ ，其中
$D=diag(a_{11},a_{22},a_{33}),$ $L=-\begin{bmatrix} 0 & 0 &0\\ a_{21} & 0&0\\ a_{31} & a_{32} &0\\ \end{bmatrix},\quad U=-\begin{bmatrix} 0 & a_{12} &a_{13}\\ 0 & 0&a_{23}\\ 0 & 0 &0\\ \end{bmatrix}.$ $\begin{aligned} (D-L-U)x &= b\\ Dx &= b+(L+U)x\\ Dx^{(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)}) \end{aligned}$

其中 $x^{(k+1)}=D^{-1}(b+Lx^{(k+1)}+Ux^{(k)})$
写成矩阵的形式
$\begin{bmatrix} x_1^{(k+1)}\\x_2^{(k+1)}\\ x_3^{(k+1)} \end{bmatrix} = \begin{bmatrix} \frac{1}{a_{11}} & &\\ & \frac{1}{a_{22}} &\\ & &\frac{1}{a_{33}}\\ \end{bmatrix}\left( \begin{bmatrix} b_1\\b_2\\ b_3 \end{bmatrix}-\begin{bmatrix} & &\\ a_{21} & &\\ a_{31} & a_{32} &\ \\ \end{bmatrix}\begin{bmatrix} x_1^{(k+1)}\\x_2^{(k+1)}\\ x_3^{(k+1)} \end{bmatrix} -\begin{bmatrix} \ & a_{12} &a_{13}\\ & &a_{23}\\ & &\\ \end{bmatrix}\begin{bmatrix} x_1^{(k)}\\x_2^{(k)}\\x_3^{(k)}\end{bmatrix}\right)$
下面是其一般形式下的算法

二、算法

$\heartsuit$ Gauss-Seidel迭代法

主要思路

输入 : $A,b,x^{(0)}$
输出 : $x$
$\begin{aligned} x^{(0)} &=\text{ initial vector } \\ x^{(k+1)}&=(D-L)^{-1}(b+Ux^{(k)}) \end{aligned}$

添加一些限制

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

当残差 <e_tol 或迭代步数 $\geq N$ 时，都会停止迭代，输出结果

实现步骤

步骤 $1$ : $k = 0$ , $x=x^{(0)}$ ；
步骤 $2$ : 计算残差 $r=\|b-Ax\|$ ，
- 如果残差 $r$ >e_tol 且 $k < N$ ，转步骤 $3$ ;
- 否则, 转步骤 $5$ ；
步骤 $3$ : 更新解向量 $x=(D-L)^{-1}(b+Ux^{(0)})$
步骤 $4$ ： $x 0 = 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/2024n = 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 < Nx(:,k+2) = inv(D-L)*(b+U*x(:,k+1)); % 不同之处r = norm(b - A*x(:,k+2),2); % 残差k = k+1;endx = x(:,2:k+1); % x取迭代时的结果if k>Nfprintf('迭代超出最大迭代次数');elsefprintf('迭代次数=%i\n',k);end
end

保存为 myGS.m 文件

四、数值算例

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

$A x = b$ ,求 $x$
$\begin{pmatrix} 10 & -1 & 2 & 0 \\ -1 & 11 & -1 & 3 \\ 2 & -1 & 10 & -1 \\ 0 & 3 & -1 & 8 \\ \end{pmatrix}\quad b = \begin{pmatrix} 6 \\ 25 \\ -11 \\ 15 \\ \end{pmatrix}$

用Gauss列主元消去法,得

x = 1.0000000000000002.000000000000000-1.0000000000000001.000000000000000

下面我们用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:k1res1(i) = norm(b-A*x11(:,i),2);end for i=1:1:k2res2(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 onplot(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下每个数值解的误差变化')