利用共轭梯度法求解线性方程组

1利用共轭梯度法求解线性方程组翟莹1012205052在自然科学和工程技术中很多问题的解决常常归结为解线性方程组，而这些方程组的系数矩阵大致可分为两种:低阶稠密矩阵和大型稀疏矩阵。而求解方程组的方法通常为直接法和迭代法。直接法用于较低阶方程组的求解，效率较高；迭代法更适用于高阶方程组的求解，常用的经典迭代法有高斯-赛德尔迭代法和雅各比迭代法，但收敛效率较低；共轭梯度法（CG）以较高的收敛速度而经常被采用。从理论上讲，一个n阶方程组最多迭代n步就可求出精确解。1直接法直接法就是经过有限步算术运算，无需迭代可直接求得方程组精确解的方法。但实际计算中由于舍入误差的存在和影响，这种方法也只能得到线性方程组的近似解，该方法是求解低阶稠密矩阵方程组的有效方法。如Cramer法则，Gauss消元法及其变形（LU分解法、Cholesky分解法、QR分解法）等。Matlab中，用矩阵除法“/”或“\”直接求解线性方程组（见附录一），它是一个内部包含着许许多多的自适应算法，对超定方程用最小二乘法求解；对欠定方程因为它的解不唯一，Matlab给出所有解中范数最小的一个特解；对于三对角阵方程组，采用追赶法求解。2迭代法迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法。该方法具有对计算机的存贮单元需求少，程序设计简单、原始系数矩阵在计算过程中不变等优点，是求解大型稀疏矩阵方程组的重要方法。迭代法不是用有限步运算求精确解，而是通过迭代产生近似解逼近精确解。如Jacobi迭代、Gauss-Seidel迭代、SOR迭代法等。在科学研究和大型工程设计中提出的求解问题，经离散后，常常归结为求解形如Ax=b的大型线性方程组，此时系数矩阵A和b是通过测量或其它方法得到的，但是在多数情况下方程可能是病态的，即A的条件数非常大。此时，我们仍然用Matlab中的常规算法，得到的解则可能不是真解。通常情况下，对系数矩阵A和方程的右端b作微小的扰动，再用上述方法求解，扰动后方程组的解与原方程组的解相差甚远。因此，从解决实际问题的角度出发，此时就不能用上述的普通求解法。如果用MATLAB中线性方程组的非定常迭代法进行求解，可能得到非常好的结果。求解大型线性方程组是科学计算中的热点之一，已经有非常多的非定常迭代算法，其中有广义极小残余算法、拟最小残量法、共轭梯度法以及其它各种基于共轭梯度法设计的算法(有PCG、CGS、BICG、LSQR等)。在MATLAB中有gmres，bicg，qmr，cgs2等适合求解大型线性方程组的算法函数，它们最基本的调用格式是x=函数名,Ab。下面以共轭梯度法为例，解低阶线性方程组。3经典的共扼梯度标准算法一般的共扼梯度标准算法是考虑线性方程组Axb（1）的求解问题，其中A为系数矩阵；x为解向量；b为数据向量。具体地说，在应力应变问题的求解中，A表示为刚度矩阵，x一般为结构的应变位移；而b则表示结构的应力。定义1：设nnAR为对称正定矩阵，若向量,npqR满足,0TTApqqAppAq则称,pqA-共轭（有时也称A-正交）。定义2：若12,,,kppp…满足,0ijApp，,0iiApp，,1,2,...,ijk，且ij。则称12,,,kppp…为非零的A-共轭向量组。下面结论显然成立：推论1：若011,,,nppp…为非零的A-共轭向量组，则它构成Rn中的一组基。推论2：若011,,,nppp…为非零的A-共轭向量组且,0jApp，0,1,2,...,1in，则0p。设x为线性代数方程组Axb的解，由推论2，如果011,,,nppp…为非零的A-共轭向量组，则x可由011,,,nppp…线性表示，即存在t0，t1，…，tn-1使得011011...nnxtptptp于是011011...nntAptAptApbiii,,ibptApp，i=0，1，…，n-1可以将上面的结果写成如下的迭代格式：1kkkkxxtp，00x，k=0，1，…（2）共轭梯度法实际上提供了一个构造A-共轭的基011,,,nppp…的方法。它可以看成是由剩余向量（也即目标函数的负梯度向量）kkrbAx共轭正交化的结果。迭代终止的条件为62210kkrbAx。3在共轭梯度法中，取初始向量0x，并令000prbAx。设在迭代过程中已经产生了kx及kp。在下一步的计算中，令1kkkkxxtp其中kt为优化问题0minkktfxtp的解，即k,,kkkkrptApp（3）构造1kp，使得1,0kkApp，其中1kp取如下1kr和kp的线性组合形式11kkkkprp不难算出，1k,,kkkkArpApp（4）综上所述，可得到共轭梯度法算法如下：第一步：任取初始向量0x，计算剩余向量00rbAx，并令00pr，0k。第二步：若r(k)=0，则令kxx，算法终止。否则计算1kkkkxxtp其中kt由式（3）确定。第三步：计算11kkrbAx，并令11kkkkprp，其中1k,,kkkkArpApp第四步：令1kk，转第二步。4实例在应力应变问题的求解中，A表示为刚度矩阵，x一般为结构的应变；而b则表示结构的应力。如下图所示，根据受力计算得到刚度矩阵为160001200800040120001200120002024000800120004800，各节点处的应力为3.3333454024.1667T。分别代入附录一二中，在保留小数点后六位有效数字的情况下，得到相同结果，即T=-0.0003308040.0001274830.0001956670.00512173x，主要因为所求方程组为低阶方程组，且共轭梯度法精度更高。4xy40KN·m3(0,0,0)30KN①20KN2(2,3,4)50KN②③1(0,0,1)10KN/m4(0,0,0)2m2m1m1m附录一一般求解线性方程组方法A=[16000-1200800；04012000-1200；-120002024000；800-120004800]；b=[3.3333454024.1667]'；B=[Ab]；n=4RA=rank（A）RB=rank（B）ifRA==RB&RA==nX=A\bD=null（A，'r'）elsefprintf（'方程组无解'）endvpa（X，6）ans=-0.0003308040.0001274830.0001956670.00512173附录二共轭梯度法求解线性方程组functionx=cg（A，b）tol=1e-6；A=[16000-1200800；04012000-1200；-120002024000；800-120004800]；b=[3.3333454024.1667]'；r=b+A*b；w=-r；z=A*w；5s=w'*z；t=（r'*w）/s；x=-b+t*w；fork=1:numel（b）；r=r-t*z；if（norm（r）tol）return；endB=（r'*z）/s；w=-r+B*w；z=A*w；s=w'*z；t=（r'*w）/s；x=x+t*w；endvpa（x，6）ans=-0.0003308040.0001274830.0001956670.00512173