南理工高等数值分析课程设计张军.docx

南京理工大学
课程考核论文

课程名称： 高等数值分析
论文题目： 求解线性方程组Ax=b的极
小化方法比较
姓 名：
学 号：
成 绩：
任课教师评语：

签名：
年 月 日

摘要
本文对求解线性方程组Ax=b的极小化方法进行了理论上的阐述，并且选取三种具有代表性的方法：最速下降法、共轭梯度法（CG）、预处理共轭梯度法（PCG），使用Matlab编程并分别求解相同的线性方程组，在准确性和收敛速度方面进行比较。结果表明，如果预处理矩阵选择得当，使用预处理共轭梯度法（PCG）效果最好。

1、 极小化方法
极小化方法的基本原理是变分原理。
设A对称正定，求解的方程组为
(1.1)
其中，，。考虑2次函数，定义为
(1.2)
有如下性质
⑴ 对一切，
(1.3)
⑵ 对一切，
(1.4)
⑶ 设为(1.1)的解，则

而且对一切，有
(1.5)
则有定理：设A对称正定，则为（1.1）解的充分必要条件是满足

求，使取最小，这就是求解等价于方程组（1.1）的变分问题。求解的方法一般是构造一个向量序列，使。

1.1 最速下降法
可通过求泛函的极小点来获得式（1.1）的解。为此，可以从任一出发，沿着泛函在处下降最快的方向搜索下一个近似点，使得在该方向上达到极小值。
由多元微积分知识，在处下降最快的方向是在该点的负梯度方向，经过简单计算，得：
,
此处，也称为近似点对应的残量。取， 确定的值使取得极小值。由式（1.4），令
,
得

从上式求出并记为：
(1.6)
综合上述推导过程，可得最速下降法的算法描述：
⑴ 给定初始点，容许误差，置。
⑵ 计算，若，停算，输出作为近似解。
⑶ 按式（1.6）计算步长因子，置，，转步骤⑵。
最速下降法在理论上是可以收敛的，但是收敛可能会比较缓慢，而且由于舍入误差存在，实际算出的与理论上的最速下降方向可能有很大差别，使计算时出现数值不稳定性。最速下降法在现代科学与工程中不再具有实用价值。
1.2 共轭梯度法（CG）
当线性方程组Ax=b的系数矩阵A是对称正定矩阵，可以采用共轭梯度法对该方程组进行求解，可以证明，式(1.7)所示的n元二次函数
(1.7)
取得极小值点x*是方程Ax=b的解。共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题。从任意给定的初始点出发，沿一组关于矩阵A的共轭方向进行线性搜索，在无舍入误差的假定下，最多迭代n次（其中n为矩阵A的阶数），就可求得二次函数的极小点，也就求得线性方程组Ax = b的解。其迭代格式为公式(1.8)。
(1.8)
共轭梯度法中关键的两点是迭代格式(21)中最佳步长和搜索方向的确定。其中可以通过一元函数极小化来求得，其表达式为公式(22)；取，则，要求满足，可得到的表达式(23).
(1.9)
(1.10)
经过一系列的证明和简化，最终可得共轭梯度法的计算过程如下
⑴ 给定初始点，容许误差，计算，，置。
⑵ 计算步长因子

置，。
⑶ 若，停算，输出作为近似解。
⑷ 计算

置，，转步骤⑵。
1.3 预处理共轭梯度法（PCG）
由于计算机存在舍入误差，故前述的共轭梯度法得到的向量会逐渐失去正交性，因而不大可能在n步之内得到原方程组的精确解。此外，遇到求解大规模的线性方程组，即使能够在n步收敛，收敛速度也不尽如人意。可以采用预处理技术改善收敛性质，加快收敛速度。
预处理就是对原方程组进行显式或隐式的修正，使得修正后的方程组能够更有效地求解。即构造一个预处理矩阵M，然后用迭代法求解如下的同解方程组
(1.11)
或
(1.12)
相应得到的算法分别称为左预处理迭代法或右预处理迭代法。将预处理技术和共轭梯度法结合，可以得到预处理共轭梯度法，算法如下：
⑴ 给定初始点，容许误差，计算，，置。
⑵ 计算步长因子

置，。
⑶ 若，停算，输出作为近似解。
⑷ 计算

