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

在工程技术、自然科学和社会科学中，经常遇到的许多问题最终都可归结为解线性方程组。◎电学中的网络问题、◎最小二乘法求实验数据的曲线拟合问题、◎三次样条函数的插值问题、◎经济运行中的投入产出问题、◎大地测量、机械与建筑结构的设计计算问题等等，归结为求解线性方程组或非线性方程组的数学问题。因此，掌握线性方程组的解法对于解决实际问题是极其重要的。第5章线性方程组的数值解法

1简记为Ax=b，其中：

一般地，设b≠0,当系数矩阵A非奇异(即detA≠0)时，如下方程组有惟一解。2线性方程组的数值解法一般有两类：直接法：就是经过有限步算术运算，可求得方程组精确解的方法（若计算过程中没有舍入误差）。直接法中具有代表性的算法是高斯(Gauss)消去法。（低阶稠密矩阵）迭代法:

就是用某种极限过程去逐步逼近线性方程组的精确解的方法。也就是从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法。(一般有限步内得不到精确解)。（高阶稀疏矩阵）解线性方程组有Gram法则，但Gram法则不能用于计算方程组的解，如n＝100，1033次/秒的计算机要算10120年.3三角形方程组的解乘除法：（n＋1）n/2次运算(1)上三角方程的一般形式

5.1高斯消去法

解：

(2)下三角方程的一般形式解：4消去后两个方程中的x1得再消去最后一个方程的x2得消元结束,经过回代得解:Gauss消去法5上述求解的消元过程可用矩阵表示为:(A,b)=这是Gauss消去法的计算形式,新的增广矩阵对应的线性方程组就是上三角形方程组,可进行回代求解。现在介绍求解一般线性方程组的顺序Gauss消去法:记则,线性方程组的增广矩阵为6第一步.设,依次用乘矩阵的第1行加到第i行,得到矩阵:其中7第二步.设,依次用乘矩阵的第2行加到第i行,得到矩阵:其中8如此继续消元下去,第n-1步结束后得到矩阵:这就完成了消元过程。对应的方程组变成：对此方程组进行回代，就可求出方程组的解。9顺序Gauss消去法求解n元线性方程组的乘除运算量是：

n2-1+(n-1)2-1+…+22-1

+1+2+…+n

n=20时,顺序Gauss消去法只需3060次乘除法运算.顺序Gauss消去法通常也简称为Gauss消去法.顺序Gauss消去法中的称为主元素.主元素都不为零

Gauss消去法才能执行下去.

主元素都不为零矩阵A的各阶顺序主子式都不为零.

10步1

输入系数矩阵A，右端项b，置k:=1；步2.

消元：对k=1,2,…,n-1，计算步3.

回代对k=n-1,…,1，计算算法5.1顺序Guass消去法

11

列主元Gauss消元法解:[方法1]小主元a11回代求解x2=1.0000,x1=0.0000.计算结果相当糟糕.

原因:求行乘数时用了”小主元”0.0001作除数.例1

在F(10,4,L,U)上用Gauss消去法求解线性方程组方程组的准确解为x*=(10000/9999,9998/9999)T.12解:[方法2]若上题在求解时将1,2行交换,即回代后得到

x=(1.0000,1.0000)T与准确解x*=(9998/9999,10000/9999)T相差不大.主元a11

方法2所用的是在Gauss消去法中加入选主元过程,利用换行,避免”小主元”作除数,该方法称为Gauss列主元消去法。13由于计算机和方程组本身的限制，在Gauss消元过程中会出现零主元和小主元，而使Gauss消去法无法进行或出现较大舍入误差，为避免此类现象，可以通过行交换来选取绝对值较大的主元素.常采用的是列主元消去法和全主元消去法.给定线性方程组Ax=b,记A(1)=A,b(1)=b,列主元Gauss消去法的具体过程如下:首先在增广矩阵B(1)=(A(1),b(1))的第一列元素中,取然后进行第一步消元得增广矩阵B(2)=(A(2),b(2)).14Gauss列主元消去法因为|mi,k|≤1,列主元消去法能避免”小主元”作除数,误差分析与数值试验表明：Gauss列主元消去法“基本稳定”.再在矩阵B(2)=(A(2),b(2))的第二列元素中,取按此方法继续进行下去,经过n-1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n)).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回代求解.然后进行第二步消元得增广矩阵B(3)=(A(3),b(3)).易证,只要|A|0,列主元Gauss消去法就可顺利进行.15算法3.2（列主元Gauss消去法）：步1、输入系数矩阵A，右端项b，置k=1;步2、对于(1)按列选主元：选取l

