Chapt线性代数方程组的数值解法实用.pptx

4.1 4.1 概概 述述 求解线性代数方程组的数值计算方法很多，大致可分为两大类：消去（元）法和迭代法。消去法是直接从方程组的系数矩阵入手，经过有限步运算求出方程组的精确解（假如没有舍入误差的话）。迭代法则是将求方程组的问题化为构造一组递推计算结构，从一组近似解出发，用这组递推结构逐步算出精度更高的近似解。上述这两类算法各有其优点和缺点，消去法的计算量小，但程序复杂；迭代法计算量大，精度不高，但程序结构简单。本章将主要介绍求解线性代数方程组的简单高斯消去法和三角分解法，迭代法将在以后相关章节中介绍。第1页/共45页4.2 4.2 高斯消去法高斯消去法 4.2.1高斯消去法的基本步骤我们以三元线性代数方程组为例，叙述简单高斯消去法（以下简称为消去法）的基本步骤，这各方法是理解其它方法的基础，消去法分为消元和回代两个过程。例4-1：消去过程实际上是对增广矩阵作行初等变换。上例可表示为 第2页/共45页4.2 4.2 高斯消去法高斯消去法 这是上三角方程组，它极容易求解：由第三个方程得 ，代入第二个方程得 ，再代入第一方程 。此一过程称为回代过程。从上述消元过程可以看出，三元议程组要经过两次消元（因为不用消）才能把增广矩阵化为拟上三角矩阵。对于一般的n元线性代数方程组，经经过(n-1)次消元才能把相应的增广矩阵变换为拟上三角矩阵。n元线性代数方程组的增广矩阵的一般形式如(4.3)式。第3页/共45页4.2 4.2 高斯消去法高斯消去法 以下是消元步骤（共分n-1部）：第一步(k=1)：从第一列中消去第2n行中x1的系数，即消去a11下方元素；第二步(k=2)：从第一列中消去第3n行中x2的系数，即消去a22下方元素；(4.3)第4页/共45页4.2 4.2 高斯消去法高斯消去法第j步(k=j)：从第j列中消去第j+1n行中xj的系数；即消去ajj下方元素;第n-1步(k=n-1)：从第n-1列中消去第n行中xn-1的系数；即消去an-1n-1下方元素因此，高斯消元步骤的计算通式可写为：(4.4)第5页/共45页4.2 4.2 高斯消去法高斯消去法经过n-1步消元后，增广矩阵(4.3)式变为：(4.5)第6页/共45页4.2 4.2 高斯消去法高斯消去法高斯回代过程的通用递推计算公式可写为：(4.6)第7页/共45页4.2 4.2 高斯消去法高斯消去法4.2.2高斯消去法的计算量分析我们从(4-4)式可以看出，高斯消去法消去过程共分n-1步，第k步变换n-k行：对这些行先乘系数，再从n+1-k元素减去第k行相对应列的元素；因此共需乘、除次数为回代过程求xn需1次除法，求xn-1需1次乘法、1次除法，求x1需n-1次乘法、1次除法，因此共需乘除次数两过程共需要乘除次数为 。当n=20时，N=3060.第8页/共45页4.2 4.2 高斯消去法高斯消去法4.2.2 高斯消去法的程序框图与通用程序 STARTINPUT N,A(),BCALL eliminationCALL backOUTPUT X(N)END第9页/共45页4.2 4.2 高斯消去法高斯消去法10 OPEN(8,&file=EliminationMetod.dat,&status=new)20 DIMENSION a(n,n+1),x(n),b(n)30C 输入n,a(),b40 READ(8,(1x,i2,)n50 DO i=1,n60 READ(8,(1x,e12.6)b(I)70 END DO80 DO i=1,n90 DO j=1,n100 READ(8,(1x,e12.6)a(i,j)110 END DO 120 END DO 120C 消元过程130 CALL elimination(a,b,n)140C 回代过程150 CALL back(a,x,n)160C 输出结果x170 DO i=1,n180 WRITE(9,(1x,e12.6)x(i)190 END DO 200 END210220230240250第10页/共45页4.2 4.2 高斯消去法高斯消去法SUBROUTINE eliminationK=1,n-1i=k+1,nR=a(i,k)/a(k,k)j=k+1,n+1a(i,j)=a(i,j)-R*a(k,j)end jend jend iend kEND SUBROUTINE elimination第11页/共45页4.2 4.2 高斯消去法高斯消去法10 SUBROUTINE&elimination(a,b,n)20 DIMENSION a(n,n+1),b(n)30C 生成增广矩阵40 DO i=1,n50 a(i,n+1)=b(i)60 END DO70C 消去过程开始80 DO k=1,n-190 DO i=k+1,n100 R=a(i,k)/a(k,k)110 DO j=k+1,n+1120 a(i,j)=a(i,j)-R*a(j,k)130 END DO140 END DO150 END DO160 END SUBROUTINE elimination第12页/共45页4.2 4.2 高斯消去法高斯消去法SUBROUTINE backX(n)=a(n,n+1)/a(n,n)k=n-1,1,-1s=0j=k+1,ns=s+a(k,j)*x(j)end jX(k)=(a(k,n+1)-s)/a(k,k)end kEND SUBROUTINE back第13页/共45页4.2 4.2 高斯消去法高斯消去法10 SUBROUTINE back(a,x,n)20 DIMENSION a(n,n+1),x(n)30 x(n)=a(n,n+1)/a(n,n)40 DO k=n-1,1,-160 s=0.0 60 DO j=k+1,n70 s=s+a(k,j)*x(j)80 END DO90 X(k)=(a(k,n+1)-s)/a(k,k)100 END DO110 END SUBROUTINE back第14页/共45页4.2 4.2 高斯消去法高斯消去法4.2.3 选主元法高斯消去法消去过程中，第k步求n-k个系数aik/akk用到的除数akk,称为主元。如果akk为接近于零或零，则消元递推计算通式(4-4)将出现被零除的情况，从而使计算机溢出或无法进行下去。例4-2准确到九位小数的解是x1=0.250001875，x2=0.499998749，若在4位计算机上按高斯消去法求解，则第15页/共45页4.2 4.2 高斯消去法高斯消去法 回代得解 ，显然产生严重失真。造成这种结果的原因，就是小主元10-5的出现。用它作除数会产生大乘数，数乘大乘数变大数，大数吃小数产生舍入误差,如()，从而使 也有误差。因此，为了保证获得正确的结果，在消元过程中可在第k步的第k列的元素akk,ak+1,k,，ank中选主元，即在其中找出绝对值最大的apk，然后交换第k行和第p行，继续进行消去过程。这种消去法为列主元消去法。也可在第k步的第k行的元素akk,ak,k+1,，akn中选主元，即在其中找出绝对值最大的akq，交换第k行和第q行;或在kn行、kn列选绝对值最大的元素 apq，交换k、p行和k、p列，然后继续进行消去过程。这两种消去法为行主元消去法和全主元消去法。本节讨论应用广泛且易于在计算机上实现的列主元素消去法。列主元素消去法的计算步骤可归纳如下：选主元素；把主元素所在的行与第k行互换；消元；回代 第16页/共45页4.2 4.2 高斯消去法高斯消去法其中、两步与简单高斯消去法中的消元与回代过程完全一样，也就是说，列主元消去法是在简单消去法的基础上引入选列主远过程，这个过程并不影响获得正确结果。列主元高斯消去法的程序框图与通用程序 为了简单和说明问题，将上面的和两步单独用一个子程序。因此，列主元高斯消去法的程序框图只需在简单高斯消去法的流程图中增加一个选列主元子程序gaussp，具体如下：第17页/共45页4.2 4.2 高斯消去法高斯消去法SUBROUTINE gausspk=1,n-1 a(p,j)=ti=k+1,nc=abs(a(i,k);p=iend it=a(k,j);a(k,j)=a(p,j)end jEND SUBROUTINE gausspc=abs(a(k,k);p=kIf abs(a(i,K)c?YNj=1,n+1 end k第18页/共45页4.2 4.2 高斯消去法高斯消去法10 SUBROUTINE gaussp(a,b,n)20 DIMENSION a(n,n+1),b(n)30C 生成增广矩阵40 DO i=1,n50 a(i,n+1)=b(i)60 END DO70C 选主元开始80 DO k=1,n-190 c=abs(a(k,k);p=k100 DO i=k+1,n110 IF(abs(a(i,k)c)THEN120 c=a(i,k);p=i130 END IF140C p,k行交换开始150 DO j=1,n+1160 t=a(k,j);a(k,j)=a(p,j)170 a(p,j)=t180 END DO 190 END DO200 END SUBROUTINE gaussp第19页/共45页4.2 4.2 高斯消去法高斯消去法4.2.4 高斯-约当法前面所讲的是如何将一个线性方程组的系数矩阵改变为一个上三角矩阵，然后回代求解。本节所讲的是如何将系数矩阵改变为一个对角矩阵，即，除对角线元素外，其他元素均为零 第20页/共45页4.2 4.2 高斯消去法高斯消去法列主元高斯-约旦消元过程的递推计算通式结构为：高斯-约旦消去法消元过程的计算量比高斯消去法有所增加（约为1.3倍），但是总的计算量增加并不大，另一方面，这一方法的程序结构简单，所以在实际中得到了广泛的应用。高斯-约旦消去法的程序框图与通用程序自行设计。改进后(4.8)(4.7)第21页/共45页4.3 4.3 三角分解法三角分解法 线性方程组的另一直接解法是三角形分解法，即方程组（4-2）的系数矩阵A分解成两个形式简单的三角形矩阵L和U的乘积：A=LU。从而求解Ax=b的问题转化为三角形方程组 Lx=y Uy=b 其中(4.9)(4.10)第22页/共45页4.3 4.3 三角分解法三角分解法比较等式左边和右边乘积矩阵LU的第r行主角元右边（含主角元）的对应元素，得再比较等式两边第r列主角元以下（不含主角元）的对应元素，得当r=1时，有(4.11)(4.12)(4.13)(4.14)第23页/共45页4.3 4.3 三角分解法三角分解法假定已求出U的第1至第n-1行和 L的第1至第n-1列，由（4-11）式计算出U的第r行元素和由（4-12）式计算出L的第r列元素（4-13）式至（4-16）式所表示的矩阵分解称Doolittle分解。类似地，若U为单位上三角矩阵，而L为下三角矩阵，则有(4.15)(4.16)(4.17)第24页/共45页4.3 4.3 三角分解法三角分解法规定 。称上述分解为Crout分解。实现了系数矩阵A的Doolittle或Crout分解后，方程组Ax=b可以分解成为Ly=b和Ux=b这两个三角矩阵来进行求解。设U为上三角矩阵，而L为单位下三角矩阵，并且 ，根据公式（4-4）和（4-6）得与Doolittle分解相对应的方程组解为(4.18)(4.19)第25页/共45页4.3 4.3 三角分解法三角分解法类似地，与Crout分解相应的求解公式(4.20)(4.21)第26页/共45页4.3 4.3 三角分解法三角分解法无论Doolittle分解还是Crout分解，其运算量均小于高斯消去法。为保证它们能顺利、稳定进行，也可选列主元。用Crout分解法求解方程组Ax=b的步骤如下：用（4-17）式求lir;确定最大列主元；用（4-18）式求uir,进行分解；回代。通用公式为(4.22)(4.23)第27页/共45页4.3 4.3 追赶法追赶法求解三对角占优方程组若作Crout分解，则有第28页/共45页4.3 4.3 追赶法追赶法比较Crout解，易知回代得 。(4.24)(4.25)第29页/共45页4.3 4.3 追赶法追赶法按照以上公式求解Ax=b的方法就叫做追赶法，其中ui、li 和yi 称追，回代称赶，乘除法运算量仅为5n-4，远比一般方程组的高斯消去法或三角分解法的运算量少。不用选主元，就可保证顺利、稳定进行计算。第30页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响 本节研究线性代数方程组解的误差。衡量数的误差用绝对值。方程组的解是向量，衡量其误差，自然也会想到绝对值的概念。4.4.1 向量与矩阵范数 众所周知，数x的绝对值x是x的函数:x=(x),具有以下三性质：（1）（2）对任意实数（3）推广到向量 ，具有如下类似性质的函数：。该函数成为向量 x的范数或模或长度，它也应该具有以下性质：（1）非负性（2）齐次性：对任意实数（3）三角不等式：第31页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响几何上三角不等式表示：在以向量x、y及其和x-y构成的三角形中，一边之长不超过另两边边长之和，即从而表示以x、y和x-y为边的三角形中，一边之长不小于另两边边长之差。常用向量范数有3种，即：分别称为1范数、2范数和无穷大范数。不难验证它们具有性质（1）和（2）。对于性第32页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响性质（3），这里已2范数为例，给出证明，以供参考。最后不等式称Cauchy不等式。注意对任意实数 第33页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响由 二次式的判别式 ，便可推出Cauchy不等式，因而2范数的性质（3）得证。对于n阶方阵，具有如下4种性质的函数 称为矩阵范数。（1）（2）对任意实数（3）（4）容易验证矩阵函数具有上述4种性质，因而，它是一种范数，称为矩阵的F范数。常用矩阵范数是如下三种范数：第34页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响其中 表示矩阵 的谱半径，即 特征值绝对值的最大值。三种范数分别称为1范数或列范数，无穷大范数或行范数，2范数或谱范数。可以证明，矩阵的这三种范数还具有以下3种性质：（满足性质（1）的矩阵范数成为相容范数）（1）（2）单位矩阵I的范数（3）时 可逆且 (4.26)第35页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响例4-3 设则第36页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响它的三个根为4.4.1 误差分析 设线性方程组Ax=bAx=b中，A A为非奇异矩阵，x x为方程组的精确解。假定b b有误差b b,则解为x+x+x x,即第37页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响 即两边取范数又因由（4-30）和（4-31）式，得另外，若矩阵A A有误差A A，则解为x+x+x x，即(4.27)(4.28)(4.29)(4.30)(4.31)(4.32)(4.33)第38页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响即有虽然A A为非奇异，但不能保证A+A+A A也为非奇异，为保证A+A+A A为非奇异，需对 A A作些限制。因为由（4-26）式，若 ，则 存在，因此（4.35)式变为由（4-34）和（4-35）式，得两边取范数(4.34)(4.35)(4.36)(4.37)(4.38)第39页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响设 ，则得从（4-32）和（4-39）式可以看出，如果线性方程组Ax=bAx=b中，A或的微小变化，都将引起方程解x为 倍，若这个倍数巨大，则方程组定义为“病态方程组”，矩阵A A为“病态矩阵”，否则定义为“病态方程组”或矩阵A A为“良态矩阵”。此时，我们称 为条件数，即常用的条件数有(4.39)(4.40)(4.41)第40页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响例4-4 试计算 ，并取6位有效数字。解（1）直接计算若使用Gauss列主元消去法，经消元则可得第41页/共45页4.4 4.4 舍入误差对解的影响舍入误差对解的影响（2）间接计算，在上线性方程组两边同乘以矩阵则有故用Gauss主元消去法得第42页/共45页习习 题题4-1 用Gauss消去法求解下列方程组Ax=b4-2 求cond(A)2和cond(A).第43页/共45页上机实习题4-3 用选列主元的高斯消去法求解下列方程组：第44页/共45页感谢您的欣赏！第45页/共45页