第3章 线性方程组的数值解法

第3章线性方程组的数值解法§1高斯消去法§2高斯―约当消去法§3解实三对角线性方程组的追赶法§4矩阵的三角分解§5行列式和逆矩阵的计算§6迭代法

§7迭代法的收敛性§1高斯消去法1.1顺序消去法如果线性方程组(3―1)的系数矩阵A具备某特殊形式,例如其为上三角矩阵且aii≠0,i=1,2,…,n这时方程组(3―1)实际为(3―4)由方程组(3―4)的最后一个方程直接可得将其代入倒数第二个方程可求得如此再解出xn-2,…,x2,x1,一般有(3―5)其中规定以上讨论告诉我们,对具有上三角形系数矩阵的方程组(3―4)求解极为方便。当然,若方程组(3―1)的系数矩阵为下三角形,则求解也很方便。于是对于一般形式的方程组(3―1),我们总设法把它化为系数矩阵呈上(或下)三角形的方程组来求解。为了达到目的,可利用消去法进行。现举例如下:解方程组①②③(3―6)作②-①消去②中的x1,作③-①×4消去③中的x1,则方程组(3―6)化为①②③(3―6′)对方程组(3―6′)作③-②×,得到三角形方程组①②③(3―6")从方程组(3―6“)的方程③解出x3,将所得的结果代入方程②求出x2,再把x3、x2同时代入方程①解出x1。这样可求出方程组的解为

上述求解方程组的方法就是高斯(Gauss)消去法。从式(3―6)到(3―6")的过程称为消元过程而由(3―6")求出x3、x2、x1的过程称为回代过程。因此用高斯消去法求解性方程组要经过消元和回代两个过程。在计算机上实现时,我们常把方程组右端的常数项排于系数矩阵的第n+1列,这样顺序高斯消去法的计算步骤为:1.消元过程对于k=1,2,…,n-1,若按顺序有某一ark≠0,r≥k,则交换k与r行,然后计算

(3―11)2.回代过程对于k=n,n-1,…,2,1,计算(3―12)1.2主元素消去法前述顺序消去法是按序通过用a11,a(1)22,…,a(n-2)n-1(a(k-1)kk≠0)作为除数来达到消元目的的。在

实际计算时,由于舍入误差的影响,计算结果会改变很大,甚至于完全失真。例如用高斯消去法求解下列方程组(用四位有效数字计算):①②②-①×105得(1.000-1.000×105)x2=1.000-6.000×104化简可得

x2=0.6000回代求得

x1=105(0.6-0.6000)=0而方程组的解应为

x1=04000x2=06000显然用上述方法求出的解x1与方程组的实际解相差很大。若改变两个方程的顺序,即

x1+x2=1

①10-5x1+x2=0.6②②-①×10-5得(1000-1000×10-5)x2=0.6-1.000×10-5

化简得

x2=0.6000

回代求得

x1=(1-0.6000)=0.4000高斯主元素消去法是顺序消去法的一种改进。它的基本思想是在逐次消元时总是选绝对值最大的元素(称之为主元)做除数,按顺序消去法的步骤消元。这里主要介绍求解线性方程组最常用的列主元素消去法和全主元素消去法。1.列主元素消去法所谓列主元素消去法就是在每一步消元过程中取系数子矩阵的第一列元素中绝对值最大者作主元。对线性方程组(3―1)进行n-1次消元后,可得到上三角形方程组必须指出的是方程组(3―13)中的系数aij(i≤j)和右端的bi已经改变了,并非与原来相同。这样就可对方程组(3―13)回代求解。例1用列主元素消去法解方程组(3―13)取四位有效数字计算。解②中-18为主元,交换②和①得①②③①②③②+①×12/18,③+①×1/18得①②③第二列消元时,主元为1167,交换方程②和③得①②③③+②×1/1167得①②③回代求得

x1=1000,x2=2000,x3=3001方程组的实际解

x1=1,x2=2,x3=3在计算机上求解时,前面已经讲过,常把右端常数项b作为系数矩阵A的第n+1列,从而得到增广矩阵(A,b),仍记为A,于是,列主元素消去法的计算过程为:(1)消元过程。对于k=1,2,…,n-1进行下述运算:①选主元,确定r,使得若ark=0,说明系数矩阵为奇异,则停止计算否则进行下一步。②交换A的r、k两行③对i=k+1,k+2,…,n,j=k+1,k+2,…,n+1计算(3―14)

aij-aik·akj/akkaij

(3―15)(2)回代过程。对于k=n,n-1,…,1计算(3―16)2.全主元素消去法所谓全主元素消去法,就是每步消元时选取系数子矩阵中绝对值最大的元素作主元。经过n-1次消元后,方程组(3―1)可化为上三角形方程组(3―17)这里方程组(3―17)中的系数aij(i≤j)及bi一般改变了。特别是未知数的排列顺序,由于全主元素的消元过程中,其系数矩阵可能进行了列对调,那么未知数也相应地作了对调。例如,若矩阵第i列与第j列对调,则未知数xi与xj也相应地对调了,xi的结果实质上为xj的结果。图3.1图3.1例2用全主元素消去法求解方程组①②③解这里主元为-18,交换方程①与方程②得①②③再全选主元,主元为2333,交换x2和x3所在的两列,同时改变两未知数的排列号得①②③③-②×0944/2333得①②③已经化为三角方程组,回代求解

x1=1000,x2=3000,x3=2000

这里未知数x2与x3已对调,所以应恢复解的顺序,方程组的实际精确解为

x1=1000,x2=2000,x3=3000

在计算机上求解时,我们仍把增广矩阵(A,b)记为A,具体计算过程为(1)消元过程。对k=1,2,…,n-1进行下列运算:图3.2图3.2①选主元,确定r,t使得若art=0,则系数矩阵为奇异的,停止计算否则进行下一步。②交换A中的r、t两行及t、k两列,并记下交换的足码t、k。③对i=k+1,k+2,…,n,j=k+1,k+2,…,n+1计算(2)回代过程。对k=n,n-1,…,1计算(3―18)(3―19)(3―20)§2高斯―约当消去法前面所述的消去法均要进行两个过程,即消元过程和回代过程。但对消元过程稍加改变可以把方程组(3―1)化为对角形(3―21)此时求解就不要回代了。这种无回代过程的主元素消去法称为高斯―约当(Jordan)消去法。特别是方程组(3―21)还可化为(3―22)显然等号右端即为方程组的解。对于n阶线性方程组(3―1),其增广矩阵为首先把主元素(按列选主元或全选主元)调换到主对角线上,并化为1,再将主元素所在列的其它元素消为0,则第一次消元后增广矩阵化为显见经n步消元后即得方程组的解。其具体计算步骤为:对k=1,2,…,n①选主元(列选或全选主元)②对j=k,k+1,…,n+1计算(3―23)③对i=1,2,…,n(i≠k),j=k+1,k+2,…,n+1计算

aij-aikakjaij

(3―24)④计算

xk=a

kn+1

还得指出,若用到全选主元,最后应注意恢复解的顺序。§3解实三对角线性方程组的追赶法实三对角线性方程组是一类特殊的方程组。我们将利用所谓“追赶法”解决。设所给的实三对角线性方程组为其中｜a1｜＞｜c1｜＞0｜ak｜≥｜bk｜+｜ck｜,bkck≠0,k=2,3,…,n-1｜an｜＞｜bn｜＞0

则三对角矩阵A的各阶主子式都不等于0,也即A为非奇异。用数学归纳法证明:令A的k阶主子式为Δk,Δ2=a1a2-b2c1≠0

其中因为所以方程式(3―26)右端的k-1阶行列式满足假设的条件,由此可得Δk≠0,即矩阵A的各阶主子式都不等于0。由于实三对角线性方程组的系数矩阵A为非奇异,则有唯一的一组解。下面来求解。从方程组(3―25)的第一个方程解得令(3―27)(3―28)将式(3―28)代入(3―25)式的第二个方程解出x2(3―29)令则再将(3―29)式代入(3―25)式的第三个方程解出x3,一般有

xk=uk-vkx

k+1,k=1,2,…,n-1(3―30)

其中k=2,3,…,n-1(3―31)当k=n时,因为

bnx

n-1+anxn=rn又

xn-1=un-1-vn-1xn于是(3―32)例3解线性方程组解依公式(3―27)、(3―31)求得v3=0(一般v3不必算)再由式(3―32)、(3―30)求得方程组的解图3.3图3.3

§4矩阵的三角分解

对于线性方程组(3―1),若其系数矩阵A能够分解成两个三角形矩阵相乘,譬如,

A=LUL为下三角矩阵,U为上三角矩阵,则求解方程组(3―1)就化为求解

LUx=b令

Ux=y(3―34)(3―35)则Ly=b(3―35)、(3―36)都为三角形方程组,先求出(3―36)中的y,然后代入(3―35)就可以求出原方程组(3―1)的解x。(3―36)4.1消去法与矩阵的初等变换从矩阵的初等变换来看,前面所述的主元素消去法(假定主元素经过行列调换以后均在对角线上),就是通过一系列的初等变换将方程组(3―1)的系数矩阵A化为三角矩阵。其每一步的消元过程相当于对增广矩阵(A,b)作一次初等变换。令在消去第一列中的ai1(i=2,3,…,n)时,左乘的矩阵为这样就把A化为一个上三角矩阵U，即在需要进行行交换时，还得左乘某些排列矩阵Pij

A=LU(3―37)

称式(3―37)为矩阵A的LU分解或三角分解。由上述讨论,只要系数矩阵A非奇异,经过适当的行交换后总可以将A分解为一个下三角矩阵L和一个上三角矩阵U之积。下面首先讨论矩阵三角分解的唯一性。4.2矩阵三角分解的唯一性设非奇异矩阵A存在一种三角分解LU,一般这种分解未必唯一,这是因为任意一个同阶的非奇异对角矩阵D总有(LD)(D-1U)=L1U1

它也是矩阵A的一种三角分解。为了确保三角分解的唯一性,需对分解加以规格化。为使LU分解规格化,可以把矩阵U换成DU,其中L为单位下三角矩阵,U为单位上三角矩阵,D为对角矩阵,即我们称

A=LDU(3―38)

为矩阵A的一种LDU分解。关于矩阵A的LDU分解有如下定理:定理1对于n×n阶矩阵A存在唯一的LDU分解的必要充分条件是A的各阶主子矩阵A1,A2,…,An均非奇异。证先证必要性。设A有唯一的LDU分解,即

A=LDU

由于

Ak=LkDkUk

其中Ak、Lk、Dk、Uk分别表示A、L、D、U的k阶主子矩阵。因为L、U是单位三角矩阵,D为非奇异对角矩阵,所以Lk、Uk、Dk均非奇异,从而知Ak必为非奇异矩阵。再证充分性。用数学归纳法证明:当k=1时,A1显然有这种三角分解且是唯一的,即

A1=L1D1U1=l11d1u11

其中

l11=1,d1=a11,u11=1设对k=n-1时结论成立,将A按下面方式分块:这里

s=(a1n,a2n,…,an-1n)T

rT=(an1,an2,…,a

nn-1)a=ann再将A、D、U仿照A分块:这里设若通过A的元素能唯一地定出L中的IT,U中的u及D中的d,就可确认A唯一地有LDU分解了。因为其中r、s、a均为已知,于是有

Ln-1Dn-1u=s(Dn-1Un-1)TI=rITDn-1u+d=a由于A的各阶主子矩阵非奇异,因而Ln-1、Dn-1、Un-1

均为非奇异,于是可以唯一地求出

所以,A存在唯一的LDU分解。证毕。由定理1不难得到定理2(证明留给读者)。定理2n×n阶矩阵A非奇异的必要充分条件为A(或A经行、列变换后)存在LDU分解。推论1奇异矩阵不能进行LDU分解。推论2若方阵A有奇异主子矩阵,则A不能直接进行LDU分解。证用反证法。若A能直接进行LDU分解,即A=LDU。设A的k(k≤n)阶主子矩阵Ak奇异,则按分块矩阵的乘法有

Ak=LkDkUk

而Lk为单位下三角矩阵,Dk为非奇异对角矩阵,Uk为单位上三角矩阵,显然与推论1矛盾,所以结论成立。可验证A为非奇异矩阵,但其左上角二阶主子矩阵奇异,故由推论2不能直接进行LDU分解。因A为非奇异,由定理2,经适当行变换后存在LDU分解。若经行变换例设则LDU分解为A1=LDU若经列变换则A2的LDU分解为4.3LU分解方法我们已经知道,若矩阵A的各阶主子矩阵非奇异,则存在唯一的LDU分解,即

A=LDU

常常称(LD)U=L1U,L(DU)=LU1(3―39)

分别为矩阵A的克劳特(Crout)分解和杜利特尔(Doolittle)分解。仍记为

A=LU只不过若是克劳特分解,则L为下三角矩阵,U为单位上三角矩阵若是杜利特尔分解,此时L为单位下三角矩阵,U为上三角矩阵。显然这两种分解均是唯一的。现分别介绍求解线性方程组的两种LU分解方法。1.克劳特分解方法设A为n×n阶非奇异矩阵,且各阶主子矩阵为非奇异,则矩阵A的克劳特分解为

A=LU(3―40)

其中

L、U中元素的确定,可不用消去法,而直接由式(3―40)来求得。按矩阵乘法规则有因为u11=1,所以

ai1=li1u11=li1,i=1,2,…,n

于是可求得L的第一列元素为

li1=ai1,i=1,2,…,n(3―42)

又因为

a1j=l11u1j,j=2,3,…,n

所以

u1j=a1j/l11,j=2,3,…,n(3―43)(3―41)于是U的第一行元素也就由式(3―43)计算出来了。假定L的前k-1列以及U的前k-1行元素都已算出,由于ukk=1,于是从式(3―41)有(3―44)(3―45)这样,L、U中的元素都已求出。计算L的各列与U的各行的次序如图3。4所示。图3.4对方程组(3―1)的系数矩阵A作出LU分解后,方程组便化为

LUx=b

则求解上列方程组就化为依次解方程组

Ly=b(3―46)

Ux=y(3―47)

由于L为下三角矩阵,U为单位上三角矩阵,故方程组(3―46)、(3―47)的求解极为方便。先由(3―46)式解出y,然后代入(3―47)式求得x,亦即求得了方程组(3―1)的解。而求解(3―46)、(3―47)式的计算公式分别为用克劳特分解求解线性方程组(3―1)的计算过程为:①LU分解过程:对于k=1,2,…,n依次计算(3―48)(3―49)从上述计算过程还可看出,求解yk的公式完全可以插在计算L、U元素的两个公式之间。其计算框图见图3.5。图3.5图3.5例4用克劳特分解方法求解下列方程组解令利用矩阵乘法可得到这样原方程组就化为依次求下列两个三角形方程组代入第二个方程组可求得原方程组的解为2.杜利特尔分解方法杜利特尔分解方法是令

A=LU

这里关于L、U的元素可按以下公式计算:对于k=1,2,…,n进行(3―50)(3―51)在L与U的元素算出后就可以按

Ly=b(3―52)

Ux=y(3―53)解方程组LUx=b,从而求得方程组(3―1)的解。例5用杜利特尔分解方法求解下列方程组的解:解利用式(3―50)、(3―51)可求得再由式(3―52)、(3―53)求得原方程组的解为从上面的讨论我们可以看到,用LU分解方法解线性方程组(3―1),就是将其化为依次求解下面两个三角形方程组:

Ly=b

Ux=y显然求解Ly=b的过程即为消元过程,它只不过是对

LUx=b

两端左乘L-1而得,即

Ux=L-1b=y

而求解Ux=y的过程即为回代过程。因此三角分解方法实质上还是高斯消去法。而克劳特分解方法中的元素lii(i=1,2,…,n)和杜利特尔分解方法中的元素uii(i=1,2,…,n)实际上就是高斯消去法中的各次主元素。4.4乔累斯基(Cholesky)分解方法设A为n×n阶正定矩阵,记为A＞0。根据前节的讨论,其有唯一的分解式

A=LDU

其中L为单位下三角矩阵,D为非奇异对角矩阵,U为单位上三角矩阵。由于A是对称的,于是有

U=LT

从而A的LDU分解式可写成

A=LDLT(3―54)

又因A为正定矩阵,故D的主对角元素均为正数。记则A可以表示成

为简便起见,我们去掉L1的下标,仍把A的这种分解记作

A=LLT(3―55)

其中L是一个下三角矩阵。分解式(3―55)便是通常所说的正定矩阵A的乔累斯基分解。下面介绍两种具体的乔累斯基分解方法。1.LLT分解方法(平方根法)对于分解式由矩阵乘法规则从而得(3―56)(3―57)当i=j时也即于是这样由式(3―56)、(3―57)确定L的元素后,就可按

Ly=b(3―58)

LTx=y

(3―59)解方程组

LLTx=b(3―60)其解x即为方程组

Ax=b的解。由式(3―58)、(3―59)依次计算y和x可按下列递推公式进行:(3―61)(3―62)例6用LLT分解方法解线性方程组取四位小数计算。解用公式(3―57)、(3―56)依次计算由式(3―61)求得

y1≈1.3416,y2≈-0.4474,y3≈1.7320再由式(3―62)求得原方程组的解为

x1≈0.9993,x2≈0.9989,x3≈0.9999实际上,原方程组的准确解为

x1=1,x2=1,x3=1LLT分解方法的计算框图见图36。2.LDLT分解方法

LLT分解方法的计算中含有开方运算,而LDLT分解方法就避免了这一点。设A＞0,则其有唯一的分解式

A=LDLT

其中图3.6图3.6类似于前面的讨论,由矩阵乘法规则可得(3―63)也即因为ljj=1,从而有故有当i=j时(3―64)根据式(3―63)、(3―64)求出矩阵L、D的元素后,就可按照

Ly=b(3―65)

LTx=D-1y(3―66)解方程组LDLTx=b

其解x就是线性方程组(3―1)的解。具体计算x和y的递推公式为(3―67)(3―68)图3.7图3.7§5行列式和逆矩阵的计算5.1行列式的计算行列式的数值计算,只要将行列式化为三角形式(上三角或下三角),求解就很方便了。因此,前面所述的解线性方程组的方法都可以用来化行列式detA为三角形式。1.高斯消去法设detA为n阶行列式,即按矩阵A经若干次初等变换以后,其主元素分别为r11,r22,…,rnn,于是2.LU分解法①克劳特分解法:detA=detL·detU=detL=l11l22…lnn②杜利特尔分解法:detA=detL·detU=detU=u11u22…unn③乔累斯基分解法:这里A为正定矩阵。detA=detL·detLT=l211l222…l2nn或者detA=detL·detD·detLT=detD=d1d2…dn5.２逆矩阵的计算设A为非奇异矩阵,则A可逆,即有AA-1=I令A-1=X=(x1,x2,…,xn)于是有AA-1=A(x1,x2,…,xn)=(e1,e2,…,en)从而求逆矩阵A-1的问题就化为求解下列几个有相同系数矩阵A的线性方程组问题

Axi=ei,i=1,2,…,n(3―69)

因此,前面介绍的求解线性方程组的方法均可用来计算逆矩阵。例7设求A-1。解将方程组写成同一个增广矩阵,即先对第一列消元,将a11化为1,a21、a31化为0,有再将第二列a32化为0,得将a33化为1,并将a23、a13化为0,得最后将a13化为0,得所以§6迭代法6.1简单迭代法设所给的线性方程组为

Mx=g(3―70)

其中且系数矩阵M为非奇异,g≠0。将方程组(3―70)改写成等价形式

x=Ax+f(3―71)

这种改写的方法很多,例如将M分解为两个矩阵之差

M=B-C

其中矩阵B可逆,于是方程组(3―70)成为

Bx=Cx+gx=B-1Cx+B-1g

令

A=B-1C,f=B-1g即得(3―71)式。当然选取的B应该便于求逆,如B为单位矩阵或对角矩阵等。对任意给定的初始向量x(0),根据(3―71)构造迭代公式

x(k+1)=Ax(k)+f,k=0,1,2,…(3―72)

这里A称为迭代矩阵。用分量形式可写成(3―73)算出迭代公式(3―72)的解的各次近似值x(1),x(2),…,x(k),…。这种迭代求解的方法称为简单迭代法。式(3―72)称为简单迭代公式。特别当M的对角元素均不为零且按绝对值来说较大时,常取则

A=D-1C,f=D-1g例8用简单迭代法解下列方程组解将方程组写成等价形式取初始值x(0)=0,按迭代公式表3―1图3.86.2赛德尔(Seidel)迭代法从简单迭代法看到,已知x(k)计算x(k+1)时,需要保留x(k)和x(k+1)两个分量。逐次用前面算出的新分量来计算下一个分量,这就是赛德尔迭代法的基本思想。设方程组(3―70)的等价形式为

x=Ax+f

将矩阵A分解为

A=B+C

其中于是

x=Bx+Cx+f

(3―74)x(k+1)=Bx(k+1)+Cx(k)+f,k=0,1,2,…(3―75)例9用赛德尔迭代法解方程组解将原方程组写成等价形式并按(3―75)构造赛德尔迭代公式表3―2由式(3―75)可得(I-B)x(k+1)=Cx(k)+f

因为矩阵I-B是非奇异矩阵,则I-B可逆,所以迭代公式(3―75)可写成

x(k+1)=(I-B)-1Cx(k)+(I-B)-1f,k=0,1,2,…(3―77)

这里迭代矩阵

A=(I-B)-1C

由此看出,对方程组(3―70)用公式(3―75)作赛德尔迭代相当于对方程组(3―70)用(3―77)式作简单迭代。因此,赛德尔迭代法实际上是某种简单迭代法。6.3松驰法赛德尔迭代公式(3―75)为

x(k+1)=Bx(k+1)+Cx(k)+f

现令

Δx=x(k+1)-x(k)=Bx(k+1)+Cx(k)+f-x(k)于是

x(k+1)=x(k)+Δx(3―78)

x(k+1)可以看作在向量x(k)上加修正项Δx而得到。若在修正项的前面加上一个参数ω,便得到松驰法的迭代公式

x(k+1)=x(k)+ωΔx=(1-ω)x(k)+ω(Bx(k+1)+Cx(k)+f)k=0,1,2,…(3―79)

或者用分量形式写成(3―80)其中ω叫做松驰因子,当ω＞1时叫超松驰,ω＜1时叫低松驰,ω=1时就是赛德尔迭代法。松驰因子ω的选取直接影响到松弛法的收敛性。因为(I-ωB)-1存在,所以还可以把(3―79)改写成

x(k+1)=(I-ωB)-1((1-ω)I+ωC)x(k)+ω(I-ωB)-1f(3―81)

这是简单迭代公式,迭代矩阵为

A=(I-ωB)-1((1-ω)I+ωC)

实际上,赛德尔迭代法或松弛法都是某种简单迭代法,它们仅是迭代矩阵A不同。§7迭代法的收敛性前面介绍的例7、例8分别用简单迭代法和赛德尔迭代法都得到了较好的迭代效果。若将原方程组改写成如下等价形式:表3―3可见当k越大时,其迭代值与准确解

x1=1.1,x2=1.2,x3=1.3

相差越大,即向量序列｛x(k)｝不收敛。此例如用赛德尔迭代法迭代求解,得到的结果也令人失望。为此,很有必要讨论迭代法的收敛性。也就是研究等价方程

x=Ax+f

的收敛条件。为了达到此目的,我们首先介绍向量范数、矩阵范数等基本