使

(2)如果，交换(A(k),b(k))

的第k行与第l

行元素(3)

消元计算：3、回代计算16Gauss-Jordan列主元消去法

高斯消去法有消元和回代两个过程，消去的是对角线下方的元素。对消元过程稍加改变便可使方程组化为对角阵。这时求解就不需要回代。这种将主元素化为1,并用主元将其所在列的冗余元素全都消为0,即消去对角线上方与下方元素的方法称为Gauss-Jordan消去法。这时等号右端即为方程组的解。17例2用Gauss-Jordan消去法求方程组的解解方程组相应的增广矩阵为：

列选主元故得：

18定理5.1设A为非奇异矩阵,方程组Ax=I的增广矩阵为C=A

,I]

，如果对C应用高斯-约当消去法化为I

,T，则T

=A-1。的逆矩阵.

例3用Gauss-Jordan消去法求

解C=A

I

=

195.2矩阵的三角分解

顺序Gauss消去法的矩阵表示若令记则有A(2)=对于给定方阵A，如果能求出下三角方阵L和上三角方阵U，使A=LU，就说对方阵A进行了三角分解，三角分解也称LU分解。20记令若则有A(3)=如此进行下去,第n-1步得到:

A(n)=其中也就是:

A(n)=Ln-1A(n-1)=Ln-1Ln-2A(n-2)=…=Ln-1Ln-2…L2L1A(1)21所以有:A=A(1)=L1-1L2-1…Ln-1-1A(n)=LU其中L=L1-1L2-1…Ln-1-1.而且有式A=LU称为矩阵A的三角分解.A(n)=Ln-1Ln-2…L2L1A(1)令A(n)=U，即22

设A为n阶方阵，如果解Ax=b用Gauss消去法能够完成(即akk(k)不为零),则矩阵A可分解为单位下三角矩阵L和上三角矩阵U的乘积,且这种分解是唯一的。

证明只证唯一性,设有两种分解A=LU=L1U1,则有=I所以得定理5.2.1定理5.2.2约化的主元素

的充要条件是矩阵A的顺序主子式于是Ax=b

LUx=b

令

Ux=y

得23LU分解的直接推导直接利用矩阵乘法来计算

LU分解的元素比较等式两边的第一行得：u1j=a1j比较等式两边的第一列得：比较等式两边的第二行得：比较等式两边的第二列得：(j=1,…,n)(i=2,…,n)(j=2,…,n)(i=3,…,n)U

的第一行L

的第一列U

的第二行L

的第二列24LU分解算法（续）第k

步：此时U

的前k-1行和

L

的前k-1列已经求出比较等式两边的第k行得：比较等式两边的第k列得：直到第n

步，便可求出矩阵L

和U

的所有元素。(j=k,…,n

)(i=k+1,…,n

)25LU分解算法算法5.3：(LU

分解

)For

k=1,2,...,nEndForj=k,…,ni=k+1,…,n为了节省存储空间，通常用A

的绝对下三角部分来存放L(对角线元素无需存储)，用

A

的上三角部分来存放U。运算量：(n3-n)/3Doolittle三角分解26由可得这就是求解方程组Ax=b的Doolittle三角分解方法。27例利用三角分解方法解线性方程组

解因为所以先解

,得28,得解线性方程组Ax=b的Doolittle三角分解法的计算量约为n3/3,与Gauss消去法基本相同.其优点在于求一系列同系数的线性方程组Ax=bk,(k=1,2,…,m)时,可大大节省运算量.此外，还有另一种LU分解，称为Crout三角分解.再解29

A正定：A的各阶顺序主子式均大于零。对称正定矩阵的三角分解如果A为对称正定矩阵，则存在一个实的非奇异下三角矩阵L，使A=LLT

，且当限定的对角元素为正时，这种分解是唯一的。三、平方根法A是正定矩阵，即对于任意n维非零列向量x

