实际问题中解线性方程组的经典解法

第二章在实际问题中解线性方程组的经典解法直接法与三角形方程组的求解分析线性方程组求解问题在许多科学计算问题中都会遇到，如应力分析、电学网络、自由振动问题等。在计算机数值方法的课程中，2线性方程组求解在样条插值、数据拟合的最小二乘法以及常微分方程边值问题中都要用到.产生的线性方程组的类型有很多，如按系数矩阵含零元素多少分类，有稠密和稀疏（零元素占80%以上）线性方程组之别；如按阶数的高低分类，有高阶（阶数在1000阶以上）和低阶之别；如按系数矩阵的形状和性质分类，又有对称正定、三对角线对角占优等之别，因为在电子计算机上求解，必须要考虑算法的计算复杂性以及算法的数值稳定性问题。所以针对不同类型的线性方程组有不同的解法，但是，基本的方法可归结为两大类，即为直接法和迭代法。本章介绍的经典解法，都把原方程组化为一个或者两个三角形方程组来求解，主要包括Gauss4消去法和它的变形----直接三角分解法设有线性方程组bAx（1，1）其中nnnnnnaaaaaaaaaA212222111211nxxxx21nbbbb21根据线性代数知识，当0detA时，方程组（1.1）的解存在且唯一，对增广矩阵)(),()1()1(bAbA施行行初等变换，化)1(A为上三角形矩阵)(nA，同时)1(b化为)(nb，这时与增广矩阵),()()(nnbA相应的线性方程组为上三角形方程组)()(nnbxA（1.2）其中)()()2(2)2(2)2(22)1(1)1(1)1(12)1(11)()(000),(nnnnnnnnnbabaabaaabA)()2(2)2(22)(!1)(!12)(!11)(000)(nnnnnnaaaaaaA)()2(2)(!1)()(nnnbbbb设)3.2.1(0)(niaiii则（1.2）的解为（1.3）它便是原方程组（1.1）的解，实现上述求解过程的方法称为Gauss消去法如果方程组（1.1）的系数矩阵可分解为两个形式简单的三角形矩阵L和U的乘积，即LUA(1.4)若L为下三角形矩阵，则U为上三角形矩阵，反之亦然从而求解bAx的问题转化为解三角形方程组bLy(1.5)和yUx（1.6）其中nnnnllllllL21222111000nnnnuuuuuuU00022211211则bLy为下三角形方程组，它的第i个方程为)2.1(2211nibylylyliiiiii（1.7）假定0iil按nyyy21,的顺序解得ijijijijilbylylby111111)3,2(ni（1,8）上三角形方程组yUx的第i个方程为ininiiiiiiyxuxuxu11)2,1(ni（1.9）假定0iiu按11,,xxxnn的顺序求解得iinijijijinnnnuyxuyuyx1（1,10）直接法的求解过程，在计算过程中无舍入误差的前提下，都可以经有限步算数运算而得到精确解。由于计算机的字长有限，初始数据取浮点数以及运算过程都不断地产生舍入误差，这些误差的传播和积累，会影响计算精度，所以，如何避免舍入误差的增长是设计算法时必须考虑的问题。直接法通常需要存储系数矩阵的全部元素，当方程组的阶数很高时，需要相当大的内存空间，因此，在算法设计上应当注意节省内存，比如，对称矩阵可以只存其下三角形部分于一维数组中。综上所述，如果矩阵A非奇异，总可以通过带有行交换或不带行交换的消元过程，将A化为非奇异上三角形矩阵)(nA因此，回代求解过程（1,3）也可以进行到底，但是，在实际应用中，常常难于事先判断系数矩阵A的奇异性，因此，需要进一步考虑当A为奇异矩阵时计算过程可能发生的情况。一是消元过程的某一步找不到非零的)(iiia，于是计算中断；二是虽然消元过程能进行到底，但0)(nnna，使回代求解过程无法进行下去，因此，在计算设计中必须考虑到上述两种可能发生的情形，此时应在算法设计中给出计算中断的信息。1Gauss列主元素消去法在Gauss消元过程中位于矩阵)2,1()(nkAk的主对角线),(kk位置上的元素)(kkka称为主元素，因为在计算解的分量时，都做除数，但应当避免用小的数做除数，即避免用较)(kika)1(nik在数量级上相对小的主元)(kkka做除数，以防止舍入误差的扩大，降低解得精度，下面的例子说明“小主元”对解得精度的影响。例1用Gauss消去法解线性方程组321643.5072.12623.4712.3132103218xxx用8位十进制尾数的浮点数计算解3643.5072.122623.4712.3113210),(),(8)1()1(bAbA88)1(11)1(212110101aam88)1(11)1(3131102102aam)102712.3()(8122122)2(22flamafla999102.0102.01000000000.0981212)2(2101.0)102()(flbmbflb类似的算出)2(3)2(33)2(32,,baa我们看到，由于小主元810做除数，使得行乘数3121,mm的数量级变得很大，这样在计算)2(22a时，即)712.3(fl与)102.0(9fl在计算机上做代数和时，)712.3(fl由于阶码升为9使尾数右移变成了机器零，这里揭示了数量级相对于分子数做除数，会使Gauss消去法在计算机有限字长的环境下造成较大舍入误差的原因。经第一步消元得9999998)2()2(102.0106.0104.00101.0103.0102.0013210),(bA经第二步消元得0000101.0103.0102.0013210),(9998)3()3(bA这一结果表明原方程组有无穷多个解，在计算机上，将会给出“无唯一解”的信息，停止计算，但实际上，由于0detA，所以原方程组有唯一解，是计算过程的舍入误差使解面目全非了，从前面的分析我们看到，导致这一错误的主要原因是用“小主元”做除数。为避免小主元做除数，在Gauss消去法中加入选主元过程，即在第k步消元时，首先在)(kA的第k列主对角线元以下元素中挑选出绝对值最大者，并通过交换的第与行对应元素，使位于对角线上，仍记为，并称之为第步消元的列主元素，然后再进行消元计算，第一步都按列选主元的消去法称为列主元素消去法，由于列主元素满足条件所以行乘数皆满足，因此，列主元素消去法能避免例1中所出现的问题，误差分析与数值试验表明，用列主元消去法解“良态”线性方程组，在绝大多数情况下，都可以得到令人满意的结果。例2用列主元素消去法解例1的线性方程组bAx解),(bA3643.5072.122623.4712.31132108132102623.4712.313643.5072.12883121105.0,5.0mm,经第一步消元，得11111)2()2(101.0103.0102.005.01018015.0103176.003643.5072.12),(bA62972292.032m，经第二步消元得68513854.01018655541.0005.01018015.0103176.003643.5072.12),(111)3()3(bA回代求解得Tx)36725739.0,05088607.0,49105820.0(方程组具有9位有效数字的解Tx)367257384.0,050886075.0,491058227.0(可见用Gauss列主元素消去法计算，得到了一个高精度的近似解。2-2带有行交换的矩阵分解设第k步消元时，先交换的第行与行后再消元，用矩阵运算表示为kikkIL),(),()1()1()()(kkkkbAbA其中1101111011kikI11111,2,1nkkkkkkmmmL),1(,)()(nkiaamkkkkikik当kik时kikikkIII,称为置换矩阵。于是带行交换的消元过程用矩阵运算表示为),(),()()(1122111121nniinnbAbAILILILn由此得UAAILILILniinnn)(1122111121（2.3）比如，4n则有UAILILILiii112233123可改写为UAIIIIILIIILILiiiiiiiii))()((123321233233123322333(2.4)令323233~iiILIL3212313223~iiiiIILIIL123123iiIIIPi容易证明，当,ijk时，kijkijLILI与kL的形状相同，也是一个初等下三角形矩阵，只是交换了kL的第k列对角元以下的第i行与第j行的元素，因此，1L和2L仍然是初等下三角形矩阵，从而321LLLL为单位下三角形矩阵，于是可以写成LPAU将对4n的处理方法推广到一般情形，令121121121121nkkkknkinikikkikikinLIIILIII（2.5）12211221nnininiiPIIII称P为排列矩阵，于是(2.3)式可以改写为1221nnLLLLPAU从而11221nnPALLLLULU其中11221nnLLLLL为单位下三角形矩阵，它与(1.4)中的L区别仅在于第k列的主对角元一下元素的排序可能不同，取决于消元过程所做的行交换，若设行交换的次数为m，则还有det(1)detmAU（2.7）式表明，带行交换的消元过程产生了矩阵PA的LU分解，相当于先对系数矩阵A施行消元过程中所需的行交换后再分解，只要A非奇异，带行交换的消元过程就可以进行到底定理2.1设A为非奇异矩阵，则必存在排列矩阵P以及单位下三角形矩阵L和非奇异上三角形矩阵U，使PALU(2.8)3直接三角分解法本节以矩阵分解定理和矩阵乘法规则为基础，研究矩阵分解的计算公式以及消去法的变形直接三角分解法当实现了系数矩阵的三角分解之后，解线性方程组原则上化为解三角形方程组(1.5)和(1.6)其解法在本章§1中已经基本解决了。定理2.1表明，A在满足一定条件下必有ALU或PALU成立，同时借助于Gauss消去法得到了L和U的表达式,这一节着重研究在ALU或PALU时，L，U的元素与A的元素之间的直接关系，并在此基础上给出求解矩阵方程组AXB的方法。3-1基本的三角分解法设矩阵ijnnAa的各阶顺序主子式均不为零，A有分解式111211112121222212221212111nnnnnnnnnnnnaaauuuaaaluuLUaaallu对比等式左边和右边乘积矩阵LU的第r行主对角元右边（含主对角元）的对应元素，得1,,;1,,13.2rijrkkikaluirnrn再对比等式两边第r列主对角元以下（不含主对角元）的对应元素，得11,,;1,,1