第5章-解线性方程组的直接法 计算方法

第5章解线性方程组的直接法

在自然科学和工程技术中,很多问题归结为解线性方程组.有的问题的数学模型中虽不直接表现为含线性方程组,但它的数值解法中将问题“离散化”或“线性化”为线性方程组.因此线性方程组的求解是数值分析课程中最基本的内容之一.

线性方程组:结束常记为矩阵形式

Ax=b(5.2)

1

此时A是一个n×n方阵,x和b是n维列向量.

根据线性代数知识若|A|≠0,(5.2)的解存在且唯一.

关于线性方程组的解法一般分为两大类,一类是直接法,即经过有限次的算术运算,可以求得(5.1)的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法对高阶方程组计算量太大,不是一种实用的算法.实用的直接法中具有代表性的算法是高斯消元法,其它算法都是它的变形和应用.

另一类是迭代法,它将(5.1)变形为某种迭代公式,给出初始解x0,用迭代公式得到近似解的序列｛xk｝,k=0,1,2,,在一定的条件下xk→x*(精确解).迭代法显然有一个收敛条件和收敛速度问题.

这两种解法都有广泛的应用,我们将分别讨论,本章介绍直接法.结束

2

3§5.1消元法5.1.1三角形方程组的解

下面三种形式的线性方程组较容易求解.

若

则

运算量：1.对角形方程组

4

若

则

运算量：求xi需要i-1次乘法和1次除法，i=1,…,n2.下三角方程组

代入第2方程，可得

共1+2+…+n=n(n+1)/2=O(n2)

53.上三角方程组

若

则

运算量：求xi需要n-i次乘法和1次除法，i=n,n-1,…,1

代入倒数第2方程，可得

共1+2+…+n=n(n+1)/2=O(n2)

高斯消元法是一种古老的方法.我们在中学学过消元法,高斯消元法就是它的标准化的、适合在计算机上自动计算的一种方法.1高斯消元法的基本思想例1

解方程组

(5.3)(5.4)(5.5)第一步,将(5.3)乘-2加到(5.4)；(5.3)乘-1加到(5.5),

得到

(5.3)(5.6)(5.7)

6结束5.1.2Gauss消元法与列主元消元法第二步,将(5.6)乘-2/3加到(5.7),得到

(5.3)(5.6)(5.8)回代:解(5.8)得x3,将x3代入(5.6)得x2,将x2,x3代入(5.3)得x1,得到解x*=(2,1,-1)T

容易看出第一步和第二步相当于增广矩阵［Ａ:b］在作行变换,用ri表示增广阵［A:b］的第i行:结束

7

由此看出上述过程是逐次消去未知数的系数,将Ax=b化为等价的三角形方程组,然后回代解之,这就是高斯消元法.2高斯消元法公式

记Ax=b为A(1)x=b(1),A(1)和b(1)的元素记为和,i,j=1,2,,n.第一次消元,目的是消掉第二个方程到第n个方程中的x1项,得到A(2)x=b(2),这个过程须假定≠0.结束

8

在［A(1):b(1)］中,红方框中的元素是要化为0的部分；［A(2):b(2)］中,红方框中的元素全部已发生变化,故上标由(1)改(2),计算公式为:结束

9(i=2,3,,n)(i,j=2,3,,n)(i=2,3,,n)(i=2,3,,n)第k次消元(1≤k≤n-1)

设第k-1次消元已完成,且≠0,此时增广矩阵如下:结束

10

本次消元的目的是对框内部分作类似第一次消元的处理,消掉第k+1个方程到第n个方程中的xk项,即把到化为零.计算公式如下:(i=k+1,,n)(i,j=k+1,,n)(i=k+1,,n)(i=k+1,,n)

只要≠0,(k=1,2,,n-1)消元过程就可以进行下去.当k=n-1时,消元过程完成,得:结束

11

它的方阵部分A(n)是一个上三角形矩阵,它对应的方程组是一个上三角形方程组,只要≠0,就可以回代求解,公式为(i=n-1,n-2,,1)综合以上讨论,高斯消元法解线性方程的公式为:1消元①令

