第2章线性方程组求解的数值方法优秀课件.ppt

第2章线性方程组求解的数值方法第1页，本讲稿共85页第二章 线性方程组求解的数值方法1.高斯消元法2.矩阵分解法3.向量范数与矩阵范数4.迭代法求解5.方程组的病态问题与误差分析主要内容：第2页，本讲稿共85页第二章 线性方程组求解的数值方法1.理解各种线性方程组数值求解；2.掌握求解方法和解的误差分析方法；3.能编程实现求解算法。特别强调：遇到问题养成用计算机编程求解的习惯，不要习惯性的用笔算，而这是国内外学生的一个主要差距。教学要求：第3页，本讲稿共85页 在自然科学和工程技术中，有很多问题的解决都需要用到线性方程组的求解。因此，求解线性方程组的问题是一个在科学技术中常见的普遍问题。解线性方程组的数值解法：有直接法和迭代法两类。直接法：计算过程没有舍入误差，经过有限次四则运算可求得方程组得精确解。（实际计算有舍入误差）高斯消元法，矩阵分解法迭代法：核心是迭代求解的收敛条件和收敛速度。雅可比（Jacobi）迭代，高斯-赛德尔（Gauss-Seidel）迭代第4页，本讲稿共85页基本思想方法：由行初等变换将系数矩阵约化为三角 矩阵；用回代的方法求解方程组。例1 用消去法（高斯消元法）解方程组高斯消元法是求解方程组的古典方法。(2.1)3.1 3.1 高斯消元法高斯消元法第5页，本讲稿共85页结论：整个计算过程可分为两部分：（1）消元：把原方程组转化为系数矩阵为上三角矩阵的方程组；（2）回代：由系数矩阵为上三角矩阵的方程组求解（2）回代求解，得：解：（1）消元：第6页，本讲稿共85页对于一般情形：n阶线性方程组的高斯消元法第7页，本讲稿共85页若记增广矩阵（2.2）第8页，本讲稿共85页（1）第1步(k=1)，一般形式的高斯消元法：设 ，首先对行计算乘数对增广矩阵 进行行初等变换：(即用 乘以2.2式的第1个方程，加到第i个方程上，消去2.2式的第二个方程直到第n个方程中的未知数 )(代表第i行)第9页，本讲稿共85页得到与原方程组 等价的方程组 。记为（2）一般第k步消元()设已完成上述消元过程第设已完成上述消元过程第1 1步，第步，第2 2步，步，第，第k k-1-1步，且步，且 第10页，本讲稿共85页其中：设 ，计算乘数第11页，本讲稿共85页(即用 乘以2.2式的第k个方程，加到第i个方程上，消去2.2式的第k+1个方程直到第n个方程中的未知数 )那么第k步消元操作即：(3)继续这一过程，直到完成第n-1次消元，最后我们得到与原方程组等价的三角形方程组(2.3)消元过程结束第12页，本讲稿共85页求解三角形方程组2.3，得到求解公式这个过程称为回代过程。第13页，本讲稿共85页例题：考虑方程组 Gauss消去法中每步用来消去其他元素的 称为该步的主元素。Gauss消去法作为数值方法，主元素的选择是否会影响计算的结果呢？则该方程的精确解为而采用(，1)作为主元素，利用高斯消去法得到的解为显然结果是错误的。第14页，本讲稿共85页 错误在哪个地方呢？原因是我们在消元时，利用了小主元 ，使得约化后的方程组元素数量级大大增长，再经舍入，而计算机的有效位数有限，造成消元后的三角形方程组就不准确了。结论：在消元过程中可能出现主元素为零的情况，这时消去法将无法进行；即使不为零，在主元素很小时，用其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠。解决方法：对一般的矩阵来说，最好每一步选取系数矩阵（或消元后的低阶矩阵）的该列中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数字稳定性。(高斯列主元素消去法)第15页，本讲稿共85页1.列主元法第一列中绝对值最大是10，取10为主元n n阶方程组第阶方程组第k k轮消元时，选第轮消元时，选第k k列的后列的后(n-k+1)(n-k+1)个元素中绝对个元素中绝对值最大作主元。值最大作主元。第16页，本讲稿共85页x3=6.2/6.2=1x2=(2.5-5x3)/2.5=-1x1=(7+7x2-0 x3)/10=0 x1=0 x2=-1x3=1第二列的后两个数中选出主元第二列的后两个数中选出主元 2.52.5第17页，本讲稿共85页2 完全主元素消去法整个矩阵中绝对值最大是10，取10为主元n n阶方程组第阶方程组第k k轮消元时，选消元后元素中绝对值最大作主元。轮消元时，选消元后元素中绝对值最大作主元。第18页，本讲稿共85页x1=0 x2=-1x3=1方框中方框中6 6最大，交换行列，交换列的时候要做记录最大，交换行列，交换列的时候要做记录（即（即x3x3和和x2x2交换了位置）：交换了位置）：第19页，本讲稿共85页 完全主元素消除法与列主元素消除法的优缺点比较：优点：数值更加稳定；缺点：计算量大；第20页，本讲稿共85页 对矩阵 A实行初等行变换相当于用初等矩阵左乘 A，于是对(2.2)做第一次消元后，化 为 化 为 ，即 ，其中 3.1 矩阵的三角分解矩阵的三角分解 LULU分解分解第21页，本讲稿共85页第k步的初等矩阵为并且重复这一过程，最后得到将上三角矩阵 记为U，则 第22页，本讲稿共85页将上三角矩阵 记为U，则 ，其中则，L为单位下三角矩阵。高斯消去法实质上是产生了一个将A分解为两个三角形矩阵相乘的因式分解。如果A是非奇异阵，则LU分解是唯一的，否则分解不唯一。第23页，本讲稿共85页消元法：消元法与三角分解法间的关系：三角分解法：讨论第24页，本讲稿共85页直接三角分解法解线性方程组（）的具体流程：1.2.2.计算计算U U的第的第r r行，行，L L的第的第r r列元素列元素 r=2,3,n3.(一)LU分解第25页，本讲稿共85页再求解Ly=b,Ux=y计算公式：(二)x的计算第26页，本讲稿共85页例 用直接三角分解法解方程组 解:(一)矩阵LU分解(1)故：第27页，本讲稿共85页(2)经计算：第28页，本讲稿共85页(二)求解x：从而第29页，本讲稿共85页 3.2 3.2 矩阵的矩阵的CholeskyCholesky分解分解 对称正定矩阵A：AT=A,A的各阶顺序主子式大于零.前面指出非奇异的矩阵可以有三角分解，当 A是某些特殊矩 阵 时，它 的 L U 分 解 会 有 更 加 简 洁 的 形 式。A的LU分解(2.4)第30页，本讲稿共85页uii 0(i=1,n)由于A是对称正定的，则有将U进一步分解：第31页，本讲稿共85页根据A的对称性：故：由分解的唯一性知：故：其中P为上三角矩阵，这种对称正定矩阵的分解称为Cholesky分解。在在MatlabMatlab中函数中函数“cholchol”给出对称正定矩阵的给出对称正定矩阵的CholeskyCholesky分解。分解。第32页，本讲稿共85页经过n步可直接求得思路:(1)分解对称正定矩阵设n阶对称正定矩阵A有分解 ,先用待定系数法求L的元素Cholesky分解的具体步骤(平方根法)：(2)分解求解方程组第33页，本讲稿共85页Cholesky方法具体计算公式：分解计算：求解：第34页，本讲稿共85页 3.3 3.3 向量范数与矩向量范数与矩阵范数范数向量、矩阵与线性方程组有着密切的关系，向量、矩阵范数是解方程组以及研究与探讨方程组本身性质的工具。一、向量范数的定义定义 范数为其中的 ，最常用的范数是 ，以及第35页，本讲稿共85页容易证明向量范数满足如下性质：(1)如果 ，则 ；而且 ，当且仅当 (2)对任何的数c，都有 (3)范数也称为2-范数或Euclidean函数；范数也称为 -范数或一致范数。在Matlab中用函数 用来计算向量x的 范数。由性质(3),还容易得到：第36页，本讲稿共85页二 矩阵范数的定义在向量范数的基础上可以进一步定义矩阵范数 如果上式中的向量范数取为 范数，则可以证明定义出的矩阵范数（称为矩阵 范数或者列范数）为 如果向量范数取为2-范数，则 其中 （模），称为矩阵B的谱半径。(为B的任意特征值。)第37页，本讲稿共85页 类似地，Matlab函数norm(A,p)给出矩阵p=1,2或 范数 如果向量范数取为 ，则可以证明定义出的矩阵范数（称为矩阵的 范数或者行范数）为 容易证明矩阵范数满足如下性质：(1)如果 ，则 ；而且 ，而且仅当 (2)对任何的数 ，(3)(4)而且由矩阵范数的定义还有称之为矩阵范数与相应的向量范数是相容的。第38页，本讲稿共85页39 3.4 古典迭代法的构造古典迭代法的构造求解线性代数方程组的方法求解线性代数方程组的方法中小规模中小规模问题问题直接法直接法迭代法迭代法 大规模，大规模，超大规模问题超大规模问题 古典方法古典方法现代方法现代方法第39页，本讲稿共85页40线性方程组的一般形式为：如果 是非奇异的，则上式有唯一解。我们将通过一个具体线性方程组的例子来讲解迭代法。取：即线性方程组为：方程的精确解 (直接法计算)第40页，本讲稿共85页41 如果对该方程组进行变形，将变量 分别从三个方程中分离出来：(1)令初值第41页，本讲稿共85页所以(1)式可以表示为迭代形式：所以我们可以得到一个向量的序列 ，只要该序列当 时有极限，那么这个极限就是该线性方程组的解。上面这种迭代求解线性方程组的方法称为Jacobi迭代法。第42页，本讲稿共85页考虑线性方程组的一般形式：其中矩阵其中矩阵A为为nn阶方阵，则方程的分量形式为：阶方阵，则方程的分量形式为：改写成：改写成：第43页，本讲稿共85页从而得到迭代公式：上面式子就是一般形式的Jacobi迭代法，也可也记做：第44页，本讲稿共85页在Jacobi迭代法中充分利用新值，可得到下面迭代：上面这种迭代方法也叫做Gauss-Seidel迭代，也可以记为：第45页，本讲稿共85页方程组的精确解为方程组的精确解为x x*=(1,1,1)=(1,1,1)T T.解解 JacobiJacobi迭代法计算公式为迭代法计算公式为取初始向量取初始向量x x(0)(0)=(0,0,0)=(0,0,0)T T,迭代可得迭代可得例1：利用Jacobi法和Gauss-Seidel法求解线性方程组第46页，本讲稿共85页kx1(k)x2(k)x3(k)x(k)-x*0123456701.41.110.9290.99061.011591.0002510.998236400.51.201.0550.96450.99531.0057951.000125501.41.110.9290.99061.011591.0002510.998236410.50.20.0710.03550.011590.0057950.0017636可见可见,迭代序列逐次收敛于方程组的解迭代序列逐次收敛于方程组的解,而且迭代而且迭代7 7次得到精确到小次得到精确到小数点后两位的近似解数点后两位的近似解.计算结果列表如下:第47页，本讲稿共85页Gauss-Seidel迭代法的计算公式为:同样取初始向量同样取初始向量x x(0)(0)=(0,0,0)=(0,0,0)T T,计算结果为计算结果为 由计算结果可见由计算结果可见,Gauss-Seidel,Gauss-Seidel迭代法收敛较快迭代法收敛较快.取精确到小数取精确到小数点后两位的近似解点后两位的近似解,Gauss-Seidel,Gauss-Seidel迭代法只需迭代迭代法只需迭代3 3次次,而而JacobiJacobi迭迭代法需要迭代代法需要迭代7 7次次.kx1(k)x2(k)x3(k)x(k)-x*012301.41.06340.995104400.781.020480.9952756801.0260.9875161.0019068610.40.06340.0048956 第48页，本讲稿共85页为了进一步研究,从矩阵角度来讨论上述迭代法.对线性方程组对线性方程组Ax=b,Ax=b,记记D=diag(aD=diag(a1111,a,a2222,a,annnn)则有则有 A=D-L-UA=D-L-U于是线性方程组于是线性方程组 Ax=b Ax=b 可写成可写成 (D-L-U)x=b(D-L-U)x=b等价于等价于 Dx=(L+U)x+b Dx=(L+U)x+b 或或 x=Dx=D-1-1(L+U)x+D(L+U)x+D-1-1b b 第49页，本讲稿共85页由此建立由此建立JacobiJacobi迭代法迭代公式迭代法迭代公式 x x(k+1)(k+1)=D=D-1-1(L+U)x(L+U)x(k)(k)+D+D-1-1b k=0,1,2,b k=0,1,2,或写成或写成 x x(k+1)(k+1)=Bx=Bx(k)(k)+f k=0,1,2,+f k=0,1,2,其中其中第50页，本讲稿共85页所以Gauss-Seidel迭代法可以写成其中Gauss-Seidel迭代法迭代公式为写成矩阵形式为：第51页，本讲稿共85页 3.5 3.5 迭代法的分析迭代法的分析 构造迭代法的途径可以有很多，那么是不是所有的方法都能收敛呢？收敛速度的快慢与什么有关系呢？它的精确解为例：第52页，本讲稿共85页kx1(k)x2(k)x3(k)x(k)-x*0123401416.792.8168104.7-41.4-52.4-254.30-1.74.56-151-238.211342.41521680可见,并不是所有的方程组都收敛，即使收敛对于不同的方法(例如Jacobi与Gauss-Seidel)收敛速度也是不同的.计算结果列表如下:第53页，本讲稿共85页 设 是矩阵B任意一个特征值，是相应的特征向量，即 再假设 为任意的向量范数及与之相融的矩阵范数，则有：首先引入几个定义、引理：所以对矩阵B任意一个特征值 ，都有记 ，称之为矩阵B的谱半径，则第54页，本讲稿共85页先引入两个引理（详细证明见附录）：引理3.1 矩阵 ，则 的充分必要条件是 引理3.2 矩阵 ，若 ，则 是非奇异的。第55页，本讲稿共85页 定理3.1 对任意的 和任意的初始向量 ，迭代法收敛的充分必要条件是证明：必要性 设迭代法产生的序列 收敛，记 是该序列的极限点，所以 满足又由迭代关系 ，可以得到第56页，本讲稿共85页由 的任意性，知：由引理3.1知：充分性 由 及引理2知 （即 ）有唯一解，记为 类似于必要性的证明，得到再由第一步知 ，故第57页，本讲稿共85页 定理3.1给出了迭代收敛的充要条件，但 不宜计算，所以在实际使用中通常并不好用。由 知，只要对某种相容的矩阵范数有 那么当然 ，所以这常常是实际中很有效的收敛判别准则。严格对角占优矩阵：考虑 ，设 ，当 时的矩阵。第58页，本讲稿共85页定理3.2 如果A是严格对角占优的矩阵，则它一定非奇异。由于A是对角占优矩阵，所以 ，故矩阵是可逆的。考虑第59页，本讲稿共85页 矩阵的 范数为由A是对角占优矩阵，得 由引理3.2知 是非奇异的，又由于 是非奇异的，所以A是非奇异的。得证！引理3.2 矩阵 ，若 ，则 是非奇异的。第60页，本讲稿共85页 定理3.3 系数矩阵为严格对角占优的线性代数方程组，它的Jacobi迭代法和Gauss-Seidel迭代都是收敛的。证明：Jacobi方法的迭代矩阵为由定理3.2中的证明知，严格对角占优矩阵满足：再由定理3.1，得Jacobi迭代法收敛。(a)Jocobi的证明1第61页，本讲稿共85页再考虑Gauss-Seidel方法，其迭代矩阵为用反证法，假设Gauss-Seidel迭代不收敛，则由定理3.1知：所以存在模大于1的特征值.设有 ，使得行列式：(b)Gauss-Seidel的证明第62页，本讲稿共85页 故，只有才能使得 成立。由于A是对角占优的，所以矩阵 也是对角占优的，那么该矩阵一定非奇异，这与上面该矩阵 的行列式等于0矛盾。由于 可逆，所以得证！第63页，本讲稿共85页定理3.4 若迭代矩阵B满足 ，则迭代生成的序列 满足收敛速度问题其中 表示精确解。证明：先证明(b),然后证明(a)第64页，本讲稿共85页而所以故，本页第一个式子可写为第65页，本讲稿共85页由于 ，故(b)得证显然对于任意的正整数p，都有那么在本页第一个式子中反复利用上式就可以得到(a)(a)得证第66页，本讲稿共85页 给出的关系可以作为计算终止的判别准则。即只要 和 已足够接近，表明 与 便已足够靠近。定理3.4(a)给出的是迭代收敛速度的一个估计。显然 越接近0，生成序列 收敛得越快。定理3.14(b)定理3.4的讨论第67页，本讲稿共85页 3.6 超松弛迭代超松弛迭代(SOR)及分块迭代方法及分块迭代方法 对于给定的迭代法，每步迭代所需的工作量是确定的。如果迭代法收敛速度缓慢，则需要比较多的迭代次数，由此导致算法工作量太大而失去使用价值，因此各种迭代法的加速技术具有重要意义。假设 是已经得到的迭代值，以 为初值进行一步Gauss-Seidel迭代得：第68页，本讲稿共85页将这一迭代值与 的值组合作为新一步的迭代值其中 为待定的参数，写成矩阵形式为：这时迭代矩阵：第69页，本讲稿共85页 新的迭代每步的计算代价与Gauss-Seidel方法相差无几，如果能取得恰当的 值，使新构造的方法比Gauss-Seidel方法收敛得更快，松弛就起到了加速的作用。在实际上真正使用的 值通常的范围为 ，被统称为超松弛方法，简称SOR(Successive Over-Relaxation)方法。上面的这个方法就叫做松弛方法，它可以视为Gauss-Seidel方法的加速。显然 时，上面这个方法就是Gauss-Seidel方法。第70页，本讲稿共85页定理3.5 SOR方法收敛的必要条件是：证明：假设SOR方法收敛，所以设 为 的特征值，则另一方面，经过计算不难得到故 ，结论得证。第71页，本讲稿共85页分块矩阵以及分块迭代方法在许多实际应用中，矩阵可以写成如下分块的形式：其中 是 阶的方阵，而且 据此可以构造分块形式的迭代法，考虑 其中第72页，本讲稿共85页求解方程组 ，分块的Jacobi迭代法为：分块的Gauss-Seidel迭代法为：第73页，本讲稿共85页分块的超松弛迭代法(BSOR)为第74页，本讲稿共85页 3.7 线性方程组的条件数线性方程组的条件数 对于线性方程组对于线性方程组 如如果果解解x关关于于问问题题（即即矩矩阵阵A和和向向量量b）的的“微微小小”变变化化(来来源源可可能能是是模模型型误误差差或或数数据据上上的的误误差差，也也可可能能是是数数值值计计算算中中的的舍舍入入误误差差)不不敏敏感感，那那么么这这个个问问题题即即是是一一个个“好好”问题，反之就是问题，反之就是“坏坏”问题或者病态的问题。问题或者病态的问题。关于病态的例子，我们在第一章中已经包含了一个，这一章主要用范数作为工具来分析线性方程组的好坏。第75页，本讲稿共85页例 设有方程组已知系数测量数据待求数据其中A为 ，若测量数据为准确数据，即 ，那么很容易计算如果测量数据不够精确，即测量数据仅保留了2位有效数字，即 ，计算结果为第76页，本讲稿共85页如果测量数据和系数矩阵的数据都仅保留了2位有效数字，即，则计算结果为以上问题称为病态的问题，病态是问题本身固有的。以上问题称为病态的问题，病态是问题本身固有的。第77页，本讲稿共85页(1)考虑由于右端的扰动所引起的解的变化，已知两式相减得利用范数关系可以得到：综合上述两式有：第78页，本讲稿共85页 这个式子给出了右端的扰动可能引起解扰动的上界，显然 越小，方程右端的变化引起解的变化就越小。(2)考虑由于矩阵A有扰动时所引起的解的变化可以得到：其中利用了第79页，本讲稿共85页经过整理可以得到(3)考虑更一般的情形所引起的解的变化假设满足如下的特殊形式：又Tylor展开，可以写成第80页，本讲稿共85页对 求导得令 ，并考虑到 ，则记 ，从而有第81页，本讲稿共85页 从上面三种情况可见，问题扰动所引起解的扰动是被同一个因子 控制的 定义3.1 设A为可逆矩阵，是与某种向量范数相容的矩阵范数，称 为线性方程组 相应于该矩阵范数的条件数。由前面的分析可见，条件数大的就是病态的矩阵，反之是良态的矩阵。求解线性方程时了解它的条件数大小是必要的，它可以帮助你判断所得数值解的可信性及模型的合理性。在Matlab中可以用cond（A）来计算条件数或者用rcond(A)来估计其量级的倒数。第82页，本讲稿共85页设 ，则 ，即对于该函数，误差会被放大n倍。二、病态问题与条件数（针对问题本身）定义：输入数据的微小变动导致输出数据的较大误差，就被称为病态问题。衡量是否病态的标准：条件数 对于函数值计算问题，条件数定义为：不同的问题，条件数具体定义不同。一般情况下，条件数大于一般情况下，条件数大于1010，就认为问题病态。，就认为问题病态。通过构造特殊算法来解决第83页，本讲稿共85页条件数的重要性质：定理3.6 (a)对任意的非奇异矩阵A，cond(A)是由任意的矩阵 范数定义的条件数，则 (b)如果Q是正交矩阵，则对于任意的非奇异矩阵A，都有其中第84页，本讲稿共85页证明：(a)中的各条性质是显然的，例如，对任意的非奇异矩阵A，有 (b)由于Q是正交矩阵，所以 另外，第85页，本讲稿共85页