共轭梯度法的预处理ppt课件

1、共轭梯度法的预处理,设 A 为 n 阶对称正定矩阵, 考虑线性方程组 (1) 其中x 是未知向量,b是右端已知向量,一、知识回顾,1、共轭梯度法 如果我们对于固定的 ，在 空间中取 个线性无关的向量 使 满足 当 时， 一定是方程 (1)的精确解，这是共轭梯度法的出发点,一、知识回顾,2、预处理共轭梯度法 （1）预处理共轭梯度法的引进 共轭梯度法理论上是一种直接方法，即按照这个算法得到的 应当是方程 的精确解。但由于其数值的稳定性不好，实际上由于舍入误差的影响 不可能满足原方程。如果把它看成一种迭代法，当A条件数较大时，它收敛的很慢。因此引进了改进共轭梯度法的手段预处理共轭梯度法,一、知识回顾

2、,2）预处理共轭梯度法原理 将Ax=b变形为C-1Ax=C-1b。 令 则 要求 当 是对称正定矩阵时， 其中 是A的最大特征值， 是A的最小特征值，K(A)是A的条件数,一、知识回顾,3）预处理共轭梯度方法 由于预处理矩阵的构造方式很多，所以就形成了多种预处理方法，总的来说，大致可归为如下几种途径： 1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分（如三对角或对角线部分）； 2、通过A的分裂，尤其是线性稳定迭代法中的矩阵 A的分裂构造预处理矩阵（如矩阵分裂法）； 3、通过A的各种近似分解得到预处理矩阵（如不完全分解）； 4、通过矩阵A的多项式构造预处理矩阵； 5、子结构、区域分裂、EBE

3、预处理途径,一、知识回顾,对角预优矩阵法 如果对称正定线性方程组( 1) 的系数矩阵 A 的对角元素相差较大, 则可取预优矩阵C-1为 这是最常用的预处理方阵,SSOR分裂 将A分裂为 ，其中 , 为严格下三角矩阵，则预处理矩阵C为： 其中 是参数，而且有,二、算法,1）共轭梯度法算法 令 对 直到收敛完成 结束,二、算法,2）预处理共轭梯度算法 令 对 直到收敛完成。 结束,三、代码,1）共轭梯度法代码 function x,n,tol=cg(A,b,ep) r=b;p=r;n=0;x=zeros(length(b),1); while n10 alpha=(r*r)/(p*A*p);r1=

4、r; x=x+alpha*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/(r1*r1); p=r+beta*p; n=n+1; end,x = 7.85971307544587 0.422926408294999 -0.0735922390236594 -0.540643016894632 0.0106261628540375 n = 5 tol = 6.745229690695623e-07,2）预处理共轭梯度法代码 以对角预优矩阵为例 clc;clear A=0.2 0.1 1 1 0;0.1 4 -1

5、1 -1;1 -1 60 0 -2;1 1 0 8 4;0 -1 -2 4 700; b=1 2 3 4 5;%SSOR此处程序改为 C=chuli1(A); A1=C*A*C; b1=C*b; ep=10(-2); x,n,tol=cg1(A1,b1,ep,C); function C=chuli1(A) D=diag(diag(A); C=D(-1/2,function C=chuli(A,w) D=diag(diag(A);L=tril(A,-1); C=(D-w*L)*D(-1/2,C=chuli(A,0.009); A1=inv(C)*A*inv(C); b1=inv(C)*b,fu

6、nction x,n,tol=cg1(A,b,ep,C) r=b;p=r;n=0;x=zeros(length(b),1); while n10 alpha=(r*r)/(p*A*p);r1=r; x=x+alpha*C*p; tol=norm(r); if norm(r)ep break; end r=r-alpha*A*p; beta=(r*r)/(r1*r1); p=r+beta*p; n=n+1; end,SSOR此处程序改为 x=x+alpha*inv(C)*p,SSOR分裂 x = 7.85971307544587 0.422926408295008 -0.073592239024

7、0462 -0.540643016894626 0.0106261628540363 m= 4 tol = 6.567520686375984e-05,对角预优矩阵 x= 7.85971307544586 0.422926408295008 -0.0735922390240463 -0.540643016894626 0.0106261628540363 m= 4 tol= 4.731753657562764e-04,此矩阵并没有显示出预处理的优势，所以我们又取矩阵 A=diag(1 10 100 1000 10000)，输出结果如下,五、参考文献,1】蔡大用，白峰杉.现代科学计算.北京：科学出版社，2000