(i,j=1,2,…,n)结束

12(i=k+1,k+2,,n)

(i,j=k+1,k+2,,n)

(i=k+1,k+2,,n)(5.9)2回代,若≠0(i=n-1,n-2,,1)(5.10)结束②对k=1到n-1,若≠0,进行:

(i=k+1,k+2,,n)

133高斯消元法的条件以上过程中,消元过程要求≠0(i=1,2,,n-1),回代过程则进一步要求≠0,但就方程组Ax=b讲,是否等于0是无法事先看出的.

注意A的顺序主子式Di

(i=1,2,,n)在消元过程中不变.这是因为消元所作的变换是“将某行的若干倍加到另一行”上,据线性代数知识,此类变换不改变行列式的值.若高斯消元过程已进行了k-1步(此时当然应有≠0,i≤k-1),这时计算A(k)的顺序主子式:结束

14

有递推公式(i=2,3,,k)显然,,可知,消元过程能进行到底的充要条件是Di≠0,(i=1,2,,n-1),若要回代过程也能完成,还应加上Dn=｜A｜≠0,综合上述有:定理5.1

高斯消元法消元过程能进行到底的充要条件是系数阵A的1到n-1阶顺序主子式不为零；Ax=b能用高斯消元法解的充要条件是A的各阶顺序主子式不为零.4高斯消元法的计算量估计消元过程的工作量,参看公式(5.9),k是消元次数,k=1,2,,n-1,第k步消元时,计算lik(i=k+1,,n)需要n-k次除法；计算(i,j=k+1,,n)需要(n-k)2次乘法及(n-k)2次加减法；计算需要n-k次乘法及n-k次减法,合计:结束

15

乘除法次数

加减法的次数

回代过程的工作量,参见公式(5.10),求xk需n-k次加减法,n-k次乘法和1次除法,合计为乘除法次数

加减法次数

总的运算次数为乘除法(当n较大时)结束

16

加减法(当n较大时)

一般讲乘除法的运算比加减法占用机时多得多,往往只统计乘除法次数而称高斯消元法的运算量为次.在上节的算法中,消元时可能出现=0的情况,高斯消元法将无法继续；即使≠0,但<<1,此时用它作除数,也会导致其它元素数量级严重增加,带来舍入误差的扩散,使解严重失真.例2

线性方程组准确解是(1,1)T.结束

175列主元的高斯消元法现设我们的计算机为四位浮点数,方程组输入计算机后成为(5.11)用高斯消元法:l12=0.2×106,r2=r2-l12×r1回代:x2=0.1000×10=1,x1=0,解严重失真.

若将r1和r2交换.消元,l12=0.5×10-5,r2=r2-l12×r1结束

18回代,x1=0.1000×10,x2=0.1000×10,得到准确解.

从上例中可以看出,对方程组作简单的行交换有时会显著改善解的精度.在实际使用高斯消元法时,常结合使用“选主元”技术以避免零主元或“小主元”出现,从而保证高斯消元法能进行或保证解的数值稳定性.列主元消元法设已用列主元消元法完成Ax=b的第k-1次消元(1≤k≤n-1),此时方程组Ax=b→A(k)x=b(k),有如下形式:(5.12)结束

19全为零,此时易证｜A｜=｜A(k)｜=0,即方程组Ax=b无确定解,应给出信息退出计算.

进行第k次消元前,先进行①、②两个步骤:①在方框内的一列内选出绝对值最大者,即

,确定ik.若=0,则必有方框内的元素②若≠0,则交换ik行和k行元素即(kjn)然后用高斯消元法进行消元.

这样从k=1做到k=n-1,就完成了消元过程,只要｜A｜≠0,列主元高斯消元法必可进行下去.结束

20很多问题都需要计算方阵的逆阵.线性代数中有多种求逆方法,本节讨论求A-1的数值方法.A非奇异,求A-1

