资源描述
.,第 2 章解线性代数方程组的直接法与迭代法,.,引言,在自然科学和工程技术中，很多问题的解决常常归结为解线性代数方程组。例如: 船体数学放样中建立三次样条函数问题， 用最小二乘法求实验数据的曲线拟合问题， 非线性方程组求解的问题， 用差分法或者有限元法解常微分方程、偏微分方程边值问题等 都导致求解线性方程组，且常常为求解大型线性方程组。,.,引言,与上述这些问题有关的计算方法包括： 求解线性方程组的数值方法与计算矩阵的特征值及特征向量的数值方法。 在本章中，我们将注意力集中在线性方程组求解的直接方法和迭代方法。 也就是我们要学习的求解线性方程组的两类方法。,.,引言,线性方程组的这两类数值解法有以下特点。直接法：经过有限步算术运算，可求得方程组的精确解的方法（若在计算过程中没有舍入误差）；迭代法：用某种极限过程去逐步逼近线性方程组精确解的方法；迭代法具有占用存储单元少，程序设计简单，原始系数矩阵在迭代过程中不变等优点，但存在收敛性及收敛速度等问题。下面由我们从小熟知的消元法开始,.,2.1 高斯消元法,设线性方程组,简记为 AX=b（矩阵形式/向量形式）,.,高斯消元法,如果AX=b是n阶线性代数方程组，则,.,高斯消元法,克莱姆法则在理论上有着重大意义，但在实际应用中存在很大的困难，在线性代数中，为解决这一困难给出了高斯消元法。我们解释一下为何存在很大的困难,.,看来我们确实有必要学习方程组的求解方法，从一个简单的例子开始,Gauss消去法的计算量这是我们介绍完Gauss消去法之后，对其计算量的分析，大家先看看高斯法解n阶线性方程组需要的总乘除法次数为注1：当n = 20时约为2670次，比克莱姆法则9.71020次大大减少。注2：世界上最快的超级计算机每秒千万亿次的浮点运算能力,即：1015。注3：目前最好的台式或笔记本电脑每秒万亿次的浮点运算能力，即：1012 。请同学们估算一下，用克莱姆法则解20阶方程组各需要耗时多少？！ 超级计算机9.71020 / 1015 /3600 =2.8小时台式或笔记本电脑9.71020 / 1012 /3600 =2800小时,另外，还有存储空间相关的空间复杂性问题,.,例题,例1.用消元法解方程组,.,例题,第一步：-2 x（1）+(3)得,.,例题,第二步：1 x（2）+（4）,解方程组（1）（2）（5）得：x=1,2,3T,.,这样，一般方程组的求解问题就转化为一个上三角形方程组的求解问题。当然，如果你愿意，同样也可以转化为一个下三角形方程组的求解问题。,总结,下面我们就来看看三角形方程组的求解方法！,.,2.1.1 高斯顺序消元法,对下三角形方程的求解 设 （1）,.,高斯顺序消元法,由（1）得,.,对上三角方程组设,.,由（2）式回代得,.,一般线性代数方程组解法的思考,通过例子，我们得出了一般线性代数方程组的求解可以转化为上或下三角形方程组来求解的结论，上面我们又研究了上或下三角形方程组的求解方法，下面是讨论一般线性代数方程组求解的高斯消去法的时候了！。,.,高斯顺序消去法,对一般线性代数方程组设 Ax=b.,记A(1)=A b(1)=b,1、第一次消元。设,.,高斯顺序消去法,.,高斯顺序消去法,设第k-1次消元得A(k)x=b(k) 其中,下面进行第k次消元,.,高斯顺序消去法,第k次消元,.,高斯顺序消去法,最后,在完成n-1次消元后得到如下上三角形矩阵,这样，一般方程组就变成了所谓的三角形方程组，问题也就转化为三角形方程组的求解问题。,.,高斯顺序消去法,上述过程，也就是对于方程组AX=b系数矩阵做如下变换：,得到上三角方程组,.,高斯顺序消去法,.,高斯顺序消去法,.,高斯顺序消去法,注意：算法对系数矩阵对角元素的非零要求,.,高斯顺序消去法,.,高斯顺序消去法算法框图,高斯顺序消去法框图,.,高斯消去法的计算量,.,Gauss消去法的计算量以乘除法的次数为主(1) 消元过程：第k步时(n k) + (n k) (n k + 1) = (n k) (n k + 2)共有,.,Gauss消去法的计算量(2) 回代过程：求xi中，乘ni次，除1次，共ni+1次（i = 1，n1）共有,.,Gauss消去法的计算量(3) 高斯法解n阶线性方程组需要的总乘除法次数为注1：当n = 20时约为2670次，比克莱姆法则9.71020次大大减少。注2：世界上最快的超级计算机每秒万万亿次的浮点运算能力,即：1016。注3：目前最好的台式或笔记本电脑每秒万亿次的浮点运算能力，即：1012 。请同学们估算一下，用克莱姆法则解20阶方程组各需要耗时多少？！,.,高斯顺序消去法条件,.,进一步的考虑解决办法(1) 若消元过程中出现akk(k) = 0，则无法继续(2) 若akk(k) 0，但较小，则小主元做除数将产生大误差(3) 每次消元都选择绝对值最大者作主元，称为高斯主元消去法(4) 通常第k步选择第k列主对角元以下元素绝对值最大者作主元（该行与第k行对调），称为列主元消去法。,.,2.1.2 高斯主元素消去法,Gauss列主元消元法从第一列中选出绝对值最大的元素，如果 ，则交换第一行和第i行，即,交换,.,高斯列主元消去法算法,.,高斯列主元消去法,第k步 从 的第k列 ， ， 中选取绝对值最大项，记录所在行，即 若 交换第k行与l行的所有对应元素，再进行顺序消元。 下面，详细描述一下列选主元的高斯消去法,.,高斯列主元消去法,.,高斯列主元消去法框图,.,注意到算法选主元素的第2步，是非正常结束，即出现上述问题的原因是在某步消元过程之前的列选主元失败，无法选到合适的非零的适当大的主元素。解决办法，在更大的范围，即余下的n-k阶子阵中来选择主元素。在第k到n列和第k到n行的子阵中选元素绝对值最大者作主元，称为全主元消去法。,.,2. 全主元消去法,先看一个例子.求解方程组,.,全主元消去法,.,全主元消去法,.,全主元消去法,很明显，用列选主元法所得到的解要比顺序高斯消去法好，但仍然有较大的误差，下面介绍所谓的全选主元高斯消去法：,.,全主元消去法,.,全主元消去法,.,全主元消去法,.,全主元消去法,.,从此可以看出，全主元法与列主元相比，其最大的差别是需要根据列的交换情况，按照相反的顺序将解的分量的顺序重新调整！,全主元消去法,.,全主元消去法,全主元消去法的严格算法描述如下：,.,Gauss全主元消元算法,.,Gauss全主元消元算法,.,Gauss全主元消元算法,.,Gauss全主元消元算法,.,小结,比较而言,Gauss顺序消去法条件苛刻，且数值不稳定；Gauss全主元消去法工作量偏大，需要比较 个元素及行列交换工作，算法复杂。因此从算法优化的角度考虑， Gauss列主元消去法比较好。,.,思考,在上述讨论中，Gauss顺序消去法、Gauss列主元消去法、 Gauss全主元消去法的表示方法或算法描述过程中主要采用的是分量的形式。,但本章一开始，我们就漂亮地将分量形式的方程组表示为了矩阵的形式“AX=b”，而从上述过程来看，无论是消元过程，还是回代求解，都可以同样漂亮地表示为矩阵的形式。,也就是说，至少在形式上方程组的求解问题可以借助矩阵运算以及相关的理论、方法。下面就来探讨方程组求解的矩阵三角分解方法。首先，看看矩阵的三角分解！,.,2.2 矩阵的三角分解法,我们知道，对矩阵进行一次行的初等变换，就相当于用相应的初等矩阵去左乘原来的矩阵。因此我们用矩阵乘法来考察Gauss消元法。,即可得到求解线性方程组的另一种直接法：矩阵的三角分解法。,.,三角分解法解方程的原理若A = LU，则Ax = b LUx = b 其中L、U为三角形矩阵，则方程组易解定义：(1) 称 为单位下三角阵(2) 设L为单位下三角阵，U为上三角阵，若A = LU，则称A可三角（LU）分解，并称A = LU为A的三角分解（杜利特尔分解）(3) A也可以分解为一个下三角阵L和一个单位上三角阵U，此时， A的三角分解（LU）被称为克劳特分解。,.,三角分解可顺利进行的条件定理：(1) A = (aij)nn有唯一LU分解 A的顺序主子式k 0，k = 1，2，.，n 1(2) 若A = (aij)nn可逆，则存在置换阵P，使PA的n个顺序主子式全不等于0注：由Ax = b PAx = Pb知，也就是说，经适当行交换后，A总可三角分解,.,矩阵三角分解过程设A的各阶顺序主子式不为0，且A = LU，即(1) L阵主对角线（含）上边：当列标m 行标k时，lkm = 0 j = k，k+1，.，n；k = 1，2，.，n(2)U阵主对角线下边：当行标m 列标k时，ukm = 0 i = k+1，k+2，.，n；k = 1，2，.，n1,.,矩阵三角分解计算公式,设A的各阶顺序主子式不为0，且A = LU，即(1) 主对角线（含）上边：当列标m 行标k时，lkm = 0,(2) 主对角线下边：当行标m 列标k时，ukm = 0,欲求lik和ukj,.,矩阵三角分解计算公式-k=1时(1) 主对角线（含）上边：当列标m 行标k时，lkm = 0 j=k,k+1, .,n；k=1,2,.,n(2) 主对角线下边：当行标m 列标k时，ukm = 0 i=k+1,k+2,.,n；k=1,2,.,n1欲求lik和ukj设k = 1，那么，就可以计算L的第一行和U的第一列，有：a1j = l11u1j u1j = a1j j = 1，2，.，n（ U的第一行）ai1 = li1u11 li1 = ai1/u11 i = 2，3，.，n（ L的第一列）,.,矩阵三角分解计算过程-逐行（列）计算L、U一般地，逐行（列）计算L、U最后：,.,三角分解后求解两个三角形方程组下三角 Ly = b：上三角 Rx = y：紧凑格式：,.,事实上，上述过程就是高斯消去法的矩阵形式,.,第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分解算法,.,选主元的三角分解法,针对方程组Ax=b，矩阵A的三角分解能顺利进行的条件是A的各阶顺序主子式不等于零。同时，我们也已经指出，三角分解法实质上是高斯消去法的矩阵表示。因此，高斯消去法求解方程组时存在的数值稳定性问题（除数为零或很小），三角分解法同样存在。为此，有必要与高斯消去法一样进行主元的选择，以便能保证三角分解的顺利进行，使三角分解法具有好的数值稳定性。 因而就有了选主元的三角分解法，只需要在直接三角分解法中引入选主元的算法，就可以实现选主元的三角分解法。这里就不详细讨论了！,下面针对一类特殊的方程组，介绍相应的三角分解法,.,2.2.3 平方根法和改进的平方根法-对称矩阵的Cholesky（乔列斯基）分解,在科学与工程计算中，线性方程组大多数的系数矩阵具有对称正定的性质，因此基于对称正定矩阵这一特性的三角分解方法是求解对称正定方程组的一种有效方法。由于系数矩阵A的良好性质，分解过程无需选主元，具有良好的数值稳定性。 首先，我们讨论对称正定矩阵的三角分解，即Cholesky（乔列斯基）分解。,.,对称矩阵的Cholesky （乔列斯基）分解,A对称：AT=A A正定：A的各阶顺序主子式均大于零。即,.,对称矩阵的Cholesky分解,由Doolittle分解，对称正定矩阵A有唯一分解,.,对称矩阵的Cholesky分解,定理2.2.4 设A为对称正定矩阵，则存在唯一分解A=LDLT，其中L为单位下三角阵，D=diag(d1,d2,dn)且di0(i=1,n),.,对称矩阵的Cholesky分解,证明：,.,对称矩阵的Cholesky分解,所以，非奇异的上三角阵U1可化为,.,对称矩阵的Cholesky分解,.,对称矩阵的Cholesky分解,.,对称矩阵的Cholesky分解,推论：设A为对称正定矩阵，则存在唯一分解其中L为具有主对角元素为正数的下三角矩阵。LLT称为对称正定矩阵A的Cholesky分解。,.,对称矩阵的Cholesky分解,证明：,推论：设A为对称正定矩阵，则存在唯一分解其中L为具有主对角元素为正数的下三角矩阵。,.,Cholesky分解的求法,.,Cholesky分解的求法,.,Cholesky分解的求法,.,Cholesky分解法,Cholesky分解法缺点及优点 优点：可以减少存储单元。 缺点：存在开方运算，可能会出现根号下负数，同时计算量也大增。怎么改进？！,.,记得我们在由LDLT获得LLT分解时，人为地将D分为D1/2* D1/2，这是导致平方根法计算量大等问题的根源。因而，改进的想法首先就是避免这种分割。具体地：,Cholesky分解法，怎么改进？！,.,改进Cholesky分解法,基于分解A=LDLT，其中,.,改进的cholesky分解,.,改进的cholesky分解,.,.,例题,.,例题,.,例题,.,例题,.,改进平方根法,A=LDLT 分解（对应改进平方根法） ，既适合于解对称正定方程组，也适合求解A为对称，而各阶顺序主子式不为零的方程组而对A=LLT分解（对应平方根法）只适合于对称正定方程组上面，就是求解系数矩阵对称正定的方程组的平方根法和改进平方根法。下面针对另外一类特殊方程组求解的方法：追赶法！,.,2.2.4 三对角方程组求解的追赶法,.,三对角方程组求解的追赶法,.,三对角方程组求解的追赶法,.,三对角方程组求解的追赶法,.,三对角方程组求解的追赶法,其计算工作量为5n-4次乘除法（包括三角分解2n-2， 求解上三角方程n-1,回代2n-1）。工作量小.其实现的条件为qi不为零。有以下定理可得三对角方程组求解的充分性条件。,.,三对角方程组求解的追赶法,.,解三对角矩阵线性方程组的追赶法程序框图,.,从三对角矩阵带状矩阵：思考,如果方程组的系数矩阵的上半带宽为s、下半带宽为r,即,有如下三角分解的结果,.,从三对角矩阵带状矩阵：思考,有兴趣的同学请自己针对你在科学研究中碰到的实际问题，分析设计相应的三角分解算法,.,直接法时间复杂性,.,解线性代数方程组的小结,高斯消去法顺序高斯消去法列选主元高斯消去法全选主元高斯消去法三角分解法杜利特尔、克劳特三角分解法平方根法、改进的平方根法追赶法进一步需要考虑的问题：上述方法求解方程组可能存在的问题。我们从研究方程组性态的工具-向量范数、矩阵范数开始,.,2.3 向量和矩阵的范数,为了研究线性方程组近似解的误差估计和下一章要讲的迭代法的收敛性，我们需要对Rn(n维向量空间)中的向量或Rnxn中矩阵的“大小”引入一种度量：向量和矩阵的范数。,.,向量和矩阵的范数,在一维数轴上，实轴上任意一点x到原点的距离用|x|表示。而任意两点x1，x2之间距离用| x1-x2 |表示。,.,向量和矩阵的范数,而在二维平面上，平面上任意一点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不等式，证明如下,.,常见的向量范数,因而，2-范数是向量范数。,.,向量范数性质,向量范数是等价的。,.,向量范数性质,等价性质：,下面证明1）,还有其它形式的向量范数之间的等价关系，。,.,向量范数性质,等价性质1）的证明：,.,向量范数性质,等价性质证明：,.,2.3.2 矩阵范数,.,相容范数,.,相容范数,.,算子范数是矩阵范数,下面证明算子范数是矩阵范数！,.,算子范数是矩阵范数,.,算子范数是矩阵范数,.,算子范数是矩阵范数,.,常见的矩阵范数,.,常见的矩阵范数,.,常见的矩阵范数,.,常见的矩阵范数,.,例题,.,2.3.3 矩阵的谱半径和矩阵序列收敛性,.,谱半径和矩阵序列的收敛性,.,例题,.,例题,.,.,2.4 病态方程组与矩阵的条件数,.,2.4.1 病态方程组与扰动方程组的误差分析,.,病态方程组与扰动方程组的误差分析,.,病态方程组与扰动方程组的误差分析,.,病态方程组与扰动方程组的误差分析,.,病态方程组与扰动方程组的误差分析,.,病态方程组,扰动方程 由于计算机字长限制，在解AX=b时，舍入误差是不可避免的。因此我们只能得出方程的近似解 。 是方程组（A+A）x=b+b (1) 这就是所谓的扰动方程。,.,病态方程组,称（A+A）x=b+b为方程Ax=b的扰动方程。其中A， b为由舍入误差所产生的扰动矩阵和扰动向量。 当A， b有微小扰动，解得扰动方程（1）的解与Ax=b的解x的相对误差不大时，称为良态方程，否则为病态方程。 下面运用前面介绍的向量范数、矩阵范数来研究方程组解得误差传播规律。,.,扰动方程组的误差界-右端项扰动b,.,扰动方程组的误差界-系数矩阵扰动A,.,扰动方程组的误差界-系数矩阵扰动A 、右端项扰动b,.,总结将上述三种情况，我们就能得到在系数矩阵有扰动A 和/或右端项有扰动b时，方程组的解的误差限！在上述三种情况中都出现了一个决定误差界大小的量，即|A|A-1|。也就是这个数决定了方程组得性态，我们将它称之为条件数。,.,2.4.2 矩阵的条件数,.,矩阵的条件数的性质,.,扰动方程组的误差界-定理,.,用条件数来刻画方程组的性态,.,例题，计算矩阵的条件数,这里，我们用了条件数的性质来计算对称矩阵的条件数，也可以用条件数的定义来计算。即先求A的逆矩阵，再求A的条件数！,.,例：Hilbert 阵,注：现在用Matlab数学软件可以很方便求矩阵的条件数!,对于严重的病态线性方程组，即使原始数据A和b都没有误差，但如果在求解过程中有舍入误差，所得到的解也会有很大的相对误差。,.,什么样的线性方程组可能是病态的？,注：一般判断矩阵是否病态，并不计算A1，而由经验得出。 行列式很大或很小（如某些行、列近似相关）； 元素间相差大数量级，且无规则； 主元消去过程中出现小主元； 特征值相差大数量级。,.,缓和甚至解决线性方程组病态的手段,采用高精度的算法（不是总有效）；通过行或列的平衡来降低矩阵的条件数（平衡法）； 残差校正法；通过矩阵的预条件处理来降低矩阵的条件数。,.,残差校正法步骤：,求解Axb，得到x的近似解x1；校正x1，令r1bAx1，解A x1 r1；x2x1 x1 ；令r2bAx2，解A x2 r2；x3x2 x2 ；。令rkbAxk，解A xk rk；Xk1xk xk ；。,.,预条件技术,80年代提出，纯代数观点M 近似 A, 求解(左预条件)M y = c 易解 相当于化学反应中寻找高效、廉价的催化剂，能极大地提高迭代方法的速度。,.,预条件技术(续),代数预条件技术ILU、SPAI、SOR、多项式几何预条件技术某方向压缩粗化、网格均匀化、区域规则化近似分析预条件技术变系数常数化、Green函数稀疏近似、强化椭圆型物理预条件量纲平衡、物理参数逼近、状态方程近似、某些物理项简化、算子分裂高级预条件MG、DDM、快速变换（FFT）、基底变换法,.,解线性代数方程组的迭代解法,.,In Scientific ComputingLarge Linear SystemsAx=bas sub-problems/as intermediate steps,Conjugate Gradient methodfor symmetric systems,.,173,直接法和迭代法的比较,直接法: 经过有限次运算后可求得方程组精确解的方法(不计舍入误差!),迭代法：从解的某个近似值出发，通过构造一个无穷列去逼近精确解的方法。（一般有限步内得不到精确解）,直接法比较适用于中小型方程组。对高阶方程组，既使系数矩阵是稀疏的，但在运算中很难保持稀疏性，因而有存储量大，程序复杂等不足。 迭代法则能保持矩阵的稀疏性，具有计算简单，编制程序容易的优点，并在许多情况下收敛较快。故能有效地解一些高阶方程组。,.,解线性代数方程组的迭代法,引进迭代法的目的：解大型方程组；存储的考虑；解病态方程 。,.,2.5 解线性代数方程组的迭代法,迭代发的一般理论,Jacobi迭代法、G-S迭代法,松弛迭代法,迭代法的收敛条件,迭代法收敛的其它判定方法,解线性方程组的迭代法,解线性代数方程组的极小化方法,.,2.5 解线性代数方程组的迭代法,迭代法的一般理论,Jacobi迭代法、G-S迭代法,松弛迭代法,解线性方程组的迭代法,2.5.1 一般迭代法,.,引例,解,迭代法的基本思想,.,或写为 ，,其中,迭代法的基本思想,.,迭代法的基本思想,.,迭代法的基本思想,.,我们任取初始值，例如 ，将这些值代入（2）式右边（若（2）式为等式即求得方程组的解，但一般不满足），得到新的值,得到向量序列和一般的计算公式（迭代公式）,迭代法的基本思想,.,其中k表示迭代次数（k=0,1,2, ）.,迭代到第10次有,迭代法的基本思想,.,即建立一种从已知近似解到新的近似解的计算公式。,如何建立这样的迭代计算公式呢？,迭代法的基本思想,迭代法求方程组的解的基本思想即是：构造一串收敛到解的序列，,.,思路,基本迭代法,.,基本迭代法,.,将线性方程组 Ax=b 化为 x=Mx+g,再由此构造向量序列x (k): x(k+1)=Mx (k)+g,若x (k)收敛至某个向量x *,则可得向量x *就是所求方程组 Ax=b 的准确解.,迭代法的基本思想,我们已经回答了 “如何建立这样的迭代计算公式呢？”的问题。方法是：,下面就来构造我们这一章中将要介绍的第一种迭代法,.,Jacobi迭代法,.,Jacobi迭代法,.,Jacobi迭代法,.,Jacobi迭代法,.,求解一般线性代数方程组的Jacobi迭代法的计算公式,Jacobi迭代法,.,举 例,按照jacobi法的一般计算公式写出本例的jacobi迭代式,.,Jacobi迭代法的特点,Jacobi迭代法的特点,.,Jacobi迭代法,.,Gauss-Seidel迭代法,.,Gauss-Seidel迭代法,.,Gauss-Seidel迭代法,.,即：,Gauss-Seidel迭代法,.,用高斯-赛德尔迭代法解方程组,例4.1,Gauss-Seidel迭代法举例,.,如此继续下去,Gauss-Seidel迭代法的特点,.,简单，每迭代一次只需进行一次矩阵和向量乘法， 不改变系数矩阵的稀疏性， 仅需要一组存储单元 。 Gauss-Seidel迭代法一般比Jacobi迭代法快，但也有比 Jacobi迭代法慢的，甚至有Jacobi迭代法收敛而Gauss- Seidel迭代法不收敛的。,Gauss-Seidel迭代法特点,.,举 例,请按照下述Gauss-Seidel法的一般计算公式写出该例的Gauss-Seidel迭代式,.,上面我们由介绍了一般迭代法的构造方法，又由一般迭代方法具体引出了Jacobi迭代法，进一步考虑迭代过程中新信息的利用，得到了Gauss-seidel迭代法。当然，这两种特殊的迭代法同样可以直接由系数矩阵的和、差分解导出。下面我们在分析这两者方法存在误差的基础上，提出一种新的方法，即所谓松弛迭代法。首先来看松弛迭代法的基本思想。,松 弛 法,.,松弛法的基本思想,我们从为Gauss-Seidel 迭代法加速的角度来考虑提出松弛法。,Gauss-Seidel 迭代公式得到，于是有,这样我们就可以得到这两次近似解的差别x,即,.,松弛法的基本思想,具体地,.,可以把 看作Gauss-Seidel 迭代的修正项，即第k次,近似解 以此项修正后得到新的近似解。,我们有理由相信，新的近似解，应有更好的精度，误差更小。当然，考虑到x也是利用两次近似估计得来的，并不是真正的第k次近似解的误差，所以在实际修正时，更好的做法是对修正项x加权！具体地有：,松弛法的基本思想,具体地，我们,.,得到新的近似解，具体公式为,松弛法是将 乘上一个参数因子 作为修正项而,松弛法的基本思想,所以，松弛法就可以进一步表示为：,.,称为松弛因子,当 时称为低松弛；,时称为超松弛法。,松弛法的基本思想,上面我们运用误差估计校正法，得到了松弛法的分量形式，这已经可以方便地让我们编程计算了，但分量形式不利于分析推理，下面将松弛法表为矩阵形式：,.,选取分裂矩阵M为带参数下三角阵,其中 为可选择的松弛因子.,松弛法的矩阵表示,.,松弛法的矩阵表示,通过简单的推导，得到,.,松弛法的矩阵表示,对照前面给出的分量形式！,.,SOR方法的计算公式,SOR 方 法,对照前面给出的分量形式！,.,例3,取 用松弛法解方程组,举 例,解：由,.,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.59984,1.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 1.61.2 1.4 1.6,计算结果,松弛法的程序,.,验证： 松弛因子与迭代次数的关系,例4,举 例,.,举 例,.,举 例,.,举 例,请按照下述松弛法的一般计算公式写出该例的松弛迭代式,.,2.5.2 迭代法的收敛条件,1、向量序列收敛与矩阵序列收敛,2、迭代法收敛概念,3、迭代法收敛定理,4、举例,迭代法的收敛条件,.,不一定!,定义4.1,关注的问题及几个相关概念,.,关注的问题及几个相关概念,.,证明：,上述定理表明，向量序列和矩阵序列的收敛可归结为对应分量或对应元素序列的收敛。,矩阵有类似结果,关注的问题及几个相关概念,.,定义4.2,迭代矩阵,.,定义4.3,矩阵的谱半径,回顾我们上一章中学习过的矩阵的谱半径的定义,谱半径是我们讨论迭代法收敛性的重要工具！它具有以下性质：,.,矩阵的谱半径,.,零矩阵的充要条件,有了谱半径，就可以给出前面定义的向量序列、矩阵序列收敛的条件了！当然也就可以讨论迭代法的收敛性问题了。,.,定理4.1,证明：（必要性）,零矩阵的充要条件,.,（充分性）,零矩阵的充要条件,.,零矩阵的充要条件,有了上述向量序列、矩阵序列收敛的条件，也就等价地有了-一般迭代法收敛条件了。我们有迭代法的如下基本定理：,.,定理4.2,证明：,迭代法的收敛定理,.,零矩阵的充要条件,有了一般迭代法收敛的充要条件，从理论上解决了迭代法收敛性的判定问题。但收敛的速度、误差传播的规律或者误差估计的问题仍然没有解决。,还记得我们有 的结果吗？如果我们能将迭代矩阵的条件加强为 ，我们应该有希望得到我们想要的误差估计：,.,和,误 差 估 计,.,误 差 估 计,.,有,由此可得,误 差 估 计,.,所以,误 差 估 计,由于,.,取范数,误 差 估 计,证毕,.,误差估计式,给出了停止迭代的一个判别准则，即用相邻两次迭代结果的差的范数来判别，它越小，迭代序列与精确解越接近。,由误差估计式,误 差 估 计,.,例：方程组,解：用雅可比迭代法,误 差 估 计 举 例,.,误 差 估 计 举 例,.,迭代法收敛性的判定,2、若A为严格对角占优阵或不可约弱对角占优阵，则Jacobi迭代法和Gauss-Seidel迭代法均收敛。,小 结,.,三种迭代法收敛性判定的其它定理,前面，我们已经给出了一个一般迭代法收敛判定的充要条件和一个迭代矩阵的某种范数小于1的特定迭代法收敛的充分条件。下面，我们针对一些更特别的方程组的迭代法来给出几个收敛判定定理。,.,若n阶方阵 满足,定义1,定义2,如果矩阵A不能通过行的互换和相应的列互换成,为形式 ，其中 为方阵，则称A为,不可约。,且至少有一个i值,使得上式中不等号严格成立，则称A,为弱对角占优阵。若对所有i,上式不等号均严格成立 ，,则称A为严格对角占优阵。,定 义,例,.,定理,定 理,设有线性方程组 ，下列结论成立,1. 若A为严格对角占优阵或不可约弱对角占优阵，则,Jacobi 迭代法和Gass-Seidel迭代法均收敛。,2.若A为严格对角占优阵， 则松弛法收敛。,3.若A为对称正定阵，则松弛法收敛的充要条件为,上两例中：,A为严格对角占优阵， Jacobi 迭代法和Gass-Seidel迭代均收敛。B为非严格对角占优阵， 故松弛法收敛。,.,例5,解：因A为对称且各阶主子式皆大于零，故A为对称正定,均收敛。A不是弱对角占优阵，故不能用条件1判断。,阵。由判别条件3，Gauss-Seidel迭代法与松弛法,Jacobi迭代法迭代矩阵,举 例,.,其特征方程,举 例,.,Jacobi与Gauss-Seidel迭代的迭代矩阵分别为,改变方程组中方程的次序，即将系数矩阵作行,交换，会改变迭代法的收敛性,例,举 例,.,若交换方程的次序，得 Ax=b 的同解方程组,为严格对角占优矩阵，因而对方程组 用,Jacobi 与Gass-Seidel 迭代求解均收敛。,举 例,.,对称超松弛迭代法（SSOR法）,.,对称超松弛迭代法（SSOR法）,.,SSOR迭代法的矩阵形式：,注：（1）关于 的收敛条件和准则与SOR方法相同；（2）收敛快慢对 的选取不敏感。,.,块超松弛迭代法（BSOR法）,.,案 例 1 (热传导问题),设有一维热传导方程的初边值问题,试用数值方法求出t=0.2时刻金属杆的温度分布.,应用实例,.,解:(1)对空间进行离散.,(2)对微分算子进行离散.,t0.020.010,0 0.1 0.2 . 0.9 1 x,.,采用无条件稳定的Crank-Nicholson格式,则有,.,或,加上边界条件后有,.,加上边界条件后有,其矩阵形式为 ATn+1=dn,.,案 例 2,数值求解正方形域上的Poisson方程边值问题,.,解:(1)剖分求解域.,(2)对微分算子进行离散.,0 1 2 . N N+1 X,YN+1N：210,.,在每个点(xi,yj)上的有限差分方程为,在边界上,又称为五点差分格式,.,对非边界点进行编号: 顺序为-从下往上,从左往右,相应的解向量和右端向量分别为,.,.,.,差分方程组的矩阵形式为 Au=f,如果把每一条线上的网点看作一个组,如,.,其中,可用块迭代法求解Au=f,.,.,2.5.3解线性代数方程组的极小化方法,一、与线性方程组等价的变分问题,三、共轭斜量法(共轭梯度法),四、预条件共轭斜量法,二、最速下降法,.,共轭梯度法最早是由计算数学家Hestenes 和几何学家Stiefel 在20世纪50年代初为求解线性方程组而各自独立提出的。他们合作的文章被公认为共轭梯度法的奠基之作，详细讨论了求解线性方程组的共轭梯度法的性质以及它和其他方法的关系。在A为对称正定阵时，上述线性方程组等价于最优化问题。,解线性代数方程组的极小化方法,.,由此，Hestenes和Stiefel的共轭梯度方法也可视为求二次函数极小值的共轭梯度法。1964年Fletcher和Reeves 将此方法推广到非线性优化，得到求一般函数极小值的共轭梯度法。,解线性代数方程组的极小化方法,.,提到最优化问题，首先介绍何为最速下降法。考虑线性方程组Ax=b的求解问题，其中A是给定的n阶对称正定矩阵，b是给定的n维向量。为此我们定义二次泛函,解线性代数方程组的极小化方法,.,定理：设A对称正定，求方程组Ax=b的解，等价于求二次泛函 的极小点由定理，求解线性方程组的问题就转化为求二次泛函的极小点的问题，也可表为与线性方程组等价的变分问题：,解线性代数方程组的极小化方法,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,一、与线性方程组等价的变分问题,.,二、最速下降法,.,二、最速下降法,.,二、最速下降法,.,二、最速下降法,.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,公式化简,.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,公式化简,三、共轭斜量法 (CG) (共轭梯度法),.,三、共轭斜量法 (CG) (共轭梯度法),.,.,四、预条件共轭斜量法(PCG),.,预条件共轭斜量法,实际计算，可通过变换，转化成用原方程组的量来计算。,.,.,预优矩阵的选取,下面介绍几种选取预优矩阵的方案:,.,几种选取预优矩阵的方案,.,几种选取预优矩阵的方案,.,几种选取预优矩阵的方案,.,三种预条件共轭斜量法在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.未给定预优矩阵的共轭斜量法 r=cholinc(sa,0) x=pcg(a,b,tol,kmax,r,r,x0),.,CG(Conjugate Gradient)（对称正定）MRES(Minimal RESidual) (对称不正定）GMRES(Generilized Minimal RESidual） （不对称不正定）BiCG(不对称不正定）BiCGSTAB (不对称不正定）,几种常用的Krylov-subspace方法,.,背景,大规模线性方程组的求解很多科学工程计算问题都转化为求解方程组Ax=b.如偏微分方程组的差分格式，有限元方法离散得到刚度矩阵.大规模矩阵特征值和特征向量的计算工程计算领域十分常见。如量子物理中的Kohn-Sham方程求解化为哈密顿矩阵某些关键特征值对的计算.,Krylov-subspace方法,.,投影方法,线性方程组的投影方法方程组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),K,L,X(0),X(1),.,解方程组的投影法的矩阵表示设nm阶矩阵V=v(1), v(2), v(m)与W=w(1), w(2), w(m)的列分别构成K与L的一组基 .记z=x(1)-x(0),z=Vy,有 当非奇异时有从而得到迭代公式,.,投影方法的最优性. (误差投影)设A为对称正定矩阵, x(0)为初始近似解,且K=L,则x(1)为采用投影方法得到的新近似解的充要条件是 其中, 且x为问题的精确解.,.,. (残量投影)设A为任意方阵,x(0)为初始近似解,且L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是 其中,.,矩阵特征值的投影方法 对于特征值问题Ax=x,其中A是nn的矩阵，斜交投影法是在m维右子空间K中寻找xi和复数i满足Axi- ixi L,其中L为m维左子空间.当L=K时，称此投影方法为正交投影法.,.,Kyrlov子空间,定义为m维Krylov子空间为v的次数，即使得q(A)v的非零首一多项式的最低次数,.,Krylov子空间的标准正交基,Arnoldi方法基于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的不变子空间.,.,.,定理 如果记以v(1),v(2),v(m)为列构成的矩阵为Vm,由hij定义的(m+1)m阶上Hessenberg矩阵为 删除最后一行得到的矩阵为Hm,则在Arnoldi算法中，可能有较大舍入误差，改写,.,.,Krylov子空间方法解线性方程组,误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法,.,.,不难看出，当采用上述FOM算法时，需要存储所有的vi，(i=1,2,m)，当m增大时，存储量以O(mn)量级增大.而FOM计算量是O(m2 n).可见其代价十分高昂.因此我们考虑重启的FOM算法,.,.,Krylov子空间方法解线性方程组,误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法 2) 非对称矩阵的IOM方法和DIOM方法,.,所谓不完全正交化方法(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子空间方法解线性方程组,误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法 2 ) 非对称矩阵的IOM方法和DIOM方法）对称矩阵Lanczos算法,.,A是非对称阵时,H是Hessenberg阵，则当A是对称阵时，H也为对称阵，故应为三对角阵，记为对它采用修正Arnoldi正交化过程，就得到Lanczos方法,.,.,Krylov子空间方法解线性方程组,误差投影型方法取L=K时的正交投影法）非对称矩阵的FOM方法 2 ) 非对称矩阵的IOM方法和DIOM方法）对称矩阵Lanczos算法 4 ）正定对称阵的CG法,.,CG法是解正定对称线性方程组最有效的方法之一 该方法充分利用了矩阵A的稀疏性，每次迭代的主要计算量是向量乘法而实际上，CG方法也是一种Krylov子空间方法 后面将给出一个数值例子,.,.,Krylov子空间方法解线性方程组,残量投影型方法取L=AK时的斜交投影法）GMERS方法,.,GMERS方法把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法,.,回顾,. (残量投影)设A为任意方阵,x(0)为初始近似解,且L=AK,则x(1)为采用投影方法得到的新近似解的充要条件是 其中,.,GMERS方法在Arnoldi正交化过程之后把方程组的求解化为一个极小问题的求解，回顾之前介绍投影类方法不去直接求x(1),转而去解此极小问题，这是GMERS方法的独到之处.再由之前的定理,.,.,.,Krylov子空间方法解线性方程组,残量投影型方法取L=AK时的斜交投影法）GMERS方法 2 ) 重启型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子空间方法解矩阵特征值问题,Arnoldi方法求解非对称矩阵特征值Lanczos方法求解对称矩阵特征值,.,Arnoldi方法 再次回顾前面的定理而 所以有以下结论：,.,1）若 ,则 是A的不变子空间,从而Hm的所有特征值是A的特征值子集 2）若 ,设Hm的特征值对是 ,即 ,由上述定理可得,.,以上讨论说明,可以用Hm的特征值 (又称Ritz值)来近似A的特征值,用向量 (又称Ritz向量)来近似相应的特征向量.事实上, Hm为A在Kyrlov子空间上的正交投影. 由于Hm是上Hessenberg阵，它的特征值求解难度远小于原问题,例如可以采取分解法来求解.,.,.,Arnoldi算法构造标准正交基 1需要存储所有的基向量,当m很大时,存储量大 2理论上为了保证收敛速度,m越大越好 矛盾！,.,Lanczos 方法 A是对称阵时,Hm是三对角阵，仍然采用Hm的特征值来近似A的特征值，相应的Ritz向量来近似A的特征向量.问题转化为三对角阵特征值的求解问题，而它的求解，复杂度大大减小，有很多成熟的办法，如分而治之法，QR法等，可以在并行机上达到O(logN)的量级.,.,.,对称Lanczos算法,在没有舍入误差的情形下,得到的是标准正交基相互正交的,并且至多n步必然终止.但是在出现舍入误差的情况下，计算得到的将失去正交性.因此,长期以来被人们认为是不稳定的.直到1971年,Paige通过仔细的舍入误差分析,发现失去正交性恰与近似特征值的精度提高有关.之后,经过大量的理论分析和数值实验,人们充分认识到Lanczos方法是求解大型稀疏对称矩阵特征值问题的最有效方法之一. 目前，Lanczos方法是求大型稀疏对称矩阵部分极端特征值的最有效方法,
展开阅读全文