[理学]数值分析第六章.ppt

数值计算方法,任课教师: 徐昱 工程学院海洋工程系 1201,第五章,解线性方程组的数值解法,引言,在自然科学和工程技术中,很多问题归结为解线性方程组.例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵(例如，阶数不超过150). 另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多).,关于线性方程组的解法一般有两大类：,但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解. 本章将阐述这类算法中最基本的和具有代表性的算法就是高斯消元法,以及它的某些变形和应用.这类方法是解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组(例如，大型带状方程组)的有效方法.,1. 直接法,经过有限次的算术运算,可以求得方程组的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法对高阶方程组计算量太大,不是一种实用的算法.,就是用某种极限过程去逐步逼近方程组精确解的方法. 迭代法具有计算机的存储单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛条件和收敛速度问题.迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)的重要方法.,2. 迭代法,5.2 高斯消去法,本节介绍高斯消去法(逐次消去法)及消去法和矩阵三角分解之间的关系. 虽然高斯消去法是一种古老的求解线性方程组的方法(早在公元前250年我国就掌握了解方程组的消去法)，但由它改进、变形得到的选主元素消去法、三角分解法仍然是目前计算机上常用的有效方法.我们在中学学过消去法,高斯消去法就是它的标准化的、适合在计算机上自动计算的一种方法.,5.2.1 高斯消去法,设有线性方程组,或写成矩阵形式,简记为Ax=b.,如: 上三角矩阵所对应的线性方程组,解为,当m=n时，对三角形方程组的解非常容易，如,下三角矩阵所对应的线性方程组,计算量（乘除法的主要部分）都为 n2/2.,解为,因此，我们将一般的线性方程组化成等价的三角形方程组来求解.,首先举一个简单的例子来说明消去法的基本思想.,例1 用消去法解方程组,解 第1步，将方程(2.2)乘上-2加到方程(2.4)上去，消去(2.4)中的未知数x1，得到,第2步，将方程(2.3) 加到方程(2.5)上去，消去(2.5)中的未知数x2，得到与原方程组等价的三角形方程组,显然，方程组是(2.6)是容易求解的，解为,上述过程相当于对方程的增广阵做初等行变换,其中ri用表示矩阵的第i行.,由此看出，用消去法解方程组的基本思想是用逐次消去未知数的方法把原方程组Ax=b化为与其等价的三角形方程组，而求解三角形方程组可用回代的方法求解. 换句话说，上述过程就是用初等行变换将原方程组系数矩阵化为简单形式(上三角矩阵)，从而求解原方程组(2.1)的问题转化为求解简单方程组的问题. 或者说，对系数矩阵A施行一些行变换(用一些简单矩阵左乘A)将其约化为上三角矩阵. 这就是高斯消去法.,下面讨论求解一般线性方程组的高斯消去法.由,将(2.1)记为A(1)x=b(1)，其中,(1) 第1步(k=1).,设a11(1)0，首先计算乘数 mi1= ai1(1) /a11(1) (i=2,3,m). 用-mi1乘(2.1)的第一个方程，加到第i个(i=2,3, ,m)方程上，消去(2.1)的从第二个方程到第m个方程中的未知数x1，得到与(2.1)等价的方程组,(2.7),简记为 A(2)x=b(2)，,其中A(2), b(2)的元素计算公式为,(2) 第k次消元(k=1,2,s, s=min(m-1,n).,设上述第1步， ，第k-1步消元过程计算已经完成，即已计算好与(2.1)等价的方程组，简记为 A(k)x=b(k)，其中,设akk(k)0，计算乘数 mik= aik(k) /akk(k) (i=k+1, ,m). 用-mik乘(2.8)的第k个方程加到第 i个(i= k+1, , m)方程上，消去从第k+1个方程到第m个方程中的未知数xk，得到与(2.1)等价的方程组A(k+1)x=b(k+1)，,其中A(k+1), b(k+1)的元素计算公式为，,显然A(k+1)中从第1行到第k行与A(k)相同.,(3) 继续上述过程，且设aii(i)0 (i=1,2, ,s)，直到完成第s步消元计算. 最后得到与原方程组等价的简单方程组A(s+1)x=b(s+1) ，其中A(s+1)为上阶梯.,特别当m=n时，与原方程组等价的简单方程组为A(n)x=b(n)，即,由(2.1)约化为(2.10)的过程称为消元过程.,如果A是非奇异矩阵，且akk(k)0 (k=1,2,n-1)，求解三角形方程组(2.10)，得到求解公式,(2.10)的求解过程(2.11) 称为回代过程.,注意：设Ax=b，其中ARnn为非奇异矩阵，如果a11(1)=0，由于A为非奇异矩阵，所以A的第1列一定有元素不等于零，例如al10，于是可交换两行元素(即r1rl)，将al1 调到(1,1)位置，然后进行消元计算，这时A(2)右下角矩阵为n-1阶非奇异矩阵，继续这过程，高斯消去法照样可进行计算.,总结上述讨论即有,定理5 设Ax=b，其中ARnn.,(1) 如果akk(k)0 (k=1,2,n)，则可通过高斯消去法将Ax=b约化为等价的三角形方程组(2.10)，且计算公式为,(a) 消元计算(k=1,2, ,n-1),(b) 回代计算,(2) 如果A为非奇异矩阵，则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为等价的三角形方程组(2.10).,例子 解方程组,解：消元,回代得,5.3 高斯主元素消元法,由高斯消去法知道, 在消元过程中可能有akk(k)0的情况，这时消去法将无法进行；即使主元素akk(k)0但很小时，用其作除数，会导致其它元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠.,高斯消去法也称主元素消去法 (条件det A0),即当akk(k)=0 时，高斯消元法无法进行；或 | akk(k) |1时，带来舍入误差的扩散。（小主元）,5.3.1 列主元素消去法,设方程组(2.1)的增广矩阵为,首先在A的第1列中选取绝对值最大的元素作为主元素，例如,交换,消元法的应用,1. 求行列式,高斯消元法,2. 求逆矩阵(用高斯-若当消去法),例子 分别用高斯消元法和列主元消元法求矩阵的行列式.,解： 高 斯 消 元 法,大数“吃掉”了小数！,列 主 元 消 元 法,5.3.2 高斯-若当消去法,高斯-若当消去法是高斯消去法的一种变形. 高斯消去法是消去对角元素下方的元素. 如果同时消去对角元素上方和下方的元素，而且将对角元素化为1，这种方法就称为高斯-若当(Gauss-Jordan)消去法.,设高斯-若当消去法已完成k-1步，于是Ax=b化为等价方程组A(k)x=b(k)，其中增广阵,在第k步计算时(k=1,2,n),考虑将第k行上、下的第k列元素都进行消元计算化为零，且akk化为1.,1. 按列选主元素，即确定ik使,2. 换行(当ikk时)，交换第k行与第ik行元素,(j=k,k+1,n),3. 计算乘数 mik=-aik/akk(i=1,2,n, ik) mkk=1/akk .,(mik可保存在存放aik的单元中).,4. 消元计算 aijaij+mikakj(i=1,2,n,且ik,j=k+1,n) bibi+mikbk (i=1,2,n,且ik),5.计算主行 akjakjmkk (j=k,k+1,n)， bkbkmkk .,上述过程结束后，即当k=n时有,显然xi=bi,i=1,2,n 就是 Ax=b的解.,说明用高斯-若当消去法虽比高斯消去法略复杂些，但将A约化为单位矩阵，计算解就在常数项位置得到，因此用不着回代求解，故也称为无回代的高斯消元法. 它的计算量大约需要n3/2次乘除法，要比高斯消去法大，但用高斯-若当消去法求一个矩阵的逆矩阵还是比较合适的.,定理9(高斯-若当法求逆矩阵) 设A为非奇异矩阵，方程组AX=In的增广矩阵为C=(A|In)，如果对C应用高斯-若当方法化为(In|T) ，则A-1=T.,事实上, 求A的逆矩阵A-1, 即求n阶矩阵X使AX=In, 其中In为单位矩阵，将X按列分块 X=(x1,x2,xn)， I=(e1,e2,en)， 于是求解AX=In等价于求解n个方程组 Axj=ej (j=1,2,n). 我们可用高斯-若当方法求解AX=In.,例4 用高斯-若当方法求下列矩阵的逆矩阵.,解 由分块矩阵,所以得,小节,比较而言,Gauss顺序消去法条件苛刻,且数值不稳定; Gauss全主元消去法工作量偏大，需要比较 个元素及行列交换工作，算法复杂；对于Gauss-Jordan消去法形式上比其他消元法简单，且无回代求解，但计算量大,比Gauss顺序消去法多 计算量。因此从算法优化的角度考虑， Gauss列主元消去法比较好。,5.4 矩阵的三角分解法,我们知道对矩阵进行一次初等变换，就相当于用相应的初等矩阵去左乘原来的矩阵。因此我们这个观点来考察Gauss消元法用矩阵乘法来表示，即可得到求解线性方程组的另一种直接法：矩阵的三角分解。,5.4.1 Gauss消元法的矩阵形式,5.4.2 Doolittle分解,Doolittle分解,若矩阵A有分解：A=LU，其中L为单位下三角阵，U为上三角阵，则称该分解为Doolittle分解，可以证明，当A的各阶顺序主子式均不为零时，Doolittle分解可以实现并且唯一。,A的各阶顺序主子式均不为零，即,Doolittle分解,Doolittle分解,Doolittle分解(看看就可),Doolittle分解,Doolittle分解,Doolittle分解,例题,例题,例题,例题,例题,上机,P27 例题,第六章 解线性方程组的迭代方法,6.1 引言 6.2 基本迭代法 6.3 迭代法的收敛性,例1 求解方程组,记为Ax=b，其中,方程组的精确解是x*=(3,2,1)T. 现将改写为,或写为x=B0x+f，其中,我们任取初始值，例如取x(0)=(0, 0, 0)T. 将这些值代入(1.3)式右边(若(1.3)式为等式即求得方程组的解，但一般不满足)，得到新的值 x(1)=(x1(1), x2(1), x3(1)T =(3.5, 3, 3)T ，再将x(1)分量代入(1.3)式右边得到 x(2)，反复利用这个计算程序，得到一向量序列和一般的计算公式(迭代公式),简写为 x(k+1)=B0x(k) +f，,其中k表示迭代次数(k=0,1,2,).,迭代到第10次有,从此例看出，由迭代法产生的向量序列x(k)逐步逼近方程组的精确解是x*=(3,2,1)T. 即有,对于任何一个方程组x=Bx+f(由Ax=b变形得到的等价方程组)，由迭代法产生的向量序列x(k)是否一定逐步逼近方程组的解x*呢？回答是不一定. 请同学们考虑用迭代法解下述方程组,但 x(k)并不是所有的都收敛到解x*！,6.2 基本迭代法,其中，A=(aij)Rnn为非奇异矩阵，下面研究任何建立Ax=b的各种迭代法.,设线性方程组,Ax=b， (2.1),其中，M为可选择的非奇异矩阵，且使Mx=d容易求解，一般选择A的某种近似，称M为分裂矩阵.,将A分裂为,A=M-N. (2.2),于是，求解Ax=b转化为求解Mx=Nx+b ，即求解,可构造一阶定常迭代法,其中 B=M-1N=M-1(M-A)=I-M-1A ， f=M-1b. 称 B=I-M-1A为迭代法的迭代矩阵，选取M矩阵，就得到解Ax=b的各种迭代法.,设aii0 (i=1,2,n)，并将A写成三部分,即 A=D-L-U,6.2.1 雅可比(Jacobi)迭代法,设aii0 (i=1,2,n)，选取M为A的对角元素部分,即选取M=D(对角阵)，A=D-N，由(2.3)式得到解方程组Ax=b的雅可比(Jacobi)迭代法. 又称简单迭代法.,其中B=I-D-1A=D-1(L+U)=J, f=D-1b. 称J为解Ax=b的雅可比迭代法的迭代矩阵.,于是雅可比迭代法可写为矩阵形式,其Jacobi迭代矩阵为,等 价 方 程 组,其中 aii(i)0 (i=1,2,n),即由方程组Ax=b得到的,建立的雅可比迭代格式为,于是，解Ax=b的雅可比迭代法的计算公式为,由(2.6)式可知，雅可比迭代法计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法且计算过程中原始矩阵A始终不变.,6.2.2 高斯赛德尔迭代法,在 Jacobi 迭代中，计算 xi(k+1)(2 i n)时，使用xj(k+1)代替xj(k) (1 j i-1)，即有,建 立 迭 代 格 式,或缩写为,称为高斯塞德尔(Gauss Seidel)迭代法.,解Ax=b的高斯塞德尔迭代法的计算公式为,或,雅可比迭代法不使用变量的最新信息计算xi(k+1)，而由高斯塞德尔迭代公式(2.8)可知，计算x(k+1)的第 i个分量xi(k+1)时，利用了已经计算出的最新分量xj(k+1) (j=1,2,i-1). 可看作雅可比迭代法的一种改进. 由(2.8)可知，高斯塞德尔迭代公式每迭代一次只需计算一次矩阵与向量的乘法.,例2 用高斯塞德尔迭代法解例1的方程组(1.2).,解 用高斯塞德尔迭代公式：取x(0)=(