设X=A-1,AX=I,将X和I列分块,得n个线性方程组:

Axi=ei,i=1,2,,n.解这n个线性方程组就求得A-1,原则上,所有解的方程组的方法都能求A-1.实际计算中使用较多的是高斯-若当消元法.高斯-若当消元法是高斯消元法的一种变形.高斯消元法是消去对角元下方的元素.若同时消去对角元上方和下方的元素,而且将对角元化为1,就是高斯-若当消元法.

设高斯-若当消元法已完成k-1步,于是Ax=b化为等价方程组A(k)x=b(k),增广阵结束

215.1.3高斯-若当(Gauss-Jordan)消元法在第k步计算时,考虑将第k行上下的第k列元素都化为零,且akk化为1.对k=1,2,,n1按列选主元,确定ik使(5.13)2换行,交换增广阵第k行和第ik行(j=k,k+1,,n),结束

223计算乘数mik=-aik／akk(i=1,2,,n,i≠k)

mkk=1／akk.(5.14)4消元

aij=aij+mikakj(i=1,2,,n,且i≠k,j=k+1,,n)

bi=bi+mikbk

(i=1,2,,n,且i≠k)(5.15)5主行计算

akj=akj×mkk

(j=k,k+1,,n)

bk=bk×mkk

(5.16)当k=n时,

显然xi=bi,i=1,2,,n就是Ax=b的解.结束

23

高斯-若当消元法的消元过程比高斯消元法略复杂,但省去了回代过程,它的计算量约为n3/2,大于高斯消元法.也称为无回代的高斯消元法.结束

24§5.2直接分解法（LU分解）

25若记，则

26若记，则

27于是同理：结束

285.2.1Dolittle分解由上式得

29而令则：

30令：则：称为单位下三角阵称为上三角阵称为矩阵的LU分解或Dolittle分解定理

矩阵An×n,只要A的各阶顺序主子式非零,则A可以分解为一个单位下三角阵L和一个上三角阵U的乘积,即A=LU,且这种分解是唯一的.若记第k步消元则有:结束

31当A进行LU分解后,Ax=b就容易解了.即Ax=b等价于:结束这样,Ax=b分解为解两个三角形方程组,三角形方程组是极易求解的.(5.18)即是高斯消元法的矩阵表述.

32

33L,U的求法：1.U的第一行元素得到：

342.L的第一列元素得到：……，若已计算出U的前k-1行和L的前k-1列元素，则

353.U的第k行元素：ukk,uk,k+1,…,ukn当列标j>=行标k时，有得到：(5.19)

364.L的第k列元素：lk+1,k,lk+2,k,…,ln,k当行标i>=列标k时，有得到：(5.20)1.解Ly=b2.解Ux=y结束

37杜利特算法实际上就是高斯消元法的另一种形式.它的计算量与高斯消元法一样.但它不是逐次对A进行变换,而是一次性地算出L和U的元素.L和U的元素算出后,不必另辟存贮单元存放,可直接存放在A的对应元素的位置,节省存贮单元,因此也称为紧凑格式法.AX=b的解：杜利特算法可用文字描述如下:将A的元素划分为形如“┏”的n框,如A的第一行和第一列元素为第一框,紧靠第一框内部的第二行和第二列元素为第二框,依此类推,ann一个元素为第n框.为便于描述我们称原A矩阵中元素为a元素,计算后每框中的每行元素为u元素(包含对角元素),每列元素为l元素(不包含对角元素).这样我们可以在A上直接修改计算,作如下的操作:(1)第1框:u元素等于对应的a元素;l元素等于对应的a元素除以本框对角元素;(2)第2到n框:u元素等于对应的a元素减去已经算过的各框中同行同列元素的乘积;l元素除进行以上操作外,还要除以本框对角元素结束

38例用Dolittle分解求解方程组结束

39

得到的矩阵由L和U的元素拼合而成,它的上三角部分(包含对角元素)即是LU分解后的U矩阵;它的下三角部分(不包含对角元素