，恒有xTAx

>0，A对称：AT=A

定理5.2.3

若A为对称正定矩阵，则A可唯一分解为A=LDLT。其中L为单位下三角，D为元素全为正数的对角阵。30U=uij=u11uij/uii111u22unn记为

A对称且分解是唯一的即记D1/2=则仍是下三角阵设A为对称正定矩阵,则有唯一分解A=LU,且ujj>0.分解A=LLT称为对称正定矩阵的Cholesky分解.如何求L？31一般有(第j列)称为平方根法，因为带了开方运算。对于k=1,2,…,n计算（一列一列算），k=1时32Cholesky分解法特点：

2.对于对称正定阵A，从可知对任意ki

有。即L

的元素不会增大，误差可控，不需选主元。缺点：存在开方运算。优点：1.可以减少存储单元。为一般矩阵LU分解计算量的一半。For

j=1,2,...,n算法5.4：(Cholesky

分解，按列)EndFori=k+1,…,n运算量：n3/6

+n2/2+n/3

33改进Cholesky分解法改进的cholesky分解法,为了避免开方运算，分解A=LDLT现确定计算L的第j列的元素。考察A的元素aij，ajj34算法5.5：(改进的Cholesky

分解法

)步1

输入对称正定矩阵A，右端项b；步2.

改进的Cholesky分解：对j=2,3,…,n，计算：步3.解下三角形方程组LDy=b步4.解上三角形方程组LTx=y改进后的LDLT分解乘除法运算量约为n3/4+O(n2),但存储量增加了。35追赶法Ax=f其中：36用矩阵的Doolitte分解，有得追：赶：计算量为5n-4次乘除法追赶法是数值稳定的,追赶法具有计算程序简单,存贮量少,计算量小的优点.37例

用追赶法求解三对角线性方程组Ax=f，其中解按公式分解有385.3向量和矩阵的范数

1．向量的范数定义1：设x

R

n，x表示定义在Rn上的一个实值函数，称之为x的范数，它满足下列性质：（3）三角不等式：即对任意两个向量x、y

R

n，恒有(1)非负性：即对一切x

R

n，x

0,x>0(2)齐次性：即对任何实数a

R,有（4）||

-x||=|

-1

|

×

||

x||=||

x||;（5）对任意的x、yÎRn，恒有

|

||

x||-||

y||

|≤

||

x-y||.有性质：39设x=(x1,x2,…,xn)T

三个常用的范数：

在Rn上定义的任一向量范数

都与范数

等价，成立。结论：Rn上定义的任何两个范数都是等价的。

对常用范数，容易验证下列不等式：

即存在正数M与m(M>m)对一切xRn，不等式40定义2：设给定Rn中的向量序列{}，即若对任何i(i=1,2,…,n)都有则向量

或者说向量序列{}依坐标收敛于向量，记为，称为向量序列{

}的极限，定理：向量序列{x(k)}依坐标收敛于x*的充要条件是向量序列依范数收敛与依坐标收敛是等价的。依范数收敛412矩阵的范数

定义3设‖‖是以n阶方阵为变量的实值函数,且满足条件:

(1)非负性:‖A‖0,且‖A‖=0当且仅当A=0(2)齐次性:‖kA‖=|k|‖A‖,kR

(3)三角不等式:‖A+B‖‖A‖+‖B‖(4)相容性:‖AB‖‖A‖‖B‖则称‖A‖为矩阵A的范数.

记A=(aij),常用的矩阵范数有:

矩阵的1-范数:‖A‖1

,也称矩阵的列范数.

矩阵的2-范数:‖A‖2

,也称为谱范数.（4'）对任意向量xRn，和任意矩阵A，有

矩阵的-范数:‖A‖

,也称为矩阵的行范数.42例设矩阵求矩阵A的范数‖A‖p,p=1,2,.

解

‖A‖1=

‖A‖=5,4，43在有限维线性空间上任何两种矩阵范数是等价的。即

m

‖A‖‖A‖

M

‖A‖,ARnn对常用范数，容易验证下列不等式都成立矩阵的范数与矩阵的特征值之间也有密切的联系.设是矩阵A的特征值,

