第六章-方程组的数值解法

《计算方法》方程组的数值解法6★1、高斯消去法☆2、选主元素的高斯消去法☆3、矩阵的三角分解※4、方程组的迭代解法2008年10月6日12时7分沈阳航空工业学院飞机设计教研室主要知识点第2页，共42页引言（一）2008年10月6日12时7分沈阳航空工业学院飞机设计教研室快速、高效地求解线性方程组是数值线性代数研究中的核心问题，也是目前科学计算中的重大研究课题之一。各种各样的科学和工程问题，往往最终都要归结为求解一个线性方程组。线性方程组的数值解法有：直接法和迭代法。直接法：在假定没有舍入误差的情况下，经过有限次运算可以求得方程组的精确解；迭代法：从一个初始向量出发，按照一定的迭代格式，构造出一个趋向于真解的无穷序列。第3页，共42页在没有舍入误差的情况下，经过有限次运算可以得到方程组的精确解的方法。求解,nnAxbAR0det()ACramer法则:12,,,iiDxinD所需乘除法的运算量大约为(n+1)!+nn=20时，每秒1亿次运算速度的计算机要算30多万年！直接法引言（二）2008年10月6日12时7分沈阳航空工业学院飞机设计教研室第4页，共42页举例（一）2008年10月6日12时7分沈阳航空工业学院飞机设计教研室解：例：直接法解线性方程组1231231232222334463xxxxxxxxx1222(,)23344163Ab61121702128006178921122200131x23871xx1232222xxx第5页，共42页高斯消去法2008年10月6日12时7分沈阳航空工业学院飞机设计教研室高斯消去法的主要思路：将系数矩阵A化为上三角矩阵，然后回代求解。11112211211222221122.........nnnnnnnnnnaxaxaxbaxaxaxbaxaxaxb考虑n阶线性方程组：Axb矩阵形式=通过初等变换逐步消去未知元，将原方程组化为同解的三角方程组。第6页，共42页第一步：消去第一列依次将增广矩阵的第i行+mi1第1行，得)1(1)1(1)1(12)1(11...baaan)2(A(2)b记(1)1(1)(1)(1)(1)(),ijnnnbAaAbbb，即。(1)(1),ijijiiaabb(1)110a(1)(1)1111(2,...,)iiianma设，计算(2)(1)(1)11(2)(1)(1)11ijijijiiiaamabbmb其中(,2,...,)ijn第二步：消去第二列依次将上述矩阵的第i行+mi2第2行，得(3)(2)(2)22(3)(2)(2)22ijijijiiiaamabbmb其中(,3,...,)ijn(2)(2)2222(3,...,)iiianma设，计算(2)220a)1(1)1(1)1(12)1(11...baaan(3)A(3)b(2)(2)(2)2222...naab2008年10月6日12时7分沈阳航空工业学院飞机设计教研室第7页，共42页高斯消去法2008年10月6日12时7分沈阳航空工业学院飞机设计教研室第k步：消去第k列依此类推，直到第n-1步，原方程化为()()(1,...,)kkikikkkmakain设，计算()0kkka(1)(1)(1)(1)1112111(2)(2)(2)22222()()nnnnnnnnaaabxxaabxba回代求解：()()nnnnnnxba计算(1)()()(1)()()kkkijijikkjkkkiiikkaamabbmb(i=k+1,…,n)()()()1()niiiiiijjiijixbaxa(i=n-1,…,1)第8页，共42页例：用基本Gauss消元法求解下列方程组123123123223347712457xxxxxxxxx解：增广矩阵223347712457()Ab2233031506842008年10月6日12时7分沈阳航空工业学院飞机设计教研室第9页，共42页223347712457()Ab223303150684223347712457121121112()()ala131131111()()ala232232222()()ala223323151684223323151266321122,,xxx乘除运算量2008年10月6日12时7分沈阳航空工业学院飞机设计教研室乘除运算量：由于计算机中做乘除运算的时间远远超过做加减运算时间，故估计运算量时，往往只估计乘除的次数。第k步：消去第k列()()(1,...,)kkikikkkmakain设，计算()0kkka计算(1)()()(1)()()kkkijijikkjkkkiiikkaamabbmb(i=k+1,…,n)回代求解：()()nnnnnnxba()()()1()niiiiiijjiijixbaxa(i=k+1,…,n)n–k次(n–k)2次n–k次n(n+1)/2次第11页，共42页基本Gauss消元法的工作量消元过程：11111()()()nnkknknknk回代过程：12()nn325326nnn332333nnnn加减法的次数32353263nnnn乘除法的次数3()On2008年10月6日12时7分沈阳航空工业学院飞机设计教研室第12页，共42页列主元高斯消去法2008年10月6日12时7分沈阳航空工业学院飞机设计教研室高斯消去法有效的条件是：()0kkka主元如果某个主元为0，则可能导致高斯消去法求解失败。列主元高斯消去法：在第k步消元时①先选取列主元：()||=kkika()max||0kikkina②ifikkthen交换第k行和第ik行；列主元高斯消去法比普通高斯消去法要多一些比较运算，但比普通高斯消去法稳定。③消元第13页，共42页全主元高斯消去法2008年10月6日12时7分沈阳航空工业学院飞机设计教研室全主元高斯消去法：第k步消元时选A(k)中绝对值最大的元素为主元，即①先选取全主元：()||=kkkija(),max||0kijkijna②ifikkthen交换第k行和第ik行；ifjkkthen交换第k列和第jk列；③消元列交换改变了xi的顺序，须记录交换次序，解完后再换回来。全主元高斯消去法具有很好的稳定性，但选全主元比较费时，故在实际计算中很少使用。第14页，共42页举例（二）2008年10月6日12时7分沈阳航空工业学院飞机设计教研室解：例：采用十进制八位浮点数，分别用Gauss消去法和列主元Gauss消去法求解线性方程组:121212xxxx9(10)精确解为，...1000...00.1101191x8个...8999...99.0212xx8个Gauss消去法：911212110/aam1010102221111100.100.0...00.110...010...00am9个101010221212100.100.0...00.110...010...00bm999101101010211,0xx小主元可能导致计算失败。第15页，共42页举例（二）续2008年10月6日12时7分沈阳航空工业学院飞机设计教研室列主元Gauss消去法：21111109112011911210111211xx例：但列主元高斯消去法有时也会导致求解失败。9911010112999911010010101201xx第16页，共42页矩阵三角分解2008年10月6日12时7分沈阳航空工业学院飞机设计教研室将一个矩阵分解成结构简单的三角形矩阵的乘积称为矩阵的三角分解。高斯消去过程其实就是一个矩阵的三角分解过程。第17页，共42页则A(k)与A(k+1)之间的关系式可以表示为：(1)()kkkALA其中：1,,1111kknkkmmL()()kkikikkkmaa(i=k+1,…,n),将Gauss消去过程中第k-1步消元后的系数矩阵记为：(1)(1)(1)1111()()()()()knkkkkkknkknknnaaaAaaaa(k=1,…,n-1)2008年10月6日12时7分沈阳航空工业学院飞机设计教研室第18页，共42页LU分解2008年10月6日12时7分沈阳航空工业学院飞机设计教研室记：，则111()12,nkLLLLUAALU其中：L---单位下三角矩阵，U---上三角矩阵LU分解于是有：()(1)121nnALLLA(1)()1121()nnALLLAA容易验证：1,,11111kknkkmmL(k=1,…,n-1)(杜利脱尔Doolittle分解)第19页，共42页LU分解的存在唯一性2008年10月6日12时7分沈阳航空工业学院飞机设计教研室LU分解存在()0kkka普通高斯消去法不被中断?定理普通高斯消去法求解方程组Ax=b时，主元的充要条件是：A的所有顺序主子式不为零。()0kkka定理若A的所有顺序主子式，则A存在唯一的LU分解。0(LU分解的唯一性)证：存在性由上面的定理可得；唯一性可用反证法证明。第20页，共42页类似分解2008年10月6日12时7分沈阳航空工业学院飞机设计教研室定理若A的所有顺序主子式，则0(1)A存在唯一的LDR分解：A=LDR，其中：L是单位下三角矩阵，D是对角矩阵，R是单位上三角矩阵ALU(2)A存在唯一的克洛脱(Crout)分解：，其中：L是下三角矩阵，U是单位上三角矩阵1diag(),DURDU1,LLDUDU第21页，共42页LU分解分解。行直接进否对矩阵因此，关键问题在于能个三角方程：就等价于解两由此解线性方程组LUAyUxbLybUxLbAx)(1112131112132122232122233132333132331111212111212111313111313111,31111(1,2,3)jjjjLUnaaauuuaaaluuaaallukauuajaaulluaaullu此分解在于如何算出的各元素，以为例时：由得；由得LU分解)(322332133133333323321331332212313232233212313213212323231321231221222222122122ululauuululakuulalululaulauuulaulauuulak得时：由得由；得由；得时：LU分解紧凑方式2008年10月6日12时7分沈阳航空工业学院飞机设计教研室直接利用矩阵乘法来计算LU分解111211112121222212221,112111nnnnnnnnnnnnnuuuaaaluuaaalluaaaLUA比较等式两边的第一行得：u1