第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 =
16
3. 高斯-塞德尔(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