共轭梯度法的预处理

共轭梯度法的预处理设A为n阶对称正定矩阵,考虑线性方程组(1)其中x是未知向量,b是右端已知向量。Axb一、知识回顾1、共轭梯度法如果我们对于固定的，在空间中取个线性无关的向量使满足当时，一定是方程(1)的精确解，这是共轭梯度法的出发点。knRk12,,,kpppkx12(,,,)(x)min(x)kkSpanpppknnx一、知识回顾2、预处理共轭梯度法（1）预处理共轭梯度法的引进共轭梯度法理论上是一种直接方法，即按照这个算法得到的应当是方程的精确解。但由于其数值的稳定性不好，实际上由于舍入误差的影响不可能满足原方程。如果把它看成一种迭代法，当A条件数较大时，它收敛的很慢。因此引进了改进共轭梯度法的手段——预处理共轭梯度法。nxAxbnx一、知识回顾（2）预处理共轭梯度法原理将Ax=b变形为C-1Ax=C-1b。令则要求当是对称正定矩阵时，其中是A的最大特征值，是A的最小特征值，K(A)是A的条件数~~~AxbAxb~22()()condAcondAA12()()ncondAKA1nbCbxCxACCATT11,~,~一、知识回顾（3）预处理共轭梯度方法由于预处理矩阵的构造方式很多，所以就形成了多种预处理方法，总的来说，大致可归为如下几种途径：1、取预处理矩阵（也称预优矩阵）为A的一个小带宽部分（如三对角或对角线部分）；2、通过A的分裂，尤其是线性稳定迭代法中的矩阵A的分裂构造预处理矩阵（如矩阵分裂法）；3、通过A的各种近似分解得到预处理矩阵（如不完全分解）；4、通过矩阵A的多项式构造预处理矩阵；5、子结构、区域分裂、EBE预处理途径。一、知识回顾对角预优矩阵法如果对称正定线性方程组(1)的系数矩阵A的对角元素相差较大,则可取预优矩阵C-1为这是最常用的预处理方阵nnijnnaAaaadiagDDC)(),,,,(,22112/11SSOR分裂将A分裂为，其中,为严格下三角矩阵，则预处理矩阵C为：其中是参数，而且有TLLCCDA1122(,,)nnDdiagaaa)10(TLLTCDDCDCCA)()()2(11LC2/1)()2(1DCDCL二、算法（1）共轭梯度法算法令对直到收敛完成结束。000000,,xrbAxpr0,1,2,k1111111(r,r)/(p,Ap),,,(r,r),(r,r),kkkkkkkkkkkkkkkkkkkkkkxxprrApprp二、算法（2）预处理共轭梯度算法令对直到收敛完成。结束。010000000,,,,xrbAxzMrpz0,1,k1111111,1111(r,r)/(p,Ap),,,,(zz)/(r,z),kkkkkkkkkkkkkkkkkkkkkkkkxxprrApzMrpzp三、代码（1）共轭梯度法代码•function[x,n,tol]=cg(A,b,ep)•r=b;p=r;n=0;x=zeros(length(b),1);•whilen10•alpha=(r'*r)/(p'*A*p);r1=r;•x=x+alpha*p;•tol=norm(r);•ifnorm(r)ep•break;•end•r=r-alpha*A*p;•beta=(r'*r)/(r1'*r1);•p=r+beta*p;•n=n+1;•endx=7.859713075445870.422926408294999-0.0735922390236594-0.5406430168946320.0106261628540375n=5tol=6.745229690695623e-07（2）预处理共轭梯度法代码以对角预优矩阵为例clc;clearA=[0.20.1110;0.14-11-1;1-1600-2;11084;0-1-24700];b=[12345]‘;%SSOR此处程序改为C=chuli1(A);A1=C*A*C;b1=C*b;ep=10^(-2);[x,n,tol]=cg1(A1,b1,ep,C);functionC=chuli1(A)D=diag(diag(A));C=D^(-1/2);functionC=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;function[x,n,tol]=cg1(A,b,ep,C)r=b;p=r;n=0;x=zeros(length(b),1);whilen10alpha=(r'*r)/(p'*A*p);r1=r;x=x+alpha*C*p;tol=norm(r);ifnorm(r)epbreak;endr=r-alpha*A*p;beta=(r'*r)/(r1'*r1);p=r+beta*p;n=n+1;end%SSOR此处程序改为x=x+alpha*inv(C')*pSSOR分裂x=7.859713075445870.422926408295008-0.0735922390240462-0.5406430168946260.0106261628540363m=4tol=6.567520686375984e-05对角预优矩阵x=7.859713075445860.422926408295008-0.0735922390240463-0.5406430168946260.0106261628540363m=4tol=4.731753657562764e-04方法迭代次数xk误差共轭梯度法5[7.85971307544587,0.422926408294999-0.0735922390236594,-0.5406430168946320.0106261628540375]6.745229690695623e-07对角预优矩阵4[7.859713075445865;0.422926408295008;-0.073592239024046;-0.540643016894626;0.010626162854036]6.567520686375984e-05SSOR分裂4[7.859713075445860;0.422926408295008;-0.073592239024046;-0.540643016894626;0.010626162854036]4.731753657562764e-04此矩阵并没有显示出预处理的优势，所以我们又取矩阵A=diag([110100100010000])，输出结果如下方法迭代次数xk误差共轭梯度法5[1;0.200000000000000;0.030000000000001;0.003999999999562;5.000000000000054e-04]3.895966164939117e-06对角预优矩阵1[1;0.200000000000000;0.030000000000000;0.004000000000000;5.000000000000000e-04]2.775557561562891e-17SSOR分裂1[1;0.200000000000000;0.030000000000000;0.004000000000000;5.000000000000001e-04]1.268840834339045e-17五、参考文献【1】蔡大用，白峰杉.现代科学计算.北京：科学出版社，2000【2】RichardL.Burden,J.DouglasFaires著.数值分析.冯烟利，朱海燕译.北京：高等教育出版社，2005【3】曾喆昭，黄创霞，周富照.数值计算方法与应用.北京：科学出版社，2013