x是对应的特征向量,则有

Ax=x利用向量和矩阵范数的相容性,则得

||‖x‖=‖x‖于是||‖A‖设n阶矩阵A的n个特征值为1,2,…,n,则称为矩阵A的谱半径.对矩阵的任何一种相容范数都有

(A)‖A‖=‖Ax‖‖A‖‖x‖44

结论：定义5.3.4设有Rnn中的矩阵序列{A(k)},其中A(k)=(aij(k))nn,若对任何i

和j(i,j=1,2,…,n)，都有则矩阵A*=(aij*)nn,称为矩阵序列{A(k)}的极限，记为依范数收敛455.4线性方程组的性态和解的误差分析方程组的精确解为x=(2,0)T.设有线性方程组Ax=b,A为非奇异方阵，x为精确解。由于A和b元素是测量或观测得到的，常带有某些观测误差，或A,b的元素是计算的结果，因而包含有舍入误差。

实际处理矩阵是A+δA，b+δb，其中δA是矩阵A的微小误差矩阵或扰动；δb是向量b的微小误差或扰动。需要研究方程组数据A,b的微小误差（扰动）对解x的影响，即考虑方程组的解y和x的差的估计问题。例5.4.146考察常数项有微小变动则方程组为即.常数项相对误差

,由于常数项微小误差，而引起解的相对误差较大，为常数项相对误差的104倍，也就是说，此方程组的解对方程组的数据A,b非常敏感。这种由于原始数据微小变化而导致解严重失真的线性方程组称为病态方程组,相应的系数矩阵称为病态矩阵.而引起解的相对误差为47设一般线性方程组Ax=b1)系数矩阵是精确的,常数项有误差δb,此时记解为x+δx

,则

A(x+δx)=b+δb于是

Aδx=δb所以又由于因此即‖A-1‖‖δb‖‖δx

‖=‖A-1δb‖‖A‖‖x‖‖b‖=‖Ax‖‖A‖‖A-1‖‖δb‖‖x‖‖δx

‖‖b‖48

2)再设b是精确的,A有误差δA,此时记解为x+δx

,则

(A+δA)(x+δx)=b则有所以

δx

=-A-1δA(x+δx)于是也就是记Cond(A)=‖A‖‖A-1‖,称为矩阵A的条件数.

Aδx

+δA(x+δx)=0

‖δx‖‖A-1‖‖δA‖‖x+δx‖49条件数的性质：

1）cond(A)≥12）cond(kA)=cond(A)(k为非零常数)3）若，则经常使用的条件数有：Condp(A)=‖A‖p‖A-1‖pp=1,2,。记Cond(A)=‖A‖‖A-1‖,称为矩阵A的条件数.是良态方程组，或A是良态的。例5.4.2

对于例5.4.1的方程组，计算A的条件数。.50由于计算条件数运算量较大,实际计算中若遇到下述情况之一,方程组就有可能是病态的:行列式很大或很小（如某些行、列近似相关）；元素间相差大数量级，且无规则；主元消去过程中出现小主元；特征值相差大数量级。Hilbert矩阵病态方程组的计算1)用双精度或更高精度计算2)采用数值稳定性较好的算法，如全主元法3)用预处理方法。515.5

解线性方程组的迭代法迭代法的基本思想是构造一串收敛到解的序列，即建立一种从已有近似解计算新的近似解的规则。思路将Ax=b等价改写为x=Bx+f,建立迭代格式：x(k+1)=Bx(k)+f,从初值x(0)出发，得到迭代序列{x(k)}(k=0,1,...).如果其迭代序列的极限存在,即则称迭代法收敛,

x*即是线性方程组Ax=b的解.否则，称迭代法发散.矩阵B称为迭代矩阵，x(k)为第k次迭代近似解，称e(k)=x(k)-x*为第k次迭代误差向量.研究内容：如何建立迭代格式？

迭代序列收敛速度？

迭代序列的收敛条件？

近似解的误差估计？525.5.1Jacobi迭代法若系数矩阵非奇异，且