置，，转步骤⑵。
预处理共轭梯度法成功与否，关键在于预处理矩阵M 是否选得合适。一个
好的预处理矩阵应该具有如下特征：
⑴ M 对称正定；
⑵ M 与A 的稀疏性相近；
⑶ 方程组MSI = rI
⑷ 易于求解；
⑸ 的特征值分布比较“ 集中”，使较小。
预处理矩阵的构造方法有很多，这里介绍本次实验使用的两种方法。
⑴ 对角预处理法
若A 是严格对角占优的矩阵，则可以选择为预处理矩阵，若A是严格分块对角占优矩阵，则可以选择为预处理矩阵；但是这样简单的选择，往往不能带来很好的加速收敛效果。
⑵ 矩阵分裂方法
此方法的基本思想是将A 分裂为：；通过矩阵和来构造预处理矩阵M，一般的做法是取线性稳定迭代法中相应的A 的分裂。比如：Jacobi 分裂（），Gauss - Seidei 分裂（，L 是A 的严格下三角部分），SSOR 分裂等，下面介绍SSOR分裂。将A分裂为，其中，为严格下三角矩阵，它的元素是由A 相应部分元素取负号以后构成的。取：
(1.13)
其中（）是参数，一般取为1，本次实验中也取为1。
则预处理矩阵为：
(1.14)
从理论上讲，此预处理矩阵可以使等于的平方根。
2 实验与分析
对于线性方程组，首先使用matlab编写最速下降法程序，并且随机生成的对称正定系数矩阵A和的矩阵b。多次运行之后得到可以求出收敛解的矩阵并保存，则分别使用最速下降法、共轭梯度法、预处理共轭梯度法解决同一个线性方程组。以下是矩阵部分截图，由于矩阵太大，不方便放在论文中，但是生成矩阵的程序在第三章中。

图1 系数矩阵A

图2 矩阵b

2.1 最速下降法的实验结果
由于结果矩阵x为的矩阵，在本文无法展示出截取计算出的最终收敛的结果，故截取前5个和后5个数值如下表：
次数
1
2
3
4
5
值
14.87552312
-0.800348952
7.756162756
-13.51561924
2.269187643
次数
996
997
998
999
1000
值
-0.365809486
-3.156700386
-1.693021867
-3.510366905
28.94290631

表1 最速下降法的部分结果

图3 Matlab显示的工作区
可以看出，使用最速下降法，迭代1821次之后得到收敛的结果。
2.2 共轭梯度法的实验结果
用共轭梯度法解同样的线性方程组，截取部分结果如下：
次数
1
2
3
4
5
值
14.8755299
-0.800346351
7.75615781
-13.5156269
2.269183581
次数
996
997
998
999
1000
值
-0.365803328
-3.156700613
-1.693023993
-3.510370066
28.94290597
表2 共轭梯度法的部分结果

图4 Matlab显示的工作区
可以看出，使用共轭梯度法，迭代112次之后得到收敛的结果。
2.3 预处理共轭梯度法的实验结果
首先将预处理矩阵M设为矩阵A的对角阵，得到数据如下：
次数
1
2
3
4
5
值
14.87553095
-0.800346018
7.75615762
-13.51562869
2.269182867
次数
996
997
998
999
1000
值
-0.365803808
-3.156700618
-1.693025329
-3.510369959
28.94290566
表3 M为A的对角阵时预处理共轭梯度法的部分结果

图5 Matlab显示的工作区
再用第二种方法：用SSOR分裂法得到预处理矩阵M，部分结果如下：
次数
1
2
3
4
5
值
14.87552976
-0.800346309
7.756158157
-13.51562743
2.269183319
次数
996
997
998
999
1000
值
-0.36580421
-3.156700787
-1.69302531
-3.510369623
28.94290593
表4 用SSOR分裂法得到M的预处理共轭梯度法的部分结果

图6 Matlab显示的工作区
可以看出，使用预处理共轭梯度法，若预处理矩阵M取A的对角阵，需迭代123次之后得到收敛的结果；若预处理矩阵M用SSOR分裂法取得，需迭代90次之后得到收敛的结果。

