第6章 解线性方程组的迭代法(基于MATLAB)
前面我们已经知道对于线性方程组，一般有两种数值解法：直接法和迭代法。直接法前面已经写过了，没看的同学可以移步阅读： 直接法 。本次主要讲述迭代法及其相应的MATLAB代码。
考虑线性方程组
\mathbf{Ax=b}
当 \mathbf{A} 为低阶稠密阵时一般采取直接法进行求解。但是我们在工程技术中经常遇到的是一些大型稀疏矩阵方程组( \mathbf{A} 的阶数很大，但 0 元素较多，解PDE(Partial Differential Equation)时常遇到)，此时利用迭代法进行求解是合适的。在计算机内存和运算两方面，迭代法通常都可利用 \mathbf{A} 中有大量 0 元素的特点。
本文主要介绍 迭代法的基本思路 和 雅可比(Jacobi)迭代法 、 高斯-塞德尔(Gauss-Seidel)迭代法 、 (逐次)超松弛(SOR, Successie over Relaxation)迭代法 、 共轭梯度(CG, Conjugate Gradient)法 。
1. 迭代法的基本思路
对于线性方程组
\mathbf{Ax=b}
为了使用迭代法对其进行求解，我们首先需要构造 迭代序列 ，我们将上式改写为
\mathbf{x}=\mathbf{Bx+f}
其中 \mathbf{B_0} 和 \mathbf{f} 都是已知的，从而我们就可以构造出一个迭代序列为
\mathbf{x}^{(k+1)}=\mathbf{Bx}^{(k)}+\mathbf{f},\qquad k = 0,1,2,\cdots
对向量序列 \{\mathbf{x}^{(k)}\} 一组初始值 \mathbf{x}^{(0)} ，如果上述迭代序列收敛(有课本定理知其充要条件为 \rho(\mathbf{B})<1 )，最终向量序列 \{\mathbf{x}^{(k)}\} 将逐步逼近原方程组的精确解 \mathbf{x}^* 。迭代次数越多，计算误差就越小，因此，任意给定一个精度，都可以通过迭代计算出满足精度要求的解(理论上，实际情况还需考虑计算机存储字长)。
这就是我们构造迭代法的基本思路。MATLAB代码为
function [x,t,it] = Iteration(A,b,I,eps)
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% 输出：
% x: 解矩阵
% t: 时间
% it: 迭代次数
% 迭代初值默认为 0
if nargin < 4
eps = 1e-6;
[n,~] = size(A);
B = zeros(size(A));
x = zeros(n,1);
f = zeros(n,1);
for r = 1:n
for l = 1:n
if l == r
B(r,r) = 0;
B(r,l) = -A(r,l)/A(r,r);
f(r) = b(r)/A(r,r);
x_exact = A\b;
it = 1;
for k = 1:I-1
x = B*x+f;
if norm(x-x_exact)>eps
it = it+1;
t = toc;
示例：
clear
A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36];
I = 1000;
[x,t,it] = Iteration(A,b,I)运行结果
x =
0.016490086000000
it =
16可见只是迭代 22 次，就满足默认精度( 10^{-6} )。
2. 雅可比(Jacobi)迭代法
其实上面第1节的代码使用的方法就是最原始的雅可比迭代法，只是在形式上看起来不那么规范。下面我们来看规范的解 \mathbf{Ax=b} 的雅可比迭代法的形式：
\mathbf{x}^{(0)},\qquad 初始向量，
\mathbf{x}^{(k+1)}=\mathbf{Bx}^{(k)}+\mathbf{f},\qquad k=0,1,\cdots,
其中 \mathbf{B=I-D}^{-1}\mathbf{A}=\mathbf{D}^{-1}(\mathbf{L+U})\equiv\mathbf{J} ， \mathbf{f=D}^{-1}\mathbf{b} ，这里 \mathbf{D,L,U} 分别为对角阵、下三角阵、上三角阵。并称 \mathbf{J} 为解 \mathbf{Ax=b} 的雅可比迭代法的迭代矩阵。
其分量计算公式为
\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)})^\mathrm{T}, x_i^{(k+1)}=\frac{b_i-\sum\limits_{\substack{j=1\\j\ne i}}^n a_{ij}x_j^{(k)}}{a_{ii}},\qquad i=1,2,\cdots,n;k=0,1,\cdots表示迭代次数。
但实际编程计算一般采取简洁且高效的矩阵形式，其MATLAB代码为
function [x,t,it] = Jacobi(A,b,I,eps)
% 雅可比迭代法
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% 输出：
% x: 解矩阵
% t: 时间
% it: 迭代次数
% 迭代初值默认为 0
if nargin < 4
eps = 1e-6;
[n,~] = size(A);
x = zeros(n,1);
D = diag(diag(A)); %求 A 的对角矩阵
L = -tril(A,-1); %求 A 的下三角矩阵,不带对角线
U = -triu(A,1); %求 A 的上三角矩阵
J = D\(L+U);
f = D\b;
x_exact = A\b;
it = 1;
for k = 1:I-1
x = J*x+f;
if norm(x-x_exact)>eps
it = it+1;
t = toc;
仍然计算第一节的例子，只需把函数名改一下：
clear
A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36];
I = 1000;
[x,t,it] = Jacobi(A,b,I)运行结果也与第1节一致：
x =
0.037686548000000
it =
163. 高斯-塞德尔(Gauss-Seidel)迭代法
高斯-塞德尔迭代法与雅可比迭代法很相似，可以看成是其的一种改进。其形式为：
\mathbf{x}^{(0)},\qquad 初始向量，
\mathbf{x}^{(k+1)}=\mathbf{Bx}^{(k)}+\mathbf{f},\qquad k=0,1,\cdots,
其中 \mathbf{B=I-(D-L)}^{-1}\mathbf{A}=\mathbf{(D-L)}^{-1}\mathbf{U}\equiv\mathbf{G} ， \mathbf{f=(D-L)}^{-1}\mathbf{b} ，这里 \mathbf{D,L,U} 分别为对角阵、下三角阵、上三角阵。并称 \mathbf{G} 为解 \mathbf{Ax=b} 的高斯-塞德尔迭代法的迭代矩阵。
其分量计算公式为
\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)})^\mathrm{T}, x_i^{(k+1)}=\frac{b_i-\sum\limits_{j=1}^{i-1} a_{ij}x_j^{(k+1)}-\sum\limits_{j=i+1}^n a_{ij}x_j^{(k)}}{a_{ii}},\qquad i=1,2,\cdots,n;k=0,1,\cdots表示迭代次数。
其MATLAB代码亦采取矩阵形式：
function [x,t,it] = GaussSeidel(A,b,I,eps)
% 高斯-赛德尔(Gauss-Seidel)迭代法
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% 输出：
% x: 解矩阵
% t: 时间
% it: 迭代次数
% 迭代初值默认为 0
if nargin < 4
eps = 1e-6;
[n,~] = size(A);
x = zeros(n,1);
D = diag(diag(A)); %求 A 的对角矩阵
L = -tril(A,-1); %求 A 的下三角矩阵,不带对角线
U = -triu(A,1); %求 A 的上三角矩阵
G = (D-L)\U;
f = (D-L)\b;
x_exact = A\b;
it = 1;
for k = 1:I-1
x = G*x+f;
if norm(x_exact-x)>eps
it = it+1;
t = toc;
同样计算前面的例子，其结果为：
x =
0.054796674000000
it =
8从结果可以看出，对于同样的精度要求，高斯-塞德尔迭代法只需要 8 次迭代(雅可比迭代法需要 16 次)，收敛明显要更快。
3. (逐次)超松弛(SOR)迭代法
简单来说，逐次超松弛迭代法就是在前面高斯-塞德尔迭代法的基础上加了松弛因子 \omega(\omega>0) ，所以解 \mathbf{Ax=b} 的SOR方法为
\mathbf{x}^{(0)} ，初始向量，
\mathbf{x}^{(k+1)}=\mathbf{L_\omega x}^{(k)}+\mathbf{f},\qquad k=0,1,\cdots,
其中 \mathbf{L_\omega}=(\mathbf{D-\omega L})^{-1}((1-\omega)\mathbf{D+\omega U}) ， \mathbf{f}=\omega(\mathbf{D-\omega L})^{-1}\mathbf{b} 。
其分量计算公式为
\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_n^{(0)})^\mathrm{T}, x_i^{(k+1)}=x_i^{(k)}+\omega \frac{b_i-\sum\limits_{j=1}^{i-1} a_{ij}x_j^{(k+1)}-\sum\limits_{j=i}^n a_{ij}x_j^{(k)}}{a_{ii}},\qquad i=1,2,\cdots,n;k=0,1,\cdots\quad\omega为松弛因子。
显然，当 \omega=1 时，SOR方法即为高斯-塞德尔迭代法。当 \omega>1 时，称为 超松弛法 ；当 \omega<1 时，称为 低松弛法 。
其MATLAB代码为
function [x,t,it,w] = SOR(A,b,I,eps,w)
% 逐次超松驰迭代法(successive over relaxation method)迭代法
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% w: 松弛因子(w=1 时即为 Gauss-Seidel 迭代法)
% 输出：
% x: 解矩阵
% t: 时间
% it: 迭代次数
% w: 松弛因子
% 迭代初值默认为 0
[n,~] = size(A);
x = zeros(n,1);
D = diag(diag(A)); %求 A 的对角矩阵
L = -tril(A,-1); %求 A 的下三角矩阵,不带对角线
U = -triu(A,1); %求 A 的上三角矩阵
w_opt = 2/(1+sqrt(1-(vrho(D\(L+U)))^2)); % 最佳松弛因子
if nargin < 4
eps = 1e-6;
w = w_opt;
if nargin < 5
w = w_opt;
Lw = (D-w*L)\((1-w)*D+w*U);
f = w*((D-w*L)\b);
x_exact = A\b;
it = 1;
for k = 1:I-1
x = Lw*x+f;
if norm(x-x_exact)>eps
it = it+1;
t = toc;
还是计算前面那个例子，选择精度为 10^{-6} ，松弛因子 \omega=1.2 ：
clear
A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36];
I = 1000;
[x,t,it,w] = SOR(A,b,I,1e-6,1.2)运行结果为
x =
3.000000000000000
2.000000000000000
1.000000000000000
0.003814113000000
it =
1.200000000000000如果采取默认松弛因子(自动选择最佳松弛因子)：
clear
A = [8 -3 2; 4 11 -1; 6 3 12]; b = [20; 33; 36];
I = 1000;
[x,t,it,w] = SOR(A,b,I)其结果为
x =
3.000000000000000
2.000000000000000
1.000000000000000
0.009697537000000
it =
1.034531942537068此时迭代次数减少为 8 次，松弛因子也选择为接近 1 的数，与高斯-塞德尔迭代法效果基本一致。
4. 共轭梯度(CG)法
也称 共轭斜量法 ，它是一种变分方法，对应于求一个二次函数的极值。CG方法是一种求解大型稀疏对称正定方程组十分有效的方法。
其算法描述为：
（1）任取 \mathbf{x}^{(0)}\in \mathbb{R}^n ，计算 \mathbf{r}^{(0)}=\mathbf{b-Ax}^{(0)} ，取 \mathbf{p}^{(0)}=\mathbf{r}^{(0)} 。
（2）对 k=0,1,\cdots ，计算
\alpha_k=\frac{(\mathbf{r}^{(k)},\mathbf{r}^{(k)})}{(\mathbf{p}^{(k)},\mathbf{Ap}^{(k)})}
\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\alpha_k\mathbf{p}^{(k)}
\mathbf{r}^{(k+1)}=\mathbf{r}^{(k)}-\alpha_k\mathbf{Ap}^{(k)},\quad \beta_k=\frac{(\mathbf{r}^{(k+1)},\mathbf{r}^{(k+1)})}{(\mathbf{r}^{(k)},\mathbf{r}^{(k)})}
\mathbf{p}^{(k+1)}=\mathbf{r}^{(k+1)}+\beta_k\mathbf{p}^{(k)}
（3）若 \mathbf{r}^{(k)}=0 ，或 (\mathbf{p}^{(k)},\mathbf{Ap}^{(k)})=0 ，计算停止，则 \mathbf{x}^{(k)}=\mathbf{x}^* 。由于 \mathbf{A} 正定，故当 (\mathbf{p}^{(k)},\mathbf{Ap}^{(k)})=0 时， \mathbf{p}^{(k)}=\mathbf{0} ，而 (\mathbf{r}^{(k)},\mathbf{r}^{(k)})=(\mathbf{r}^{(k)},\mathbf{p}^{(k)})=0 ，也即 \mathbf{r}^{(k)}=\mathbf{0} 。
写成MATLAB代码为
function [x,t,it] = CG(A,b)
% 共轭梯度法 CG(conjugate gradient)
% 针对大型稀疏对称正定矩阵方程组
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% 输出：
% x: 解矩阵
% t: 时间
% it: 迭代次数
[n,~] = size(A);
x = zeros(n,1);
r = b - A*x;
p = r;
it = 0;
while (sum(r.*r) ~= 0) && (sum(p.*(A*p)) ~= 0)
it = it+1;
rr = sum(r.*r);
alpha = rr/sum(p.*(A*p));
x = x + alpha*p;
r = r - alpha*A*p;
beta = sum(r.*r)/rr;
p = r + beta*p;
t = toc;
用此CG方法求解一个 6 阶Hilbert矩阵：
clear
n = 6; A = zeros(n,n);
% Hilbert 矩阵
for i = 1:n
for j = 1:n
A(i,j) = 1/(i+j-1);
x_star = ones(n,1); b = A*x_star;
[x,t,it] = CG(A,b)运行结果为：
x =
0.999999999999897
1.000000000004002
0.999999999968048
1.000000000092161
0.999999999890764
1.000000000045352
0.006673698000000
it =
350而如果采用Gauss-Seidel迭代法，精度要求为 10^{-6} ，其结果为
x =
0.999999999999871
1.000000000004329
0.999999999967693
1.000000000089720
0.999999999896190
1.000000000042388