(i=1,2,…,n)，将方程组改成53然后写成迭代格式（5.5.5）（5.5.6）式也可以简单地写为Jacobi迭代格式。54写成矩阵形式：A=－L－UDBJacobi迭代阵(5.5.8)Jacobi迭代格式简单，每迭代一次只需计算一次矩阵与向量乘法，计算机中只需要两组工作单元用来保存x(k)及x(k+1),用来控制迭代终止。55算法5.7（Jacobi迭代法）1.输入矩阵A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次数N，置k=0;

2.对i=1,2,…,n,计算3.若‖x(k+1)-x(k)‖<ε,则停算，输出x(k+1)为近似解；4.若k<N,置k=k+1,转步2；否则，停算，输出迭代失败信息.

56

例5.5.1

用雅可比迭代法解方程组解：雅可比迭代格式为取经14次迭代,近似解为误差为57为了加快收敛速度，同时为了节省计算机的内存，作如下的改进：

5.5.2Gauss-Seidel迭代法逐一写出来即为每算出一个分量的近似值，立即用到下一个分量的计算中去，即用迭代格式：58…………写成矩阵形式：BGauss-Seidel

迭代阵(5.5.11)(5.5.13)只保存一组向量59算法5.8（Gauss-Seidel迭代法）：1.输入矩阵A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次数N，置k=0;

2.计算对i=2,…,n-1,计算3.若‖x(k+1)-x(k)‖<ε,则停算，输出x(k+1)为近似解；4.若k<N,置k=k+1,转步2；否则，停算，输出迭代失败信息.

60

例用Gauss-Seidel迭代法解方程组取经6次迭代,近似解为误差为解：Gauss-Seidel迭代格式61松弛（5.5.18）是松驰因子(0<<2),当0<<1时叫低松弛，>1时叫超松弛，=1时，就是Gauss-Seidel迭代法。5.3

超松弛法逐次超松弛迭代法（SuccessiveOverRelaxationMethod.简称为SOR法）是Gauss-Seidel法的一种加速方法。这种方法将前一步的迭代结果xi(k)与Gauss-Seidel法的迭代值,

适当进行线性组合，以构成一个收敛速度较快的近似解序列。62对Gauss－Seidel迭代格式有故SOR的迭代格式用矩阵语言描述为：Bωfω63算法5.9（SOR迭代法）1.输入矩阵A=(aij),b=(b1,…,bn),初始向量x(0)=(x1(0),…,xn(0)),精度要求ε，最大的迭代次数N，参数ω；置k=0;

2.计算对i=2,…,n-1,计算3.若‖x(k+1)-x(k)‖<ε,则停算，输出x(k+1)为近似解；4.若k<N,置k=k+1,转步2；否则，停算，输出迭代失败信息.

64例5.5.3用SOR法法解方程组初值x(0)=0，写出计算格式。解：SOR迭代格式为

选取ω=1.3,第11次迭代结果为：

且满足:65记误差向量e(k)=x(k)-x*,则迭代法收敛就是e(k)0.

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

k=0,1,2,…

x*=Bx*+f

所以

e(k+1)=Be(k),

k=0,1,2,…递推可得

e(k)=Bke(0),

k=0,1,2,…可见,当k时,

e(k)0Bk

O.对任意初始向量x(0),迭代法收敛

(B)<1.定理5.6.15.6迭代法的收敛性分析定理5.6.2若‖B‖=q<1,则对任意x(0),迭代法收敛,而且

由定理5.3.2知，

，再由定理5.6.1知，迭代法收敛。证明：先证收敛性。66证由于

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

,x(k)=Bx(k-1)+f,且

x*=Bx*+f；

所以

x(k+1)-x(k)=B(x

(k)-x(k-1)),x(k+1)–x*=B(x

(k)–x*)于是有

‖x(k+1)-x(k)‖q‖x

(k)-x(k-1)‖,‖x(k+1)–x*‖q‖x

(k)–x*‖

‖x(k+1)-x(k)‖=‖(x

(k+1)–x*)-(x(k)–x*)‖

‖x

(k)–x*‖-‖x(k+1)–x*‖(1-q)‖x(k)–x*‖所以定理5.6.2只是收敛的充分条件,并不必要,如则‖M‖1=1.2,‖M‖=1.3,‖M‖2=1.09,但(M)=0.8<1,所以迭代法是收敛的.可见,q越小收敛越快,且当‖x(