南理工高等数值分析课程设计张军Word下载.docx
- 文档编号:14163360
- 上传时间:2022-10-19
- 格式:DOCX
- 页数:14
- 大小:228.55KB
南理工高等数值分析课程设计张军Word下载.docx
《南理工高等数值分析课程设计张军Word下载.docx》由会员分享,可在线阅读,更多相关《南理工高等数值分析课程设计张军Word下载.docx(14页珍藏版)》请在冰豆网上搜索。
考虑2次函数,定义为
有如下性质
对一切,
设为(1.1)的解,则
而且对一切,有
则有定理:
设A对称正定,则为(1.1)解的充分必要条件是满足
求,使取最小,这就是求解等价于方程组(1.1)的变分问题。
求解的方法一般是构造一个向量序列,使。
1.1最速下降法
可通过求泛函的极小点来获得式(1.1)的解。
为此,可以从任一出发,沿着泛函在处下降最快的方向搜索下一个近似点,使得在该方向上达到极小值。
由多元微积分知识,在处下降最快的方向是在该点的负梯度方向,经过简单计算,得:
此处,也称为近似点对应的残量。
取,确定的值使取得极小值。
由式(1.4),令
得
从上式求出并记为:
综合上述推导过程,可得最速下降法的算法描述:
给定初始点,容许误差,置。
计算,若,停算,输出作为近似解。
按式(1.6)计算步长因子,置,,转步骤。
最速下降法在理论上是可以收敛的,但是收敛可能会比较缓慢,而且由于舍入误差存在,实际算出的与理论上的最速下降方向可能有很大差别,使计算时出现数值不稳定性。
最速下降法在现代科学与工程中不再具有实用价值。
1.2共轭梯度法(CG)
当线性方程组Ax=b的系数矩阵A是对称正定矩阵,可以采用共轭梯度法对该方程组进行求解,可以证明,式(1.7)所示的n元二次函数
取得极小值点x*是方程Ax=b的解。
共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小化的问题。
从任意给定的初始点出发,沿一组关于矩阵A的共轭方向进行线性搜索,在无舍入误差的假定下,最多迭代n次(其中n为矩阵A的阶数),就可求得二次函数的极小点,也就求得线性方程组Ax=b的解。
其迭代格式为公式(1.8)。
共轭梯度法中关键的两点是迭代格式(21)中最佳步长和搜索方向的确定。
其中可以通过一元函数极小化来求得,其表达式为公式(22);
取,则,要求满足,可得到的表达式(23).
经过一系列的证明和简化,最终可得共轭梯度法的计算过程如下
给定初始点,容许误差,计算,,置。
计算步长因子
置,。
若,停算,输出作为近似解。
计算
置,,转步骤。
1.3预处理共轭梯度法(PCG)
由于计算机存在舍入误差,故前述的共轭梯度法得到的向量会逐渐失去正交性,因而不大可能在n步之内得到原方程组的精确解。
此外,遇到求解大规模的线性方程组,即使能够在n步收敛,收敛速度也不尽如人意。
可以采用预处理技术改善收敛性质,加快收敛速度。
预处理就是对原方程组进行显式或隐式的修正,使得修正后的方程组能够更有效地求解。
即构造一个预处理矩阵M,然后用迭代法求解如下的同解方程组
或
相应得到的算法分别称为左预处理迭代法或右预处理迭代法。
将预处理技术和共轭梯度法结合,可以得到预处理共轭梯度法,算法如下:
预处理共轭梯度法成功与否,关键在于预处理矩阵M是否选得合适。
一个
好的预处理矩阵应该具有如下特征:
M对称正定;
M与A的稀疏性相近;
方程组MSI=rI
易于求解;
的特征值分布比较“集中”,使较小。
预处理矩阵的构造方法有很多,这里介绍本次实验使用的两种方法。
对角预处理法
若A是严格对角占优的矩阵,则可以选择为预处理矩阵,若A是严格分块对角占优矩阵,则可以选择为预处理矩阵;
但是这样简单的选择,往往不能带来很好的加速收敛效果。
矩阵分裂方法
此方法的基本思想是将A分裂为:
;
通过矩阵和来构造预处理矩阵M,一般的做法是取线性稳定迭代法中相应的A的分裂。
比如:
Jacobi分裂(),Gauss-Seidei分裂(,L是A的严格下三角部分),SSOR分裂等,下面介绍SSOR分裂。
将A分裂为,其中,为严格下三角矩阵,它的元素是由A相应部分元素取负号以后构成的。
取:
其中()是参数,一般取为1,本次实验中也取为1。
则预处理矩阵为:
从理论上讲,此预处理矩阵可以使等于的平方根。
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最速下降法的部分结果
图3Matlab显示的工作区
可以看出,使用最速下降法,迭代1821次之后得到收敛的结果。
2.2共轭梯度法的实验结果
用共轭梯度法解同样的线性方程组,截取部分结果如下:
14.8755299
-0.800346351
7.75615781
-13.5156269
2.269183581
-0.365803328
-3.156700613
-1.693023993
-3.510370066
28.94290597
表2共轭梯度法的部分结果
图4Matlab显示的工作区
可以看出,使用共轭梯度法,迭代112次之后得到收敛的结果。
2.3预处理共轭梯度法的实验结果
首先将预处理矩阵M设为矩阵A的对角阵,得到数据如下:
14.87553095
-0.800346018
7.75615762
-13.51562869
2.269182867
-0.365803808
-3.156700618
-1.693025329
-3.510369959
28.94290566
表3M为A的对角阵时预处理共轭梯度法的部分结果
图5Matlab显示的工作区
再用第二种方法:
用SSOR分裂法得到预处理矩阵M,部分结果如下:
14.87552976
-0.800346309
7.756158157
-13.51562743
2.269183319
-0.36580421
-3.156700787
-1.69302531
-3.510369623
28.94290593
表4用SSOR分裂法得到M的预处理共轭梯度法的部分结果
图6Matlab显示的工作区
可以看出,使用预处理共轭梯度法,若预处理矩阵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最速下降法程序
最速下降法主程序
loaddate_A.mat
loaddate_b.mat
[x,iter]=mgrad(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分别为近似解向量和迭代次数
ifnargin<
5,N=1e12;
end
4,ep=1e-6;
3,x=zeros(size(b));
foriter=1:
N
r=b-A*x;
ifnorm(r)<
ep,break;
a=r'
*r/(r'
*A*r);
x=x+a*r;
iter
x
3.3共轭梯度法程序
共轭梯度法主程序
loaddate_A.mat
[x,iter]=mcg(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
%分别为近似解向量和迭代次数
ifnargin<
5,
N=10000;
4,
ep=1e-6;
r=b-A*x;
rr=r'
*r;
ifiter==1
p=r;
else
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 理工 高等 数值 分析 课程设计