用MATLAB实现雅克比迭代、高斯-赛德尔迭代以及超松弛迭代

---

前言

近期在上南方科技大学何炳生老师的数值分析课程时，学习了解线性方程组的三种迭代方法，即：雅克比（Jacobi）迭代；高斯-塞德尔（Gauss-Seidel）迭代以及超松弛（SOR）迭代。并按要求使用MATLAB对这三种迭代方法进行了仿真设计。通过仿真设计也对这三种迭代法有了一定的认识，希望通过此贴以通俗的形式分享自己对于这类方法的理解，文中一些不够准确或错误的表达，还望大家指证。

---

一、解线性方程组的方法

求解线性方程组在很多工程领域都是很重要的内容，无论是机器学习也好亦或是其他的工程优化也好。通常来说，对于一个线性方程组$AX=b$，对其的数值解法一般有两类：直接法和迭代法。

• 直接法实际就是通过使用各种线性代数知识在有限步算数运算，计算出该线性方程组精确解的办法（在计算过程中没有舍入误差的情况下），例如本科线性代数中的克拉默法则以及考研线性代数大题最常用的高斯消去法。
• 迭代法则是用某种极限过程去逐步逼近线性方程组的精确解的方法，也就是最后得到的是一个无限逼近与精确解的近似解，这个解只有在迭代无穷次下可以看做是精确解（一般有限步内得不到精确解）。

对于熟悉直接法的人而言，可以发现这种方法的求解过程在计算机下编程实际是较为复杂的，例如：高斯消去法常规的求解过程就要求先构建增广矩阵，而后做初等行变换消元将系数矩阵化为一个上三角矩阵，然后再回代求解出解向量$X$，并且从计算量而言也是较为复杂的，尽管之后变形出了例如三角分解（LU分解）这样的变形解法，但从编程的角度考虑依然是不直观的。
而迭代法，其基本思想是将线性方程组转化为便于迭代的等价方程组，然后基于初始$X$值即$X=x_i^{(0)}=0(i=1,2,3,...n)$，然后按照一定的计算规则，不断地修正得到新的解$X$，最终得到满足精度要求的方程组的近似解。这显然是更加直观且易于编程设计的一种方法。而这其中就有三种最常见的迭代方式，即：

1. 雅克比（Jacobi）迭代
2. 高斯-塞德尔（Gauss-Seidel）迭代
3. 超松弛（SOR）迭代

二、解线性方程组的迭代法及其代码实现

1. 迭代法的收敛性

在直接引入三种迭代法及其MATLAB实现之前，需要先提到一点就是迭代法的收敛性，简单的理解其实就是迭代法并不适用于所有线性方程组，即对于一类方程组迭代时会收敛，近似解会不断逼近精确解，但是对于另外一类方程组则会发散，近似解不断远离精确解。而决定迭代法是否收敛，是否可用于该方程组，主要需要考虑的其实就是系数矩阵$A$。这里直接给出使上述三种迭代法收敛的充分条件：系数矩阵$A$按行（或列）严格对角占优或满足弱对角占优不可约。

（由于本人能力有限，具体的关于于迭代法收敛的充要条件，请参照清华大学出版，李庆阳、王能超、易大义老师编的《数值分析》第五版内容）

2. 基本参数设置

基本参数：矩阵阶数，计算误差，最大迭代次数。其中计算误差和最大迭代次数是迭代终止的判断要求，达到精度或者迭代到最大次数时停止迭代并输出结果。此处对计算误差进行说明，当计算精度满足$\underset{1⩽i⩽n}{\max }\left|x_i^{(k+1)}-x_i^{(k)}\right|<eps$ 时结束迭代。

%==基本参数
%矩阵阶数n
n = 4;
%计算误差
eps = 0.000005;
%最大迭代次数
Nmax = 10000;

矩阵设置：由于以下三种迭代法收敛的充分条件为系数矩阵A满足严格对角占优或按行(或列)满足若对角占优不可约，故在参数设置上设计了两种方法，为方便对比计算结果，之后采用的均是方法一所设置的矩阵。
方法一：手动输入系数矩阵$A$及结果向量$b$

 % 方法一 ： 自行设置矩阵
 A = [-4,1,1,1;
 1,-4,1,1;
 1,1,-4,1;
 1,1,1,-4];
 B = [1,1,1,1];

方法二：根据设置的矩阵阶数，自动生成强对角占优矩阵

 % 方法二：随机生成矩阵
 % 系数矩阵(默认生成一个对角占优矩阵)
 A = 100 * rand(n) - 50;
 for i=1:n
 A(i,i) = sum(abs(A(i,:)))+25 * rand(1); 
 end
 %结果向量
 B = 100 * rand([n 1]);

3. 雅克比（Jacobi）迭代

雅克比迭代实际是将系数矩阵A分解为了三个矩阵的组合，如下图

将其记做$A=D-L-U$，则原方程组 $AX=b$ 等价于 $(D-L-U)X=b$，即
$$DX=(L+U)X+b$$
因为$a_{ii}\ne$