数值计算方法第三章

28第三章线性代数方程组的解法针对有惟一解的非齐次线性代数方程组求解问题，本章主要主要介绍Gauss(高斯)顺序消去法、主元素消去法、三角分解法等直接求解算法，以及Jacobi(雅可比)法、Gauss-Seidel(赛德尔)法和逐次超松弛等迭代求解算法．§3.1引言大量的科学与工程实际问题常常可以归结为求解含有多个未知量x1，x2，…，xn的线性代数方程组nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111(3.1)其对应的矩阵形式为bAx，其中n阶非奇异矩阵A以及n维列向量x和b分别定义如下1112111212222212nnnnnnnnaaaxbaaaxbxbaaaAxb线性代数方程组数值解法可以分为直接法和迭代法两类．所谓直接法，是指在没有舍入误差的假设下，经过有限步运算就能得到方程组精确解的一类方法；而迭代法则是从一个初始向量(0)x出发，按照一定的迭代格式产生一个向量序列()0{}kkx，当该序列收敛时，其极限就是线性方程组的解．§3.2Gauss消去法Gauss消去法的基本思想是通过逐步消元将线性方程组(3.1)化为系数矩阵为三角形矩阵的同解方程组，然后用回代算法解此三角形方程组，并得到原方程组的解．一、三角形方程组解法对于如下形式的下三角形方程组，29nnnnnnbxaxaxabxaxabxa221122221211111(3.2)若aii0(1,2,,)in，则(3.2)的解为11111122,11/()/(2,3,...,)kkkkkkkkkxbaxbaxaxaxakn(3.3)此方法称为前推算法．类似地，对于如下形式的上三角形方程组，1111221122222.........nnnnnnnnaxaxaxbaxaxbaxb(3.4)若aii0(1,2,,)in，则(3.4)的解为,11()/(1,2,,1)nnnnkkkkkknnkkbxaxbaxaxaknn(3.5)此方法称为回代算法．由于求解三角形方程组的过程很简单，因此只要把方程组化为等价的三角形方程组，求解过程就很容易完成．二、Gauss顺序消去法设线性方程组(3.1)的系数矩阵nnija)(A非奇异，为描述方便，记该线性方程组的增广矩阵为)1(1,)1()1(2)1(1)1(1,2)1(2)1(22)1(21)1(1,1)1(1)1(12)1(11nnnnnnnnnnaaaaaaaaaaaabA(3.6)其中(1)(1),1,(,1,2,...,)ijijiniaaabijn(3.7)Gauss顺序消去法实际上就是在求解过程中方程次序不变(即没有换行操作)的Gauss消去法，具体过程如下：30第1步设0)1(11a，用)1(11)1(1aai乘以(3.6)第1行后加到第i行(2,3,,)in，则第i行的第j个元素化为(1)(1)(1)(2)11(1)11:iijjijaaaaa(2,3,,1)jn(3.8)此时增广矩阵(3.6)相应地化为)2(1,)2()2(2)2(1,2)2(2)2(22)1(1,1)1(1)1(12)1(11nnnnnnnnnaaaaaaaaaa(3.9)第2步设0)2(22a，用)2(22)2(2aai乘以(3.9)第2行后加到第i行(n,,,i43)，则第i行的第j个元素化为(2)(2)(2)(3)22(1)22:iijjijaaaaa(3,4,,1)jn(3.10)增广矩阵(3.9)相应地化为)3(1,)3()3(3)3(1,3)3(3)3(33)2(1,2)2(2)2(23)2(22)1(1,1)1(1)1(13)1(12)1(11nnnnnnnnnnnaaaaaaaaaaaaaaa(3.11)类似地，当完成了第1~1n步消元后，(3.11)就化为)(1,)()2(1,2)2(2)2(22)1(1,1)1(1)1(12)1(11nnnnnnnnnnaaaaaaaaa(3.12)此时已把原方程组(3.1)等价转化成了上三角方程组(3.12)，可用回代算法求解该上三角方程组，回代公式为31(),1()()(),1()11(1,2,,1)nnnnnnnnkkkknkjjkjkkkaxaxaaxknna(3.13)由上面的消元过程知，Gauss顺序消去法的消元过程能进行下去的前提条件是()0(1,2,,1)kkkakn．回代过程则要求()0(1,2,,)kkkakn．例3.1用Gauss顺序消去法解线性方程组(计算过程保留4位小数)12120.000364.0001694xxxx解将第一个方程乘以6再除以0.0003然后加到第二个方程得1220.00036.00004.0001120009.000080006.0000xxx可得6667.02x，代入方程第一个方程得01x，而方程组的精确解为311x，322x，可见用Gauss顺序消去法得到的解完全失真．例3.1说明，当()0(1,2,,)kkkakn时，Gauss顺序消去法能够进行下去，但当0)(kkka或)(kkka相对于其他元素很小时，消去过程在浮点计算过程中产生的舍入误差将导致计算结果的误差急剧增大．在这种情况下，可以采用Gauss主元消去法．三、Gauss主元消去法对于上面的例3.1，如果我们先交换两个方程的顺序，得到等价方程组12126940.000664.0002xxxx计算过程中同样保留4位小数，用Gauss顺序消去法可解得6667.02x，3333.01x．32上面的求解过程实际上就是Gauss主元消去法．根据主元选取范围的不同，Gauss主元消去法又分为列主元消去法和全主元消去法．1Gauss列主元消去法Gauss列主元消去法的执行过程如下：首先，在增广矩阵(3.6)的第1列元素中选绝对值最大的元素)1(11ia，称之为第1列的主元，即有)1(11)1(1max1iniiaa如果11i，则交换增广矩阵中第1行和第1i行的元素，交换后的增广矩阵仍记为式(3.6)，但此时)1(11a已是第1列的主元．用主元)1(11a将其下边的1n个元素),,3,2()1(1niai消元为零的过程与Gauss顺序消去过程的第一步相同，从而得增广矩阵(3.9)．接着，在式(3.9)的第2列中除第1个元素外的其余1n个元素中选主元)2(22ia，即有)2(22)2(2max2iniiaa如果22i，则交换(3.9)中第2行和第2i行的元素，此时新的)2(22a已是第2列下边1n个元素的主元，同Gauss消去法过程中的第二步得到增广矩阵(3.11)．其余消元过程可以类似进行．当完成第1~1n步的选列主元及相应高斯消去过程后，则得增广矩阵(3.12)，最后利用回代公式(3.13)求得原方程组(3.1)的解．列主元消去法除了每步需要按列选主元并作相应的行交换外，其消去过程与Gauss顺序消去法的消去过程相同．归纳起来有如下算法．算法3.1（列主元消去算法）输入方程组阶数n、二维数组A存放的增广矩阵(|)Ab；输出方程组的解(1,2,,)ixin；Step1选主元的消去过程对121n,,,k33(1)找行号n,,k,kik1，使)()(maxkiknikkkiaak；(2)若kik则交换增广矩阵A的第k行和第ki行元素；(3)消元：对n,,k,ki21，计算)()(kkkkikikaal，对121n,,k,kj，计算)()()1(kkjikkijkijalaa；Step2回代过程(1)计算)()(1,nnnnnnnaax；(2)对121,,n,nk，计算nkjjkkjknkkkkkxaaax1)()(1,)(1；Step3输出方程组的解(1,2,,.ixin)不同于Gauss顺序消去法，Gauss列主元消去法只要线性方程组的系数矩阵A非奇异，即可求出其惟一解．2Gauss全主元消去法全主元消去法选主元的范围更大．对增广矩阵(3.6)来说，第1步是在整个系数矩阵中选主元，即选绝对值最大的元素，经过行列交换将其放在11a元素的位置，然后实施第一步消元过程．第二步是在(3.9)中的1n阶矩阵)2()2(2)2(2)2(22nnnnaaaa中选主元，其余各步选主元的范围依次类推，而每步选主元后的消去过程同列主元法的消去过程．全主元消去法每步所选主元的绝对值，通常比列主元消去法所选主元的绝对值要大，因而全主元消去法的求解结果更加可靠．但全主元消去法每步选主元要花费更多的机时，并且由于对增广矩阵可能进行列交换，使得未知量nx,,x,x21的次序发生了变化，在编程时应当注意．§3.3矩阵三角分解法34前面3.2节中介绍过，对于下三角形方程组和上三角形方程组很容易通过前推算法和回代算法求解．如果能将线性方程组bAx的系数矩阵A分解成LUA(其中L是下三角阵，U是上三角阵)，则通过求解下三角方程组Lyb得到列向量y，再通过求解上三角方程组Uxy即可得到原方程组的解x．一、矩阵三角分解的基本概念Gauss顺序消去过程实际上就是对增广矩阵(|)Ab进行若干次初等行变换，使之化为(|)Ug的形式，其中U为上三角矩阵．由线性代数理论知，对一个矩阵进行一次初等行变换，相当于给这个矩阵左乘一个相应的初等矩阵．定义如下的一组单位下三角矩阵（对角线为１的下三角矩阵）213111101001nlll1L，32210101001nll2L，……，1,,11101kknkllkL，……，,11111nnln-1L其中)()(kkkkikikaal(1,2,...,1;1,2,,)knikkn则对应的消元过程依次为(1)(1)(2)(2)1LAbAb……(k)(k)(k+1)(k+1)kLAbAb……(1)(1):[]nng(n)(n)n-1LAbAbU35归纳消去过程有1221nnLLLLAbUg(3.14)由(3.14)得1221nnLLLLAU(3.15)记112()nn21LLLLL(L仍是单位下三角阵)，则有,ALULgb(3.16)定义3.1若(2)nn阶矩阵A可以分解为一个下三角阵L和一个上三角阵U的乘积，即LUA，则称这种分解为矩阵A的一个三角分解或LU分解；若L是单位下三角阵，U是上三角阵，则称为Doolittle(杜里特尔)分解；如果L是下三角阵，U是单位上三角阵，则称为Crout(克劳特)分解．定理3.1[26]如果n阶(2n)矩阵A的前1n阶顺序主子式均不为零，则A有惟一Doolittle分解和惟一Crout分解．下面我们以Doolittle分解为例说明矩阵三角分解的具体过程．二、Doolittle分解设矩阵A的所有顺序主子式不为零，并且有Doolittle分解LUA，则求解线性方程组bAx等价于求解下三角方程组Lyb及上三角方程组Uxy,即