比较上述解线性方程组的三种方法：
①最速下降法法是以残向量作为解的修正方向，收敛速度可能很慢。若对应系数矩阵A的最大、最小特征值，当时，收敛很慢；而且当很小时，因舍入误差影响，计算会不稳定。由实验也可以看出，最速下降法需要的迭代次数最多。
②共轭梯度法属于变分方法的范围，要求系数矩阵A对称正定。使用CG法求解n阶方程组，理论上最多n步便可得到精确解。实验也证明，只需要112步便可得到精确解，比最速下降法快得多。
③ PCG法的使用效果，在于预处理矩阵的选择；若预处理矩阵为系数矩阵的对角阵，用此方法效果可能不尽如人意，比如实验中迭代次数为123次，比CG法还慢些；而使用SSOR法对矩阵进行预处理，再用CG法求解方程组会比较快，实验也可看出，90步迭代后，便可得到达到精确度要求的解。
3 实验用到的matlab程序
3.1 生成随机矩阵的程序
N=1000;
D = diag(rand(N,1));
U = orth(rand(N,N));
B = U' * D * U ;
save('date_A.mat','B');
b = rand(N,1);
3.2 最速下降法程序
A 最速下降法主程序
load date_A.mat
load date_b.mat
[x,iter] = mgrad(B,b);
B 最速下降法模块
%最速下降法
function [x,iter]=mgrad(A,b,x,ep,N)
%用途：用最速下降法解线性方程组Ax=b
%格式：[x,iter]=mgrad(A,b,x,ep,N) 其中A为系数矩阵，b为右端向量，x为初始向量
%（默认零向量），ep为精度（默认 1e-6），N为最大迭代次数（默认1000次），返回参数
%x，iter分别为近似解向量和迭代次数
if nargin<5, N=1e12;end
if nargin<4,ep=1e-6;end
if nargin<3,x=zeros(size(b));end
for iter=1:N

r=b-A*x;
if norm(r)<ep,break;end
a=r'*r/(r'*A*r);
x=x+a*r;
iter
x
end

3.3 共轭梯度法程序
A 共轭梯度法主程序
load date_A.mat
load date_b.mat
[x,iter]=mcg(B,b)
B 共轭梯度法模块
function [x,iter]=mcg(A,b,x,ep,N)
%用途：用共轭梯度法解线性方程组Ax=b
%格式：[x,iter]=mcg(A,b,x,ep,N)其中A为系数矩阵，b为右端向量，x为初始向量（默认
%零向量）ep为精度（默认1e-6），N为最大迭代次数（默认10000次），返回参数x，iter
%分别为近似解向量和迭代次数
if nargin <5,
N=10000;
end
if nargin<4,
ep=1e-6;
end
if nargin<3,x=zeros(size(b));
end
r=b-A*x;
for iter=1:N
rr=r'*r;
if iter==1
p=r;
else
beta=rr/rr1;
p=r+beta*p;
end
q=A*p;
alpha=rr/(p'*q);
x=x+alpha*p;
r=r-alpha*q;
rr1=rr;
if (norm(r)<ep),break;
end
end
3.4 预处理共轭梯度法程序
A 得到预处理矩阵M的程序
① M为A的对角矩阵
load date_A.mat
M=diag(diag(B));
save('date_juzhenM1.mat','M');
② M由SSOR分裂法得到
load date_A.mat
D=diag(diag(B));
CL=tril(-B);
L=((D-CL)*D)^1/2;
M=L*L';
save('date_juzhenM.mat','M');
B 预处理共轭梯度法主程序
load date_A.mat
load date_b.mat
[x,iter]=mpcg(B,b)
C 预处理共轭梯度法模块
function [x,iter]=mpcg(A,b,x,ep,N)
%用途：用预处理共轭梯度法解线性方程组Ax=b
%格式：[x,iter]=mpcg(A,b,x,ep,N)其中A为系数矩阵，b为右端向量，x为初始向量（默认
%零向量）ep为精度（默认1e-6），N为最大迭代次数（默认10000次），返回参数x，iter
%分别为近似解向量和迭代次数
load date_juzhenM.mat
L=inv(M);
if nargin <5,
N=10000;
end
if nargin<4,
ep=1e-6;
end
if nargin<3,x=zeros(size(b));
end
r=L*(b-A*x);
for iter=1:N
rr=r'*r;
if iter==1
p=r;
else
beta=rr/rr1;
p=r+beta*p;
end
q=L*A*p;
alpha=rr/(p'*q);
x=x+alpha*p;
r=r-alpha*q;
rr1=rr;
if (norm(r)<ep),break;
end
end