《武汉大学高等工程数学课件演示文档》由会员分享，可在线阅读，更多相关《武汉大学高等工程数学课件演示文档（362页珍藏版）》请在金锄头文库上搜索。
1、第第 2 2 章章解线性代数方程组的解线性代数方程组的直接法与迭代法直接法与迭代法.引言在自然科学和工程技术中，很多问题的解决常常归结为解线性代数方程组。例如:船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解的问题，用差分法或者有限元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组，且常常为求解大型线性方程组。.引言与上述这些问题有关的计算方法包括：求解线性方程组的数值方法与计算矩阵的特征求解线性方程组的数值方法与计算矩阵的特征值及特征向量的数值方法。值及特征向量的数值方法。 在本章中，我们将注意力集中在线性方程组求在本章中，我们将注意力集中在
2、线性方程组求解的直接方法和迭代方法。解的直接方法和迭代方法。也就是我们要学习的求解线性方程组的两类方法。.引言n线性方程组的这两类数值解法有以下特点。n直接法：经过有限步算术运算，可求得方程组的精直接法：经过有限步算术运算，可求得方程组的精确解的方法（若在计算过程中没有舍入误差）；确解的方法（若在计算过程中没有舍入误差）；n迭代法：用某种极限过程去逐步逼近线性方程组精迭代法：用某种极限过程去逐步逼近线性方程组精确解的方法；确解的方法；n迭代法具有占用存储单元少，程序设计简单，原始系数矩阵在迭代过程中不变等优点，但存在收敛性及收敛速度等问题。n下面由我们从小熟知的消元法开始.2.1高斯消元法n设
3、线性方程组n简记为AX=b（矩阵形式/向量形式）.高斯消元法n如果AX=b是n阶线性代数方程组，则.高斯消元法n克莱姆法则在理论上有着重大意义，但克莱姆法则在理论上有着重大意义，但在实际应用中存在很大的困难，在线性在实际应用中存在很大的困难，在线性代数中，为解决这一困难给出了高斯消代数中，为解决这一困难给出了高斯消元法。元法。n我们解释一下为何我们解释一下为何存在很大的困难存在很大的困难.看来我们确实有必要学习方程组的求解方法，从一个简单的例子开始nGauss消去法的计算量消去法的计算量这是我们介绍完这是我们介绍完Gauss消消去法之后，对其计算量的分析，大家先看看去法之后，对其计算量的分析，
4、大家先看看n高斯法解高斯法解n n阶线性方程组需要的总乘除法次数为阶线性方程组需要的总乘除法次数为n注注1 1：当：当n n = 20 = 20时约为时约为26702670次，比克莱姆法则次，比克莱姆法则9.79.7 10102020次次大大减少。大大减少。n注注2 2：世界上最快的超级计算机每秒千万亿次的浮点运算：世界上最快的超级计算机每秒千万亿次的浮点运算能力能力, ,即：即：10101515。n注注3 3：目前最好的台式或笔记本电脑每秒万亿次的浮点运：目前最好的台式或笔记本电脑每秒万亿次的浮点运算能力，即：算能力，即：10101212 。n请同学们估算一下，用克莱姆法则解请同学们估算一下
5、，用克莱姆法则解2020阶方程组各需要耗阶方程组各需要耗时多少？！时多少？！ n超级计算机超级计算机9.79.7 101020 20 / 10/ 1015 15 /3600 =2.8/3600 =2.8小时小时n台式或笔记本电脑台式或笔记本电脑9.79.7 101020 20 / 10/ 1012 12 /3600 =2800/3600 =2800小时小时另外，还有存储空间相关的空间复杂性问题另外，还有存储空间相关的空间复杂性问题.例题n例1.用消元法解方程组.例题n第一步：-2x（1）+(3)得.例题n第二步：1x（2）+（4）n解方程组（1）（2）（5）得：x=1,2,3T.n这样，一般方
6、程组的这样，一般方程组的求解问题就转化为一求解问题就转化为一个上三角形方程组的个上三角形方程组的求解问题。当然，如求解问题。当然，如果你愿意，同样也可果你愿意，同样也可以转化为一个下三角以转化为一个下三角形方程组的求解问题形方程组的求解问题。总结n下面我们就来看看三下面我们就来看看三角形方程组的求解方角形方程组的求解方法！法！.2.1.1高斯顺序消元法n对下三角形方程的求解设（1）.高斯顺序消元法n由（1）得.n对上三角方程组设.n由（2）式回代得.一般线性代数方程组解法的思考n通过例子，我们得出了一般线性代数方程通过例子，我们得出了一般线性代数方程组的求解可以转化为上或下三角形方程组组的求解
7、可以转化为上或下三角形方程组来求解的结论，上面我们又研究了上或下来求解的结论，上面我们又研究了上或下三角形方程组的求解方法，下面是讨论一三角形方程组的求解方法，下面是讨论一般线性代数方程组求解的高斯消去法的时般线性代数方程组求解的高斯消去法的时候了！。候了！。.高斯顺序消去法n对一般线性代数方程组设Ax=b.n记A(1)=A b(1)=b1、第一次消元。设.高斯顺序消去法.高斯顺序消去法n设第k-1次消元得A(k)x=b(k)其中下面进行第下面进行第k k次消元次消元.高斯顺序消去法第k次消元.高斯顺序消去法n最后,在完成n-1次消元后得到如下上三角形矩阵n这样，一般方程组就变成了所谓的三角形
8、方程组，问题也就转化为三角形方程组的求解问题。.高斯顺序消去法n上述过程，也就是对于方程组AX=b系数矩阵做如下变换：得到上三角方程组.高斯顺序消去法.高斯顺序消去法.高斯顺序消去法注注意意：算算法法对对系系数数矩矩阵阵对对角角元元素素的的非非零零要要求求.高斯顺序消去法.高斯顺序消去法算法框图高斯顺序消去法框图.高斯消去法的计算量.nGauss消去法的计算量消去法的计算量n以乘除法的次数为主以乘除法的次数为主n(1) 消元过程：消元过程：n第第k步时步时(n k) + (n k) (n k + 1) = (n k) (n k + 2)n共有共有.nGauss消去法的计算量消去法的计算量n(2
9、) 回代过程：回代过程：n求求xi中，乘中，乘ni次，除次，除1 1次，共次，共ni+1+1次（次（i = 1，n1）n共有共有.nGauss消去法的计算量消去法的计算量n(3) (3) 高斯法解高斯法解n n阶线性方程组需要的总乘除法次阶线性方程组需要的总乘除法次数为数为n注注1 1：当：当n n = 20 = 20时约为时约为26702670次，比克莱姆法则次，比克莱姆法则9.79.7 10102020次大大减少。次大大减少。n注注2 2：世界上最快的超级计算机每秒万万亿次的：世界上最快的超级计算机每秒万万亿次的浮点运算能力浮点运算能力, ,即：即：10101616。n注注3 3：目前最好
10、的台式或笔记本电脑每秒万亿次：目前最好的台式或笔记本电脑每秒万亿次的浮点运算能力，即：的浮点运算能力，即：10101212 。n请同学们估算一下，用克莱姆法则解请同学们估算一下，用克莱姆法则解2020阶方程组阶方程组各需要耗时多少？！各需要耗时多少？！.高斯顺序消去法条件.n进一步的考虑进一步的考虑解决办法解决办法n(1) (1) 若消元过程中出现若消元过程中出现akk( (k) ) = 0 = 0，则无，则无法继续法继续n(2) (2) 若若akk( (k) ) 0 0，但较小，则小主元做，但较小，则小主元做除数将产生大误差除数将产生大误差n(3) (3) 每次消元都选择绝对值最大者作主元，
11、每次消元都选择绝对值最大者作主元，称为高斯主元消去法称为高斯主元消去法n(4) (4) 通常第通常第k步选择步选择n第第k列主对角元以下元素绝对值最大列主对角元以下元素绝对值最大者作主元（该行与第者作主元（该行与第k行对调），称为列行对调），称为列主元消去法。主元消去法。.2.1.2高斯主元素消去法nGauss列主元消元法n从第一列中选出绝对值最大的元素，如果，则n交换第一行n和第i行，即交换.高斯列主元消去法算法.高斯列主元消去法n第k步从的第k列，中选取绝对值最大项，记录所在行，即若交换第k行与l行的所有对应元素，再进行顺序消元。下面，详细描述一下列选主元的高斯消去法.高斯列主元消去法.高
12、斯列主元消去法框图.n注意到算法选主元素的第注意到算法选主元素的第2步，是非正常结束，步，是非正常结束，即即n出现上述问题的原因是在某步消元过程之前的出现上述问题的原因是在某步消元过程之前的列选主元失败，无法选到合适的非零的适当大列选主元失败，无法选到合适的非零的适当大的主元素。的主元素。n解决办法，在更大的范围，即余下的解决办法，在更大的范围，即余下的n-k阶子阶子阵中来选择主元素。阵中来选择主元素。n在第在第k到到n列和第列和第k k到到n n行的子阵中选元素绝对行的子阵中选元素绝对值最大者作主元，称为全主元消去法。值最大者作主元，称为全主元消去法。.2.全主元消去法n先看一个例子.求解方
13、程组.全主元消去法.全主元消去法.全主元消去法很明显，用列选主元法所得到很明显，用列选主元法所得到的解要比顺序高斯消去法好，的解要比顺序高斯消去法好，但仍然有较大的误差，下面介但仍然有较大的误差，下面介绍所谓的全选主元高斯消去法：绍所谓的全选主元高斯消去法：.全主元消去法.全主元消去法.全主元消去法.全主元消去法.从此可以看出，全主元法与列主元相比，其最大的差别是需要根据列的交换情况，按照相反的顺序将解的分量的顺序重新调整！全主元消去法.全主元消去法全主元消去法的严格算法描述如下：.Gauss全主元消元算法.Gauss全主元消元算法.Gauss全主元消元算法.Gauss全主元消元算法.小结n比
14、较而言,Gauss顺序消去法条件苛刻，且数值不稳定；Gauss全主元消去法工作量偏大，需要比较个元素及行列交换工作，算法复杂。因此从算法优化的角度考虑，Gauss列主元消去法比较好。.思考n在上述讨论中，Gauss顺序消去法、Gauss列主元消去法、Gauss全主元消去法的表示方法或算法描述过程中主要采用的是分量的形式。n但本章一开始，我们就漂亮地将分量形式的方程组表示为了矩阵的形式“AX=b”，而从上述过程来看，无论是消元过程，还是回代求解，都可以同样漂亮地表示为矩阵的形式。n也就是说，至少在形式上方程组的求解问题可以借助矩阵运算以及相关的理论、方法。下面就来探讨方程组求解的矩阵三角分解方法
15、。首先，看看矩阵的三角分解！.2.2矩阵的三角分解法n我们知道，对矩阵进行一次行的初等变换，就相当于用相应的初等矩阵去左乘原来的矩阵。因此我们用矩阵乘法来考察Gauss消元法。n即可得到求解线性方程组的另一种直接法：矩阵的三角分解法。.n三角分解法解方程的原理三角分解法解方程的原理n若A = LU，则Ax = b LUx = b n其中L、U为三角形矩阵，则方程组易解n定义：n(1) 称 为单位下三角阵n(2) 设L为单位下三角阵，U为上三角阵，若A = LU，则称A可三角（LU）分解，并称A = LU为A的三角分解（杜利特尔分解）n(3) A也可以分解为一个下三角阵L和一个单位上三角阵U，此
16、时， A的三角分解（LU）被称为克劳特分解。.三角分解可顺利进行的条件三角分解可顺利进行的条件n定理：(1) A = (aij)nn有唯一LU分解 A的顺序主子式k 0，k = 1，2，.，n 1n(2) 若A = (aij)nn可逆，则存在置换阵P，使PA的n个顺序主子式全不等于0n注：由Ax = b PAx = Pb知，也就是说，经适当行交换后，A总可三角分解.n矩阵三角分解过程n设A的各阶顺序主子式不为0，且A=LU，即n(1) L阵主对角线（含）上边：当列标m 行标k时，lkm = 0n j = k，k+1，.，n；k = 1，2，.，nn(2)U阵主对角线下边：当行标m 列标k时，u
17、km = 0n i = k+1，k+2，.，n；k = 1，2，.，n1.矩阵三角分解计算公式n设A的各阶顺序主子式不为0，且A=LU，即n(1) 主对角线（含）上边：当列标m 行标k时，lkm = 0 n(2) 主对角线下边：当行标m 列标k时，ukm = 0n欲求欲求lik和和ukj.n矩阵三角分解计算公式-k=1时n(1) 主对角线（含）上边：当列标m 行标k时，lkm = 0n j=k,k+1, .,n；k=1,2,.,nn(2) 主对角线下边：当行标m 列标k时，ukm = 0n i=k+1,k+2,.,n；k=1,2,.,n1n欲求lik和ukjn设k = 1，那么，就可以计算L的
18、第一行和U的第一列，有：a1j = l11u1j u1j = a1j j = 1，2，.，n（ U的第一行）nai1 = li1u11 li1 = ai1/u11 i = 2，3，.，n（ L的第一列）.n矩阵三角分解计算过程-逐行（列）计算L、Un一般地，逐行（列）计算L、Un最后：u11u12.u1n第1步l21u22.u2n第2步l31l32.ln1ln2unn第n步.n三角分解后求解两个三角形方程组三角分解后求解两个三角形方程组n下三角下三角 Ly = = b：n上三角上三角 Rx = = y：n紧凑格式：紧凑格式：.事实上，上述过程就是高斯消去法的矩阵形式事实上，上述过程就是高斯消去
19、法的矩阵形式.第2步.第n-1步完成后，就可以将A变换为上三角阵U.所以，A就有三角分解LU.一个例子，n=3的情况.Doolittle分解的一个例子，n=3的情况.Doolittle分解，一般情况.Doolittle分解（L第一列，U第一行）.Doolittle分解（L第2列，U第2行）.Doolittle分解（L第k列，U第k行）.Doolittle分解（L第k列，U第k行）.Doolittle分解,一般计算公式.Doolittle分解后求解两个三角形方程组后求解两个三角形方程组.例题.例题.例题.例题.Doolittle分解算法.Doolittle分解算法.选主元的三角分解法选主元的三角
20、分解法n针对方程组Ax=b，矩阵A的三角分解能顺利进行的条件是A的各阶顺序主子式不等于零。同时，我们也已经指出，三角分解法实质上是高斯消去法的矩阵表示。n因此，高斯消去法求解方程组时存在的数值稳定性问题（除数为零或很小），三角分解法同样存在。n为此，有必要与高斯消去法一样进行主元的选择，以便能保证三角分解的顺利进行，使三角分解法具有好的数值稳定性。因而就有了选主元的三角分解法，只需要在直接三角分解法中引入选主元的算法，就可以实现选主元的三角分解法。这里就不详细讨论了！下面针对一类特下面针对一类特殊的方程组，介殊的方程组，介绍相应的三角分绍相应的三角分解法解法.2.2.3平方根法和改进的平方根法
21、-对称矩阵的Cholesky（乔列斯基）分解n在科学与工程计算中，线性方程组大多数的系数矩阵具有对称正定对称正定的性质，因此基于对称正定矩阵这一特性的三角分解方法是求解对称正定方程组的一种有效方法。n由于系数矩阵A的良好性质，分解过程无需选主元，具有良好的数值稳定性。首先，我们讨论对称正定矩阵的三角分解，即Cholesky（乔列斯基）分解。.对称矩阵的Cholesky（乔列斯基）分解nA对称：AT=AA正定：A的各阶顺序主子式均大于零。即.对称矩阵的Cholesky分解n由Doolittle分解，对称正定矩阵A有唯一分解.对称矩阵的Cholesky分解n定理2.2.4设A为对称正定矩阵，则存在
22、唯一分解A=LDLT，其中L为单位下三角阵，D=diag(d1,d2,dn)且di0(i=1,n).对称矩阵的Cholesky分解n证明：.对称矩阵的Cholesky分解所以，非奇异的上三角阵U1可化为.对称矩阵的Cholesky分解.对称矩阵的Cholesky分解.对称矩阵的Cholesky分解推论：设A为对称正定矩阵，则存在唯一分解其中L为具有主对角元素为正数的下三角矩阵。LLT称为对称正定矩阵称为对称正定矩阵A的的Cholesky分解。分解。.对称矩阵的Cholesky分解n证明：证明：推论：设A为对称正定矩阵，则存在唯一分解其中L为具有主对角元素为正数的下三角矩阵。.Cholesky分
23、解的求法.Cholesky分解的求法.Cholesky分解的求法.Cholesky分解法Cholesky分解法缺点及优点 优点：可以减少存储单元。 缺点：存在开方运算，可能会出现根号下负数，同时计算量也大增。怎么改进？！怎么改进？！.记得我们在由LDLT获得LLT分解时，人为地将D分为D1/2* D1/2，这是导致平方根法计算量大等问题的根源。因而，改进的想法首先就是避免这种分割。具体地：Cholesky分解法，怎么改进？！怎么改进？！.改进Cholesky分解法n基于分解A=LDLT，其中.改进的cholesky分解.改进的cholesky分解.例题.例题.例题.例题.改进平方根法nA=LD
24、LT分解（对应改进平方根法），既适合于解对称正定方程组，也适合求解A为对称，而各阶顺序主子式不为零的方程组n而对A=LLT分解（对应平方根法）只适合于对称正定方程组n上面，就是求解系数矩阵对称正定的方程组的平方根法和改进平方根法。下面针对另外一类特殊方程组求解的方法：追赶法！.2.2.4三对角方程组求解的追赶法.三对角方程组求解的追赶法.三对角方程组求解的追赶法.三对角方程组求解的追赶法.三对角方程组求解的追赶法n其计算工作量为5n-4次乘除法（包括三角分解2n-2，求解上三角方程n-1,回代2n-1）。工作量小.n其实现的条件为qi不为零。有以下定理可得三对角方程组求解的充分性条件。.三对角
25、方程组求解的追赶法.解三对角矩阵线性方程组的追赶法程序框图.从三对角矩阵带状矩阵：思考如果方程组的系数矩阵的上半带宽为s、下半带宽为r,即有如下三角分解的结果.从三对角矩阵带状矩阵：思考有兴趣的同学请自己针对你在科学研究中碰到的实际问题，分析设计相应的三角分解算法.Cramer ruleGauss eliminationLU factorizationGauss-Jordan eliminationSquare root /improved square root追赶法(n+1)!n3/3n3/3n3/2n3/65n-4Row/column/complete PivotingOnly elim
26、inating elements in a column below the diagonal oneNo pivotingDirectly factorization column pivotingEliminating elements in a row except for the diagonal elementA symmetric and positive definiteNo pivotingA三对角，弱对角占优直接法时间复杂性直接法时间复杂性.解线性代数方程组的小结n高斯消去法n顺序高斯消去法n列选主元高斯消去法n全选主元高斯消去法n三角分解法n杜利特尔、克劳特三角分解法n平方
27、根法、改进的平方根法n追赶法n进一步需要考虑的问题：上述方法求解方程组进一步需要考虑的问题：上述方法求解方程组可能存在的问题。我们从研究方程组性态的工可能存在的问题。我们从研究方程组性态的工具具-向量范数、矩阵范数开始向量范数、矩阵范数开始.2.3 向量和矩阵的范数向量和矩阵的范数n为了研究线性方程组近似解的误差估计和下一章要讲的迭代法的收敛性，我们需要对Rn(n维向量空间)中的向量或Rnxn中矩阵的“大小”引入一种度量：n向量和矩阵的范数。.向量和矩阵的范数n在一维数轴上，实轴上任意一点x到原点的距离用|x|表示。而任意两点x1，x2之间距离用|x1-x2|表示。.向量和矩阵的范数n而在二维
28、平面上，平面上任意一点P(x,y)到原点的距离用表示。而平面上任意两点P1(x1,y1),P2(x2,y2)的距离用表示。推广到n维空间，则称为向量范数。首先引入向量范数的定义.向量范数.常见的向量范数可以验证上述三种定义的向量空间的常用范数都满足向量范数的定义中的三个条件，均为向量范数。分别称为1-范数，2-范数，-范数。 .常见的向量范数参见教材61页，验证2-范数满足向量范数的三个条件首先是非负性，对x0,有2-范数0,对x=0,有2-范数=0；其次是齐次性，对任意c,有cx的2-范数等于|c|乘x的2-范数；然后是三角不等式，这里要用到Cauchy-Schwarz不等式，证明如下.常见
29、的向量范数因而，因而，2-范数是向量范数。范数是向量范数。.向量范数性质向量范数是等价的。向量范数是等价的。.向量范数性质等价性质：下面证明下面证明1）还有其它形式的向量范数之间的等价关系，。还有其它形式的向量范数之间的等价关系，。.向量范数性质等价性质1）的证明：.向量范数性质等价性质证明：.2.3.2矩阵范数.相容范数.相容范数.算子范数是矩阵范数下面证明算子范数是矩阵范数！下面证明算子范数是矩阵范数！.算子范数是矩阵范数.算子范数是矩阵范数.算子范数是矩阵范数.常见的矩阵范数.常见的矩阵范数.常见的矩阵范数.常见的矩阵范数.例题.2.3.3矩阵的谱半径和矩阵序列收敛性.谱半径和矩阵序列的
30、收敛性.例题.例题.2.4病态方程组与矩阵的条件数.2.4.1病态方程组与扰动方程组的误差分析.病态方程组与扰动方程组的误差分析.病态方程组与扰动方程组的误差分析.病态方程组与扰动方程组的误差分析.病态方程组与扰动方程组的误差分析.病态方程组n扰动方程由于计算机字长限制，在解AX=b时，舍入误差是不可避免的。因此我们只能得出方程的近似解。是方程组（A+A）x=b+b(1)这就是所谓的扰动方程。.病态方程组称（A+A）x=b+b为方程Ax=b的扰动方程。其中A，b为由舍入误差所产生的扰动矩阵和扰动向量。当A，b有微小扰动，解得扰动方程（1）的解与Ax=b的解x的相对误差不大时，称为良态方程，否则
31、为病态方程。 下面运用前面介绍的向量范数、矩阵范数来研下面运用前面介绍的向量范数、矩阵范数来研究方程组解得误差传播规律。究方程组解得误差传播规律。.扰动方程组的误差界-右端项扰动b.扰动方程组的误差界-系数矩阵扰动A.扰动方程组的误差界-系数矩阵扰动A、右端项扰动b.总结将上述三种情况，我们就能得到在系数矩阵有扰动A和/或右端项有扰动b时，方程组的解的误差限！在上述三种情况中都出现了一个决定误差界大小的量，即|A|A-1|。也就是这个数决定了方程组得性态，我们将它称之为条件数。.2.4.2矩阵的条件数.矩阵的条件数的性质.扰动方程组的误差界-定理.用条件数来刻画方程组的性态.例题，计算矩阵的条
32、件数这里，我们用了条件这里，我们用了条件数的性质来计算对称数的性质来计算对称矩阵的条件数，也可矩阵的条件数，也可以用条件数的定义来以用条件数的定义来计算。即先求计算。即先求A的逆的逆矩阵，再求矩阵，再求A的条件的条件数！数！.例：例：Hilbert 阵阵cond (H2) = 27cond (H3) 748cond (H6) = 2.9 106注：注：现在用现在用MatlabMatlab数学软件可以很方便数学软件可以很方便求矩阵的条件数求矩阵的条件数! !对于严重的病态线性方程组，即使原始数据对于严重的病态线性方程组，即使原始数据A和和b都没有都没有误差，但如果在求解过程中有舍入误差，所得到的
33、解也误差，但如果在求解过程中有舍入误差，所得到的解也会有很大的相对误差。会有很大的相对误差。.什么样的线性方程组可能是病态的？什么样的线性方程组可能是病态的？注：注：一般判断矩阵是否病态，并不计算一般判断矩阵是否病态，并不计算A 1，而由经验得出。，而由经验得出。 行列式很大或很小（如某些行、列行列式很大或很小（如某些行、列近似相关）；近似相关）； 元素间相差大数量级，且无规则；元素间相差大数量级，且无规则； 主元消去过程中出现小主元；主元消去过程中出现小主元； 特征值相差大数量级。特征值相差大数量级。.缓和甚至解决线性方程组病态的手段缓和甚至解决线性方程组病态的手段1.1.采用高精度的算法（
34、不是总有效）；采用高精度的算法（不是总有效）；采用高精度的算法（不是总有效）；采用高精度的算法（不是总有效）；2.2.通过行或列的平衡来降低矩阵的条件数（平衡法）；通过行或列的平衡来降低矩阵的条件数（平衡法）；通过行或列的平衡来降低矩阵的条件数（平衡法）；通过行或列的平衡来降低矩阵的条件数（平衡法）；3.3. 残差校正法；残差校正法；残差校正法；残差校正法；4.4.通过矩阵的预条件处理来降低矩阵的条件数。通过矩阵的预条件处理来降低矩阵的条件数。通过矩阵的预条件处理来降低矩阵的条件数。通过矩阵的预条件处理来降低矩阵的条件数。.残差校正法步骤：残差校正法步骤：1.求解求解Axb，得到，得到x的近似
35、解的近似解x1；2.校正校正x1，令，令r1bAx1，解，解A x1 r1；3.x2x1 x1 ；4.令令r2bAx2，解，解A x2 r2；5.x3x2 x2 ；6.。7.令令rkbAxk，解，解A xk rk；8.Xk1xk xk ；。；。.预条件技术n80年代提出，纯代数观点nM 近似 A, 求解(左预条件)nM y = c 易解n n相当于化学反应中寻找高效、廉价的催化高效、廉价的催化剂剂，能极大地提高迭代方法的速度。.预条件技术(续)n代数预条件技术代数预条件技术nILU、SPAI、SOR、多项式n几何预条件技术几何预条件技术n某方向压缩粗化、网格均匀化、区域规则化近似n分析预条件技
36、术分析预条件技术n变系数常数化、Green函数稀疏近似、强化椭圆型n物理预条件物理预条件n量纲平衡、物理参数逼近、状态方程近似、某些物理项简化、算子分裂n高级预条件高级预条件nMG、DDM、快速变换（FFT）、基底变换法.解线性代数方程组的迭代解法.In Scientific ComputingLarge Linear SystemsAx=bas sub-problems/as intermediate stepsGauss-Seidel methodJacobi methodSOR methodConjugate Gradient methodfor symmetric systemsGau
37、ssian eliminationLU factorizationCholesky factorizationGMRESGCRBi-CGCGSBi-CGSTABBi-CGSTAB2GPBi-CGBi-CGSTAB(L).173直接法和迭代法的比较直接法: 经过有限次运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷列去逼近精确解的方法。（一般有限步内得不到精确解） 直接法比较适用于中小型方程组。对高阶方程组，既使系数矩阵是稀疏的，但在运算中很难保持稀疏性，因而有存储量大，程序复杂等不足。 迭代法则能保持矩阵的稀疏性，具有计算简单，编制程序容易的优点，
38、并在许多情况下收敛较快。故能有效地解一些高阶方程组。.解线性代数方程组的迭代法解线性代数方程组的迭代法引进迭代法的目的引进迭代法的目的引进迭代法的目的引进迭代法的目的：1.1.解大型方程组；解大型方程组；解大型方程组；解大型方程组；2.2.存储的考虑；存储的考虑；存储的考虑；存储的考虑；3.3.解病态方程解病态方程解病态方程解病态方程 。.2.5 解线性代数方程组的迭代法解线性代数方程组的迭代法迭代发的一般理论迭代发的一般理论Jacobi迭代法、迭代法、G-S迭代法迭代法松弛迭代法松弛迭代法迭代法的收敛条件迭代法的收敛条件迭代法收敛的其它判定方法迭代法收敛的其它判定方法解线性方程组的迭代法解线
39、性代数方程组的极小化方法解线性代数方程组的极小化方法.2.5 解线性代数方程组的迭代法解线性代数方程组的迭代法迭代法的一般理论迭代法的一般理论Jacobi迭代法、迭代法、G-S迭代法迭代法松弛迭代法松弛迭代法解线性方程组的迭代法 2.5.1 一般迭代法一般迭代法.引例引例解解迭代法的基本思想.或写为或写为 ，其中其中迭代法的基本思想.迭代法的基本思想.迭代法的基本思想.我们任取初始值我们任取初始值，例如例如 ，将这些值代入将这些值代入（2）式右边式右边（若若（2）式为等式即式为等式即求得方程组的解求得方程组的解，但一般不满足但一般不满足），），得到新的值得到新的值得到向量序列和一般的计算公式（
40、迭代公式）得到向量序列和一般的计算公式（迭代公式）迭代法的基本思想.其中其中k k表示迭代次数表示迭代次数（k k=0,1,2, ）.迭代到第迭代到第1010次有次有迭代法的基本思想.即建立一种从已知近似解到新的近似解的计算公式。即建立一种从已知近似解到新的近似解的计算公式。如何建立这样的迭代计算公式呢？如何建立这样的迭代计算公式呢？迭代法的基本思想迭代法求方程组的解的基本思想即是：构造迭代法求方程组的解的基本思想即是：构造一串收敛到解的序列，一串收敛到解的序列，.思路思路基本迭代法.基本迭代法.将线性方程组将线性方程组 A Ax x=b =b 化为化为 x x=M=Mx x+g+g再由此再由
41、此构造向量序列构造向量序列 x x (k)(k): x x(k+1)(k+1)=M=Mx x (k)(k)+g+g若若 x x (k)(k) 收敛至某个向量收敛至某个向量x x *,*,则可得向量则可得向量x x * *就是就是所求方程组所求方程组 A Ax x=b =b 的准确解的准确解. .迭代法的基本思想我们已经回答了我们已经回答了 “如何建立这样的迭代计算公式如何建立这样的迭代计算公式呢？呢？”的问题。方法是：的问题。方法是：下面就来构造我们这一章中将要介绍的第一种迭代法下面就来构造我们这一章中将要介绍的第一种迭代法.Jacobi迭代法.Jacobi迭代法.Jacobi迭代法.Jaco
42、bi迭代法.求解一般线性代数方程组的求解一般线性代数方程组的JacobiJacobi迭迭代法的计算公式代法的计算公式Jacobi迭代法.举 例按照按照jacobi法的一般计算公式写出本例的法的一般计算公式写出本例的jacobi迭代式迭代式.JacobiJacobi迭代法的特点迭代法的特点Jacobi迭代法的特点.Jacobi迭代法.Gauss-Seidel迭代法.Gauss-Seidel迭代法.Gauss-Seidel迭代法.即即：Gauss-Seidel迭代法.用高斯用高斯- -赛德尔迭代法解方程组赛德尔迭代法解方程组例例4.1Gauss-Seidel迭代法举例.如此继续下去如此继续下去Ga
43、uss-Seidel迭代法的特点. 简单，每迭代一次只需进行一次矩阵和向量乘法，简单，每迭代一次只需进行一次矩阵和向量乘法，简单，每迭代一次只需进行一次矩阵和向量乘法，简单，每迭代一次只需进行一次矩阵和向量乘法， 不改变系数矩阵的稀疏性，不改变系数矩阵的稀疏性，不改变系数矩阵的稀疏性，不改变系数矩阵的稀疏性， 仅需要一组存储单元仅需要一组存储单元仅需要一组存储单元仅需要一组存储单元 。 Gauss-SeidelGauss-Seidel迭代法一般比迭代法一般比迭代法一般比迭代法一般比JacobiJacobi迭代法快，迭代法快，迭代法快，迭代法快，但也有比但也有比但也有比但也有比 JacobiJa
44、cobi迭代法慢的，甚至有迭代法慢的，甚至有迭代法慢的，甚至有迭代法慢的，甚至有JacobiJacobi迭代法收敛而迭代法收敛而迭代法收敛而迭代法收敛而Gauss- Gauss- Seidel Seidel迭代法不收敛的迭代法不收敛的迭代法不收敛的迭代法不收敛的。Gauss-Seidel迭代法特点迭代法特点.举 例请按照下述请按照下述Gauss-Seidel法的一般计算公式写出该法的一般计算公式写出该例的例的Gauss-Seidel迭代式迭代式.上面我们由介绍了一般迭代法的构造方法，又上面我们由介绍了一般迭代法的构造方法，又由一般迭代方法具体引出了由一般迭代方法具体引出了Jacobi迭代法，进
45、迭代法，进一步考虑迭代过程中新信息的利用，得到了一步考虑迭代过程中新信息的利用，得到了Gauss-seidel迭代法。迭代法。当然，这两种特殊的迭代法同样可以直接由系当然，这两种特殊的迭代法同样可以直接由系数矩阵的和、差分解导出。数矩阵的和、差分解导出。下面我们在分析这两者方法存在误差的基础上，下面我们在分析这两者方法存在误差的基础上，提出一种新的方法，即所谓松弛迭代法。首先提出一种新的方法，即所谓松弛迭代法。首先来看来看松弛迭代法的基本思想。松弛迭代法的基本思想。松 弛 法.松弛法的基本思想我们从为我们从为Gauss-Seidel 迭代法加速的角度来考虑迭代法加速的角度来考虑提出松弛法。提出
46、松弛法。Gauss-Seidel 迭代公式得到，于是有迭代公式得到，于是有这样我们就可以得到这两次近似解的差别这样我们就可以得到这两次近似解的差别x,x,即即.松弛法的基本思想具体地具体地.可以把可以把 看作看作Gauss-Seidel 迭代的修正项，即第迭代的修正项，即第k次次近似解近似解 以此项修正后得到新的近似解。以此项修正后得到新的近似解。我们有理由相信，新的近似解，应有更好的精度，误我们有理由相信，新的近似解，应有更好的精度，误差更小。当然，考虑到差更小。当然，考虑到xx也是利用两次近似估计得来也是利用两次近似估计得来的，并不是真正的第的，并不是真正的第k k次近似解的误差，所以在实
47、际修次近似解的误差，所以在实际修正时，更好的做法是对修正项正时，更好的做法是对修正项xx加权！具体地有：加权！具体地有：松弛法的基本思想具体地，我们具体地，我们.得到新的近似解，具体公式为得到新的近似解，具体公式为松弛法松弛法是将是将 乘上一个参数因子乘上一个参数因子 作为修正项而作为修正项而松弛法的基本思想所以，松弛法就可以进一步表示为：所以，松弛法就可以进一步表示为：.称为称为松弛因子松弛因子 当当 时称为时称为低松弛低松弛；时称为时称为超松弛法超松弛法。是是Gauss-Seidel迭代迭代松弛法的基本思想上面我们运用误差估计上面我们运用误差估计校正法，得到了松弛法的分校正法，得到了松弛法
48、的分量形式，这已经可以方便地让我们编程计算了，但分量形式，这已经可以方便地让我们编程计算了，但分量形式不利于分析推理，下面将松弛法表为矩阵形式：量形式不利于分析推理，下面将松弛法表为矩阵形式：.选取分裂矩阵选取分裂矩阵M M为带参数下三角阵为带参数下三角阵其中其中 为可选择的松弛因子为可选择的松弛因子.松弛法的矩阵表示.松弛法的矩阵表示通过简单的推导，得到通过简单的推导，得到.松弛法的矩阵表示对照前面给出的分量形式！对照前面给出的分量形式！.SORSOR方法的方法的计算公式计算公式SOR 方 法对照前面给出的分量形式！对照前面给出的分量形式！.例例3取取 用松弛法解方程组用松弛法解方程组举 例
49、解：由解：由.1. 1. 1.1. 1. 1.561. 1.392 1.5321.2744 1.3528 1.6021.1372 1.40376 1.591641.22775 1.39396 1.601081.18467 1.40106 1.598891.20687 1.39917 1.600241.19667 1.40021 1.599841.20148 1.39988 1.600041.19932 1.40004 1.599981.2003 1.39998 1.600011.19987 1.40001 1.61.20006 1.4 1.61.19998 1.4 1.61.20001 1.4
50、 1.61.2 1.4 1.6计算结果计算结果松弛法的程序.验证验证： 松弛因子与迭代次数的关系松弛因子与迭代次数的关系例例4举 例.举 例.松弛因子松弛因子 满足误差满足误差松弛松弛因子因子 满足误差满足误差1.0221.5171.1171.6231.2121.7331.3111.8531.4141.9109举 例.举 例请按照下述松弛法的一般计算公式写出该例的请按照下述松弛法的一般计算公式写出该例的松弛迭代式松弛迭代式. 2.5.2 迭代法的收敛条件迭代法的收敛条件1、向量序列收敛与矩阵序列收敛、向量序列收敛与矩阵序列收敛2、迭代法收敛概念、迭代法收敛概念3、迭代法收敛定理、迭代法收敛定理
51、4、举例、举例迭代法的收敛条件.思思考考不一定不一定!定义定义4.1关注的问题及几个相关概念.定理定理关注的问题及几个相关概念.证明证明：上述定理表明，向量序列和矩阵序列的收敛可上述定理表明，向量序列和矩阵序列的收敛可归结为对应分量或对应元素序列的收敛。归结为对应分量或对应元素序列的收敛。矩阵有类矩阵有类似结果似结果关注的问题及几个相关概念.定义定义4.2迭代矩阵.定义定义4.3矩阵的谱半径回顾我们上一章中学习过的矩阵的谱半径的定义回顾我们上一章中学习过的矩阵的谱半径的定义谱半径是我们讨论迭代法收敛性的重要工具！它谱半径是我们讨论迭代法收敛性的重要工具！它具有以下性质：具有以下性质：.矩阵的谱
52、半径.零矩阵的充要条件有了谱半径，就可以给出前面有了谱半径，就可以给出前面定义的向量序列、矩阵序列收定义的向量序列、矩阵序列收敛的条件了！敛的条件了！当然也就可以讨论迭代法的收当然也就可以讨论迭代法的收敛性问题了。敛性问题了。.定理定理4.14.1证明：（必要性）证明：（必要性）零矩阵的充要条件.零矩阵的充要条件有了上述向量序列、矩阵序列收有了上述向量序列、矩阵序列收敛的条件，也就等价地有了敛的条件，也就等价地有了-一般迭代法收敛条件了。我们有一般迭代法收敛条件了。我们有迭代法的如下基本定理：迭代法的如下基本定理：.定理定理4.2证明证明：迭代法的收敛定理.零矩阵的充要条件有了有了一般迭代法收
53、敛的充要条件，从理论一般迭代法收敛的充要条件，从理论上解决了迭代法收敛性的判定问题。但收上解决了迭代法收敛性的判定问题。但收敛的速度、误差传播的规律或者误差估计敛的速度、误差传播的规律或者误差估计的问题仍然没有解决。的问题仍然没有解决。还记得我们有还记得我们有 的结果吗？如果我的结果吗？如果我们能将迭代矩阵的条件加强为们能将迭代矩阵的条件加强为 ，我，我们应该有希望得到我们想要的误差估计：们应该有希望得到我们想要的误差估计：.和和误 差 估 计.误 差 估 计.有有由此可得由此可得误 差 估 计.所以所以误 差 估 计由于由于.取范数取范数移项得移项得代入得代入得误 差 估 计证毕证毕.误差估
54、计式误差估计式 给出了停止迭代的一个判别准则，即用相邻两给出了停止迭代的一个判别准则，即用相邻两次迭代结果的差的范数来判别，它越小，迭代序列与次迭代结果的差的范数来判别，它越小，迭代序列与精确解越接近。精确解越接近。由误差估计式由误差估计式根据事先给定的精度根据事先给定的精度可估计出迭代的次数可估计出迭代的次数误 差 估 计.例：方程组例：方程组解：用雅可比迭代法解：用雅可比迭代法误 差 估 计 举 例.即需要迭代即需要迭代13次才能满足精度。次才能满足精度。则有则有误 差 估 计 举 例.迭代法收敛性的判定迭代法收敛性的判定2 2、若、若A A为严格对角占优阵或不可约弱对角占优阵，则为严格对
55、角占优阵或不可约弱对角占优阵，则JacobiJacobi迭代法和迭代法和Gauss-SeidelGauss-Seidel迭代法均收敛。迭代法均收敛。小 结.三种迭代法收敛性判定的其它定理三种迭代法收敛性判定的其它定理n前面，我们已经给出了一个一般迭代法前面，我们已经给出了一个一般迭代法收敛判定的充要条件和一个迭代矩阵的收敛判定的充要条件和一个迭代矩阵的某种范数小于某种范数小于1的特定迭代法收敛的充分的特定迭代法收敛的充分条件。下面，我们针对一些更特别的方条件。下面，我们针对一些更特别的方程组的迭代法来给出几个收敛判定定理。程组的迭代法来给出几个收敛判定定理。.若若n阶方阵阶方阵 满足满足定义定
56、义1 1定义定义2 2如果矩阵如果矩阵A不能通过行的互换和相应的列互换成不能通过行的互换和相应的列互换成为形式为形式 ，其中，其中 为方阵，则称为方阵，则称A为为不可约。不可约。且至少有一个且至少有一个i值值,使得上式中不等号严格成立，则称使得上式中不等号严格成立，则称A为为弱对角占优阵弱对角占优阵。若对所有若对所有i,上式不等号均严格成立上式不等号均严格成立 ，则称则称A为为严格对角占优阵严格对角占优阵。定 义例例.定理定理定 理设有线性方程组设有线性方程组 ，下列结论成立，下列结论成立1. 若若A为严格对角占优阵或不可约弱对角占优阵，则为严格对角占优阵或不可约弱对角占优阵，则Jacobi
57、迭代法和迭代法和Gass-Seidel迭代法均收敛。迭代法均收敛。2.若若A为严格对角占优阵，为严格对角占优阵， 则松弛法收敛。则松弛法收敛。3.若若A为对称正定阵，则松弛法收敛的充要条件为为对称正定阵，则松弛法收敛的充要条件为上两例中：上两例中：A为严格对角占优阵，为严格对角占优阵， Jacobi 迭代法和迭代法和Gass-Seidel迭代均收敛。迭代均收敛。B为非严格对角占优阵，为非严格对角占优阵， 故松弛法收敛。故松弛法收敛。.例例5 5解：因解：因A为对称且各阶主子式皆大于零，故为对称且各阶主子式皆大于零，故A为对称正定为对称正定均收敛。均收敛。A不是弱对角占优阵，故不能用条件不是弱对
58、角占优阵，故不能用条件1判断。判断。阵。由判别条件阵。由判别条件3，Gauss-Seidel迭代法与松弛法迭代法与松弛法Jacobi迭代法迭代矩阵迭代法迭代矩阵举 例.其特征方程其特征方程得得因而因而迭代法不收敛迭代法不收敛举 例.Jacobi与与Gauss-Seidel迭代的迭代矩阵分别为迭代的迭代矩阵分别为改变方程组中方程的次序，即将系数矩阵作行改变方程组中方程的次序，即将系数矩阵作行交换，会改变迭代法的收敛性交换，会改变迭代法的收敛性谱半径分别是谱半径分别是均不收敛均不收敛例例举 例.若交换方程的次序，得若交换方程的次序，得 Ax=b 的同解方程组的同解方程组为严格对角占优矩阵，因而对方
59、程组为严格对角占优矩阵，因而对方程组 用用 Jacobi 与与Gass-Seidel 迭代求解均收敛。迭代求解均收敛。举 例.对称超松弛迭代法（对称超松弛迭代法（SSOR法）法）.对称超松弛迭代法（对称超松弛迭代法（SSOR法）法）.SSORSSOR迭代法的矩阵形式：迭代法的矩阵形式：注：注：（1）关于）关于 的收敛条件和准则与的收敛条件和准则与SOR方法相同；方法相同；（2）收敛快慢对）收敛快慢对 的选取不敏感。的选取不敏感。.块超松弛迭代法（块超松弛迭代法（BSORBSOR法）法）.案案 例例 1 (1 (热传导问题热传导问题) )设有一维热传导方程的初边值问题设有一维热传导方程的初边值问
60、题试用数值方法求出试用数值方法求出t=0.2时刻金属杆的温度分布时刻金属杆的温度分布.应用实例应用实例.解解:(1)对空间进行离散对空间进行离散.(2)对微分算子进行离散对微分算子进行离散.t0.020.0100 0.1 0.2 . 0.9 1 x.采用无条件稳定的采用无条件稳定的Crank-Nicholson格式格式,则有则有.或或加上边界条件后有加上边界条件后有.加上边界条件后有加上边界条件后有其矩阵形式为其矩阵形式为 ATn+1=dn.案案 例例 2 2 数值求解正方形域上的数值求解正方形域上的Poisson方程边值问题方程边值问题.解解:(1)剖分求解域剖分求解域.(2)对微分算子进行
61、离散对微分算子进行离散.0 1 2 . N N+1 XYN+1N：210.在每个点在每个点(xi,yj)上的有限差分方程为上的有限差分方程为在边界上在边界上又称为五点差分格式又称为五点差分格式(i,j)(i,j+1)(i,j-1)(i-1,j)(i+1,j).对非边界点进行编号对非边界点进行编号: 顺序为顺序为-从下往上从下往上,从左往右从左往右相应的解向量和右端向量分别为相应的解向量和右端向量分别为 .差分方程组的矩阵形式为差分方程组的矩阵形式为 Au=f如果把每一条线上的网点看作一个组如果把每一条线上的网点看作一个组,如如.其中其中可用块迭代法可用块迭代法求解求解Au=f.2.5.32.5
62、.3解线性代数方程组的极小化方解线性代数方程组的极小化方法法一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题三、共轭斜量法三、共轭斜量法(共轭梯度法共轭梯度法)四、预条件共轭斜量法四、预条件共轭斜量法二、最速下降法二、最速下降法.n共轭梯度法最早是由计算数学家Hestenes和几何学家Stiefel在20世纪50年代初为求解线性方程组而各自独立提出的。他们合作的文章被公认为共轭梯度法的奠基之作，详细讨论了求解线性方程组的共轭梯度法的性质以及它和其他方法的关系。n在A为对称正定阵时，上述线性方程组等价于最优化问题。解线性代数方程组的极小化方法解线性代数方程组的极小化方法.n由此，He
63、stenes和Stiefel的共轭梯度方法也可视为求二次函数极小值的共轭梯度法。n1964年Fletcher和Reeves将此方法推广到非线性优化，得到求一般函数极小值的共轭梯度法。解线性代数方程组的极小化方法解线性代数方程组的极小化方法.n提到最优化问题，首先介绍何为最速下降法。n考虑线性方程组Ax=b的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量。n为此我们定义二次泛函解线性代数方程组的极小化方法解线性代数方程组的极小化方法.n定理：设A对称正定，求方程组Ax=b的解，等价于求二次泛函的极小点n由定理，求解线性方程组的问题就转化为求二次泛函的极小点的问题，也可表为与线性方程
64、组等价的变分问题：解线性代数方程组的极小化方法解线性代数方程组的极小化方法.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问题一、与线性方程组等价的变分问题.一、与线性方程组等价的变分问
65、题一、与线性方程组等价的变分问题.二、最速下降法二、最速下降法.二、最速下降法二、最速下降法.二、最速下降法二、最速下降法.二、最速下降法二、最速下降法.三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG)
66、 (共轭梯度法共轭梯度法) .公式化简公式化简.三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .公式化简公式化简三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .三、共轭斜量法三、共轭斜量法 (CG) (共轭梯度法共轭梯度法) .四、预条件共轭斜量法四、预条件共轭斜量法(PCG).预条件共轭斜量法预条件共轭斜量法 实际计算，可通过变实际计算，可通过变换，转化成用原方程组换，转化成用原方程组的量来计算。的量来计算。.预优矩阵
67、的选取预优矩阵的选取下面介绍几种选取预优矩阵的方案下面介绍几种选取预优矩阵的方案:.几种选取预优矩阵的方案几种选取预优矩阵的方案.几种选取预优矩阵的方案几种选取预优矩阵的方案.几种选取预优矩阵的方案几种选取预优矩阵的方案.三种预条件共轭斜量法在三种预条件共轭斜量法在MATLAB中的调用格式中的调用格式:1.不用预优矩阵的共轭斜量法不用预优矩阵的共轭斜量法 x=pcg(a,b,tol,kmax)2.用预优矩阵的共轭斜量法用预优矩阵的共轭斜量法 (1)x=pcg(a,b,tol,kmax,m) (2) r=chol(m) x=pcg(a,b,tol,kmax,r,r,x0)3.未给定预优矩阵的共轭
68、斜量法未给定预优矩阵的共轭斜量法 r=cholinc(sa,0) x=pcg(a,b,tol,kmax,r,r,x0).nCG(ConjugateGradient)（对称正定）nMRES(MinimalRESidual)(对称不正定）nGMRES(GenerilizedMinimalRESidual）（不对称不正定）nBiCG(不对称不正定）nBiCGSTAB(不对称不正定）几种常用的几种常用的Krylov-subspace方法方法.背景n大规模线性方程组的求解很多科学工程计算问题都转化为求解方程组Ax=b.如偏微分方程组的差分格式，有限元方法离散得到刚度矩阵.n大规模矩阵特征值和特征向量的计
69、算工程计算领域十分常见。如量子物理中的Kohn-Sham方程求解化为哈密顿矩阵某些关键特征值对的计算.Krylov-subspace方法.投影方法n线性方程组的投影方法方程组Ax=b,A是nn的矩阵.给定初始x(0),在m维空间K(右子空间）中寻找x的近似解x(1)满足残向量r=b-Ax(1)与m维空间L(左子空间）正交，即b-Ax(1)L此条件称为Petrov-Galerkin条件.当空间K=L时，称相应的投影法为正交投影法，否则称为斜交投影法.K(L)X(0)X(1)KLX(0)X(1).n解方程组的投影法的矩阵表示设nm阶矩阵V=v(1),v(2),v(m)与W=w(1),w(2),w(
70、m)的列分别构成K与L的一组基.记z=x(1)-x(0),z=Vy,有当非奇异时有从而得到迭代公式.n投影方法的最优性.(误差投影)设A为对称正定矩阵,x(0)为初始近似解,且K=L,则x(1)为采用投影方法得到的新近似解的充要条件是其中,且x为问题的精确解.(残量投影)设A为任意方阵, x(0)为初始近似解,且L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是其中.n矩阵特征值的投影方法对于特征值问题Ax=x,其中A是nn的矩阵，斜交投影法是在m维右子空间K中寻找xi和复数i满足Axi-ixiL,其中L为m维左子空间.当L=K时，称此投影方法为正交投影法.Kyrlov子空间定义为m
71、维Krylov子空间为v的次数，即使得q(A)v的非零首一多项式的最低次数.Krylov子空间的标准正交基nArnoldi方法基于Gram-Schmit正交化方法首先，选取一个Euclid范数为1的向量v(1),对,通常可取在已知v(1),v(2),v(j)的情况下,不妨设v(1),v(2),v(j),Av(j)线性无关(否则构造完毕),则可求出与每个都正交的向量.而不难看出，再记,得到与v(1),v(2),v(j)都正交的向量重复此过程，即可得到一组标准正交基.若期间某个j使得hj+1,j=0,则说明v的次数是j,且Kj是A的不变子空间.n.定理定理 如果记以v(1),v(2),v(m)为列
72、构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为删除最后一行得到的矩阵为Hm,则在Arnoldi算法中，可能有较大舍入误差，改写.Krylov子空间方法解线性方程组n误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法.不难看出，当采用上述FOM算法时，需要存储所有的vi，(i=1,2,m)，当m增大时，存储量以O(mn)量级增大.而FOM计算量是O(m2n).可见其代价十分高昂.因此我们考虑重启的FOM算法.Krylov子空间方法解线性方程组n误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法2)非对称矩阵的IOM方法和DIOM方法.所谓不完全正交
73、化方法(IOM),是指在正交化过程中,v(j+1)仅与最近k个v(j-k+1),v(j-1),v(j)正交,这样做虽然破坏的正交性,但是降低了计算量.当然k选得越小，对每个j对应的计算量也越小,但可能要选更大的m才能取得满足精度要求的近似解.IOM算法仅仅是把FOM算法的第三步改为(iii)对i=max(1,j-k+1),j,计算hij与w(j).但采用IOM后，仍然需要存储v(1),v(2),v(m),因为在第(vi)步中仍然需要这些向量.解决这个问题可以考虑采用H的LU分解，通过自身分解的迭代更新以减少每一步的存储量.Krylov子空间方法解线性方程组n误差投影型方法取L=K时的正交投影法
74、）非对称矩阵的FOM方法2)非对称矩阵的IOM方法和DIOM方法）对称矩阵Lanczos算法.A是非对称阵时,H是Hessenberg阵，则当A是对称阵时，H也为对称阵，故应为三对角阵，记为对它采用修正Arnoldi正交化过程，就得到Lanczos方法.Krylov子空间方法解线性方程组n误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法2)非对称矩阵的IOM方法和DIOM方法）对称矩阵Lanczos算法4）正定对称阵的CG法.CG法是解正定对称线性方程组最有效的方法之一该方法充分利用了矩阵A的稀疏性，每次迭代的主要计算量是向量乘法而实际上，CG方法也是一种Krylov子空间方法后面
75、将给出一个数值例子.Krylov子空间方法解线性方程组n残量投影型方法取L=AK时的斜交投影法）GMERS方法.GMERS方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法.回顾.(残量投影)设A为任意方阵, x(0)为初始近似解,且L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是其中.GMERS方法在Arnoldi正交化过程之后把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法不去直接求x(1),转而去解此极小问题，这是GMERS方法的独到之处.再由之前的定理.Krylov子空间方法解线性方程组n残量投影型方法取L=AK时的斜交投影法）GMERS方法2)重
76、启型GMERS方法、QGMERS、DGMERS.GMERS算法具有很好的收敛性，在不考虑舍入误差的情况下，该算法能在在有限步得到精确解.但是，和FOM算法类似，它需要保存大量正交基，因此我们可以采取类似重启型FOM的算法实现重启型GMERS算法.同样可以采取不完全正交的方法，在正交化过程中,v(j+1)仅与最近k个v(j-k+1),v(j-1),v(j)正交就能得到QGMERS方法.为了避免存储v(1),v(2)v(m-k)，可以对进行QR分解.这种算法被称为DQGMERS.Krylov子空间方法解矩阵特征值问题nArnoldi方法求解非对称矩阵特征值nLanczos方法求解对称矩阵特征值.n
77、Arnoldi方法再次回顾前面的定理而所以有以下结论：.1）若,则是A的不变子空间,从而Hm的所有特征值是A的特征值子集2）若,设Hm的特征值对是,即,由上述定理可得.以上讨论说明,可以用Hm的特征值(又称Ritz值)来近似A的特征值,用向量(又称Ritz向量)来近似相应的特征向量.事实上,Hm为A在Kyrlov子空间上的正交投影.由于Hm是上Hessenberg阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解.Arnoldi算法构造标准正交基1需要存储所有的基向量,当m很大时,存储量大2理论上为了保证收敛速度,m越大越好矛盾！.nLanczos方法A是对称阵时,Hm是三对角阵，仍
78、然采用Hm的特征值来近似A的特征值，相应的Ritz向量来近似A的特征向量.问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，QR法等，可以在并行机上达到O(logN)的量级.对称Lanczos算法,在没有舍入误差的情形下,得到的是标准正交基相互正交的,并且至多n步必然终止.但是在出现舍入误差的情况下，计算得到的将失去正交性.因此,长期以来被人们认为是不稳定的.直到1971年,Paige通过仔细的舍入误差分析,发现失去正交性恰与近似特征值的精度提高有关.之后,经过大量的理论分析和数值实验,人们充分认识到Lanczos方法是求解大型稀疏对称矩阵特征值问
79、题的最有效方法之一.目前，Lanczos方法是求大型稀疏对称矩阵部分极端特征值的最有效方法.改进技术n预条件法Krylov子空间方法的收敛性一般都和矩阵的特征值分布有关，特征值分布越集中，收敛速度越快，若特征值分布较分散，则收敛的很慢，此时可以采用预条件子来改善矩阵的特征值分布.选取矩阵使得A的特征值比A集中，并且M-1容易求得，于是将原问题化为问题M-1Ax=M-1b,这是左预条件法，还可采用右预条件法和分裂预条件法.nLanczos双正交化方法对于非对称线性方程组，也可采用类似于Lanczos算法将非对称阵化为三对角阵，它的好处在于可以可以形成短迭代，而FOM、IOM、GMERS的计算量和
80、存储量都随着m的增大而增大在双正交化过程中，取容易看出A和在其中扮演对偶的角色，此方法特别适用于需要求解两个系数矩阵分别是A和AT的方程组.n多项式加速技术在求解矩阵特征值问题时，可以利用chebeshev多项式来加快收敛所谓多项式加速,就是采用作为迭代初始向量,其中为多项式,如果,其中pi为A的特征值对应的特征向量则将特征值按实部递减顺序排列,其中为r个“想要”的特征值,将“不想要”的特征值点集用S表示.如果多项式,能使得S尽可能小,则相当于让求解过程产生加速.而Chebyshev多项式由于它的特殊性质,恰好能满足这种要求.不过，至今,仍无有效方法确定Chebyshev多项式的某些参数.n块
81、方法、调和投影法当矩阵的特征值分布较密集或有重特征值时，常规的Arnoldi方法或Lanczos方法就必须取m很大时才能收敛，此时可以采用块方法取可以使得无须很大的m就可让迭代收敛，并且此方法适合改为并行算法。当特征值密集时，还可采用调和投影法，这也是最近研究的一个热点。当左右子空间的基满足Wm=AVm时，称为调和投影，它除了可以通过选取位移点改善原矩阵的特征值分布外，还可以将内部特征值问题化为外部特征值问题.尚待解决的问题n在Arnodli过程和Lanczos双正交化过程中，会出现崩溃，怎样的矩阵容易出现崩溃。目前已有一些处理此崩溃的办法如前瞻技术，但并不十分成熟.n在投影类方法中Krylo
82、v方法是否是最优的,K是否有更好的选择.n怎样将预条件技术运用到求解矩阵特征值问题.解线性方程组迭代法总结n对任意的向量对任意的向量b，迭代法，迭代法 xk1Gxkb 收敛的充分必要条收敛的充分必要条件是件是(G)1n如果如果|G|1，则迭代法，则迭代法 xk1Gxkb 收敛收敛n经典的三种迭代法：经典的三种迭代法：Jacobi迭代法，迭代法，Gauss-Seidel迭代法迭代法和和SOR迭代法迭代法 nSOR迭代法中，松驰因子迭代法中，松驰因子 02n对SOR迭代法，目前迭代法，目前还没有一般通用的方法来没有一般通用的方法来计算最算最优的的松驰因子松驰因子n从并行能力来讲：从并行能力来讲：J
83、acobi迭代法优于迭代法优于Gauss-Seidel迭代法迭代法和和SOR迭代法迭代法 n实际应用中共轭梯度法及其变体（实际应用中共轭梯度法及其变体（KS方法）比经典的三种方法）比经典的三种迭代法优越迭代法优越 .In Scientific ComputingLarge Linear SystemsAx=bas sub-problems/as intermediate stepsGauss-Seidel methodJacobi methodSOR methodConjugate Gradient methodfor symmetric systemsGaussian eliminationLU factorizationCholesky factorizationGMRESGCRBi-CGCGSBi-CGSTABBi-CGSTAB2GPBi-CGBi-CGSTAB(L).
