解线性方程组的直接方法

解线性方程组的直接方法2023/6/61第一页，共一百五十五页，编辑于2023年，星期二这一章讨论线性方程组引言与预备知识的数值解法.2023/6/62第二页，共一百五十五页，编辑于2023年，星期二在自然科学和工程技术中,很多问题归结为解线性方程组.例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵(例如，阶数不超过150).另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多).一、引言2023/6/63第三页，共一百五十五页，编辑于2023年，星期二有的问题的数学模型中虽不直接表现为含线性方程组,但它的数值解法中将问题“离散化”或“线性化”为线性方程组.因此线性方程组的求解是数值分析课程中最基本的内容之一.

关于线性方程组的解法一般有两大类：2023/6/64第四页，共一百五十五页，编辑于2023年，星期二但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解.本章将阐述这类算法中最基本的和具有代表性的算法就是高斯消去法,以及它的某些变形和应用.这类方法是解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组(例如，大型带状方程组)的有效方法.

1.直接法经过有限次的算术运算,可以求得方程组的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法对高阶方程组计算量太大,不是一种实用的算法.2023/6/65第五页，共一百五十五页，编辑于2023年，星期二就是用某种极限过程去逐步逼近方程组精确解的方法.迭代法具有计算机的存储单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛条件和收敛速度问题.迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)的重要方法(见第6章).为了讨论线性方程组数值解法，需复习一些基本的矩阵代数知识.

2.迭代法2023/6/66第六页，共一百五十五页，编辑于2023年，星期二二、向量和矩阵

基本概念：用Rm×n表示全部m×n实矩阵的向量空间；用Cm×n表示全部m×n复矩阵的向量空间.(由数排成的矩阵表，称为m行n列矩阵).2023/6/67第七页，共一百五十五页，编辑于2023年，星期二其中aj为A的第j列的m维列向量.同理(m维列向量).其中biT为A的第i行的n维行向量.2023/6/68第八页，共一百五十五页，编辑于2023年，星期二

矩阵的基本运算：

(1)矩阵加法

(2)矩阵与标量的乘法

(3)矩阵与矩阵的乘法

(4)转置矩阵2023/6/69第九页，共一百五十五页，编辑于2023年，星期二(5)单位矩阵其中

(6)非奇异矩阵则称B是A的逆矩阵，记为A-1，且(A-1)T=(AT)-1.如果A-1存在，则A称为非奇异矩阵.如果A、B均为非奇异矩阵，则有(AB)-1=B-1A-1.2023/6/610第十页，共一百五十五页，编辑于2023年，星期二

(7)矩阵的行列式设A∈Rn×n，则A的行列式可按任一行(列)展开,其中Aij为aij的代数余子式，Aij=(-1)i+jMij，为Mij元素aij的余子式.

行列式性质：2023/6/611第十一页，共一百五十五页，编辑于2023年，星期二三、矩阵的特征值与谱半径

定义1设A=(aij)Rn×n，若存在数(实数或复数)和非零向量x=(x0,x1,…,xn)TRn，使Ax=x,(0.1)则称为A的特征值,x为A对应的特征向量,A的全体特征值称为A的谱,记作σ(A),即σ(A)={1,2,…,n}.记称为A的谱半径.由(0.1)式可知可使齐次线性方程组(I-A)x=0有非零解，故系数行列式det(I-A)=0，记2023/6/612第十二页，共一百五十五页，编辑于2023年，星期二p()称为矩阵A的特征多项式，方程(0.3)称为A的特征方程.因为n次代数方程p()在复数域中有n个根1,2,…,n，故2023/6/613第十三页，共一百五十五页，编辑于2023年，星期二由(1.3)式中的行列式展开可得故A=(aij)Rn×n的n个特征值1,2,…,n是它的特征方程(1.3)的n个根，并有称trA为A的迹.及2023/6/614第十四页，共一百五十五页，编辑于2023年，星期二

A的特征值和特征向量x还有以下性质：⑴AT与A有相同的特征值及特征向量x.⑵若A非奇异,则A-1的特征值为-1,特征向量为x.⑶相似矩阵B=S-1AS有相同的特征多项式.

例1求矩阵的特征值及谱半径

解矩阵A的特征方程为得A的特征值为1=2=2,3=7，谱半径(A)=7.2023/6/615第十五页，共一百五十五页，编辑于2023年，星期二特殊矩阵设A=(aij)∈Rn×n.

(1)对角矩阵如果当i≠j时，aij=0.

(2)三对角矩阵如果当|i-j|>1时，aij=0.

(3)上三角矩阵如果当i>j时，aij=0.

(4)上海森伯格(Hessenberg)矩阵如果当i>j+1时，aij=0.

(5)对称矩阵如果AT=A.

(6)埃尔米特(Hermit)矩阵设A∈Cn×n，如果当AH=A(AH是A的共轭转置矩阵).2023/6/616第十六页，共一百五十五页，编辑于2023年，星期二

(7)对称正定矩阵如果AT=A且对任意非零向量，x∈Rn，(Ax,x)=xTAx>0.

(8)正交矩阵如果A-1=AT.(9)酉矩阵设A∈Cn×n，如果A-1=AH.(10)初等置换阵由单位矩阵I交换第i行与第j行(或第i列与第j列)，得到的矩阵记为Iij，且.

IijA=B(为交换A第i行与第j行得到的矩阵)；

AIij=C(为交换A第i列与第j列得到的矩阵).

(11)置换阵由初等置换阵的乘积得到的矩阵.2023/6/617第十七页，共一百五十五页，编辑于2023年，星期二

定理1设A∈Rn×n，则下述命题等价：(1)对任何b∈Rn，方程组Ax=b有唯一解.(2)齐次方程组Ax=0只有唯一解零解x=0.(3)det(A)≠0.(4)A-1存在.(5)A的秩rank(A)=n.2023/6/618第十八页，共一百五十五页，编辑于2023年，星期二

定理2设A∈Rn×n为对称正定矩阵，则(1)A为非奇异矩阵，且A-1亦是对称正定矩阵.(2)设Ak为A的顺序主子阵，则Ak(k=1,2,,n)亦是对称正定矩阵，其中(3)A的特征值i>0

(i=1,2,,n).(4)A的顺序主子式都大于零，即Ak的行列式det(Ak)>0(k=1,2,,n).2023/6/619第十九页，共一百五十五页，编辑于2023年，星期二

定理3设A∈Rn×n为对称矩阵，如果det(Ak)>0

(k=1,2,,n)，或A的特征值i>0

(i=1,2,,n)，则A为对称正定矩阵.有重特征值的矩阵不一定相似于对角矩阵，那么一般n阶矩阵A在相似变换下能简化到什么形状.

定理4(Jordan标准型)设A为n阶矩阵，则存在一个非奇异n阶矩阵P使得2023/6/620第二十页，共一百五十五页，编辑于2023年，星期二为若当(Jordan)块.其中2023/6/621第二十一页，共一百五十五页，编辑于2023年，星期二(1)当A的若当标准型中所有若当块Ji均为一阶时，此标准型变为对角矩阵.(2)如果A的特征值各不相同，则其若当标准型必为对角矩阵diag(1,

2,,n).2023/6/622第二十二页，共一百五十五页，编辑于2023年，星期二6.1高斯消去法

本节介绍高斯消去法(逐次消去法)及消去法和矩阵三角分解之间的关系.虽然高斯消去法是一种古老的求解线性方程组的方法(早在公元前250年我国就掌握了解方程组的消去法)，但由它改进、变形得到的选主元素消去法、三角分解法仍然是目前计算机上常用的有效方法.我们在中学学过消去法,高斯消去法就是它的标准化的、适合在计算机上自动计算的一种方法.2023/6/623第二十三页，共一百五十五页，编辑于2023年，星期二6.1.1高斯消去法设有线性方程组或写成矩阵形式简记为Ax=b.2023/6/624第二十四页，共一百五十五页，编辑于2023年，星期二解为当m=n时，对三角形方程组的解非常容易的.如:上三角矩阵所对应的线性方程组2023/6/625第二十五页，共一百五十五页，编辑于2023年，星期二下三角矩阵所对应的线性方程组计算量（乘除法的主要部分）都为n2/2.解为因此，我们将一般的线性方程组化成等价的三角形方程组来求解.2023/6/626第二十六页，共一百五十五页，编辑于2023年，星期二首先举一个简单的例子来说明消去法的基本思想.

例2用消去法解线性方程组

解第1步，将方程(1.2)乘上-2加到方程(1.4)上去，消去(1.4)中的未知数x1，得到第2步，将方程(1.3)加到方程(1.5)上去，消去(1.5)中的未知数x2，得到与原方程组等价的三角形方程组2023/6/627第二十七页，共一百五十五页，编辑于2023年，星期二显然，方程组是(1.6)是容易求解的，解为上述过程相当于对方程的增广阵做初等行变换其中ri用表示矩阵的第i行.2023/6/628第二十八页，共一百五十五页，编辑于2023年，星期二由此看出，用消去法解方程组的基本思想是用逐次消去未知数的方法把原方程组Ax=b化为与其等价的三角形方程组，而求解三角形方程组可用回代的方法求解.换句话说，上述过程就是用初等行变换将原方程组系数矩阵化为简单形式(上三角矩阵)，从而求解原方程组(1.1)的问题转化为求解简单方程组的问题.或者说，对系数矩阵A施行一些行变换(用一些简单矩阵左乘A)将其约化为上三角矩阵.这就是高斯消去法.2023/6/629第二十九页，共一百五十五页，编辑于2023年，星期二下面讨论求解一般线性方程组的高斯消去法.由2023/6/630第三十页，共一百五十五页，编辑于2023年，星期二将(1.1)记为A(1)x=b(1)，其中2023/6/631第三十一页，共一百五十五页，编辑于2023年，星期二(1)第1步(k=1).

设a11(1)0，首先计算乘数

mi1=ai1(1)/a11(1)(i=2,3,…,m).用-mi1乘(1.1)的第一个方程，加到第i个(i=2,3,,m)方程上，消去(1.1)的从第二个方程到第m个方程中的未知数x1，得到与(1.1)等价的方程组(1.7)2023/6/632第三十二页，共一百五十五页，编辑于2023年，星期二简记为A(2)x=b(2)，其中A(2),b(2)的元素计算公式为(2)第k次消元(k=1,2,,s,s=min(m-1,n)).

设上述第1步，

，第k-1步消元过程计算已经完成，即已计算好与(1.1)等价的方程组，简记为A(k)x=b(k)，其中2023/6/633第三十三页，共一百五十五页，编辑于2023年，星期二

设akk(k)0，计算乘数

mik=aik(k)/akk(k)

(i=k+1,,m).用-mik乘(1.8)的第k个方程加到第i个(i=k+1,,m)方程上，消去从第k+1个方程到第m个方程中的未知数xk，得到与(1.1)等价的方程组A(k+1)x=b(k+1)，2023/6/634第三十四页，共一百五十五页，编辑于2023年，星期二其中A(k+1),b(k+1)的元素计算公式为，显然A(k+1)中从第1行到第k行与A(k)相同.

(3)继续上述过程，且设aii(i)0(i=1,2,,s)，直到完成第s步消元计算.最后得到与原方程组等价的简单方程组A(s+1)x=b(s+1)，其中A(s+1)为上阶梯.2023/6/635第三十五页，共一百五十五页，编辑于2023年，星期二

特别当m=n时，与原方程组等价的简单方程组为A(n)x=b(n)，即由(1.1)约化为(1.10)的过程称为消元过程.

如果A是非奇异矩阵，且akk(k)0(k=1,2,,n-1)，求解三角形方程组(1.10)，得到求解公式(1.10)的求解过程(1.11)称为回代过程.2023/6/636第三十六页，共一百五十五页，编辑于2023年，星期二

注意：设Ax=b，其中A∈Rn×n为非奇异矩阵，如果a11(1)=0，由于A为非奇异矩阵，所以A的第1列一定有元素不等于零，例如al10，于是可交换两行元素(即r1↔rl)，将al1

调到(1,1)位置，然后进行消元计算，这时A(2)右下角矩阵为n-1阶非奇异矩阵，继续这过程，高斯消去法照样可进行计算.总结上述讨论即有2023/6/637第三十七页，共一百五十五页，编辑于2023年，星期二

定理5设Ax=b，其中A∈Rn×n.

(1)如果akk(k)0(k=1,2,…,n)，则可通过高斯消去法将Ax=b约化为等价的三角形方程组(1.10)，且计算公式为

①消元计算(k=1,2,,n-1)2023/6/638第三十八页，共一百五十五页，编辑于2023年，星期二

②回代计算

(2)如果A为非奇异矩阵，则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为等价的三角形方程组(1.10).2023/6/639第三十九页，共一百五十五页，编辑于2023年，星期二

算法1(高斯算法)设A∈Rm×n(m>1),s=min(m-1,n)，如果akk(k)0(k=1,2,,s)，本算法用高斯方法将A约化为上三角形矩阵，且A(k)覆盖A，乘数mik覆盖aik.

对于k=1,2,,s

(1)如果akk=0，则计算停止；

(2)对于i=k+1,,m

(a)aik←mik=aik/akk

(b)对于j=k+1,,n

aij←aij-mik·akj

.2023/6/640第四十页，共一百五十五页，编辑于2023年，星期二

显然，算法1第k步需要m-k次除法，(m-k)(n-k)次乘法，因此，本算法(从第1步到第s步消元计算总的计算量)大约需要s3/3-(m-k)s2/2+mns次乘法运算(对相当大的s).当m=n时,总共大约需要n3/3次乘法运算.

数akk(k)在高斯消去法中有着突出的作用，称为约化的主元素(简称主元).2023/6/641第四十一页，共一百五十五页，编辑于2023年，星期二

算法2(回代算法)设Ux=b，其中U=(uij)∈Rn×n为非奇异上三角矩阵，本算法计算Ux=b的解.

对于i=n,n-1,,1

(1)xi←bi

(2)对于j=i+1,,n

xi←xi-uij·xj(i>n不算)

(3)xi←xi/uii

这个算法需要n(n-1)/2次乘除法运算.2023/6/642第四十二页，共一百五十五页，编辑于2023年，星期二例子解方程组解：消元回代得2023/6/643第四十三页，共一百五十五页，编辑于2023年，星期二

高斯消去法对于某些简单的矩阵可能会失败，例如

由此，需要对算法1进行修改，首先研究原来矩阵A在什么条件下才能保证akk(k)0(k=1,2,,s).下面的定理给出了这个条件.2023/6/644第四十四页，共一百五十五页，编辑于2023年，星期二

定理6约化的主元素aii(i)0(i=1,2,,k)的充要条件是方阵A的顺序主子式Di0(i=1,2,,k)，即

证明首先利用归纳法证明充分性.显然,当k=1时，结论成立，假设结论对k-1是成立的，对k由条件设Di0(i=1,2,,k)，于是由归纳法假设我们有aii(i)0

(i=1,2,,k-1)，可用高斯消去法将A(1)约化到A(k)，即2023/6/645第四十五页，共一百五十五页，编辑于2023年，星期二利用行列式的性质，我们有2023/6/646第四十六页，共一百五十五页，编辑于2023年，星期二

由条件有Di0(i=1,2,,k)，利用(1.13)式，则有aii(i)0(i=1,2,,k)，定理6的对k成立.

必要性，由条件aii(i)0(i=1,2,,k)，利用(1.13)式亦可推出Di0(i=1,2,,k).

定理6得证.

推论如果A的顺序主子式Di0(i=1,2,,n-1),则2023/6/647第四十七页，共一百五十五页，编辑于2023年，星期二6.1.2矩阵的三角分解

下面我们借助矩阵理论进一步对消去法做些分析，从而建立高斯消去法与矩阵因式分解的关系.

设(1.1)的系数矩阵A∈Rn×n的各顺序主子式均不为零，对于A施行初等行变换相当于用初等矩阵左乘A，于是对(1.1)施行第1次消元后化为(1.7)，这时A(1)化为A(2)，b(1)化为b(2)，即L1A(1)=A(2),L1b(1)=b(2)，2023/6/648第四十八页，共一百五十五页，编辑于2023年，星期二其中

一般第k步消元，A(k)化为A(k+1)，b(k)化为b(k+1)，相当于LkA(k)=A(k+1),Lkb(k)=b(k+1)，2023/6/649第四十九页，共一百五十五页，编辑于2023年，星期二其中重复以上过程，最后得到(1.14)2023/6/650第五十页，共一百五十五页，编辑于2023年，星期二将上三角矩阵A(n)记为U，由(1.14)得到其中为单位下三角矩阵.2023/6/651第五十一页，共一百五十五页，编辑于2023年，星期二

这就是说，高斯消去法实质上产生了一个将A分解为两个三角形矩阵相乘的因式分解，于是我们得到如下重要定理，它在解方程组的直接法中起着重要作用.

定理7(矩阵的LU分解)设A为n阶矩阵，如果A的顺序主子式Di0(i=1,2,,n-1)，则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，即A=LU,且这种分解是唯一的.2023/6/652第五十二页，共一百五十五页，编辑于2023年，星期二

证明根据以上高斯消去法的矩阵分析，A=LU的存在性已经得到证明，现仅在A为非奇异矩阵的假定下来证明唯一性，当A为奇异矩阵的情况留作练习，设A有两种分解为A=LU=L1U1，其中L,L1为单位下三角矩阵，U,U1为上三角矩阵.由于L-1,U1-1存在，故有L-1L1=UU1-1.上式右边为上三角矩阵，左边为单位下三角矩阵，从而上式两边都必须等于单位矩阵I，故有

L=L1，U=U1.证毕.2023/6/653第五十三页，共一百五十五页，编辑于2023年，星期二

例3对于例2，系数矩阵由高斯消去法m21=0，m31=2，m32=-1，故从而得到求矩阵行列式的计算公式为2023/6/654第五十四页，共一百五十五页，编辑于2023年，星期二6.1.3列主元消去法

由高斯消去法知道,在消元过程中可能有akk(k)＝0的情况，这时消去法将无法进行；即使主元素akk(k)0但很小时，用其作除数，会导致其它元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠.先看一个例子.

高斯消去法也称主元素消去法(条件detA0)

即当akk(k)=0时，高斯消元法无法进行；或

|akk(k)|<<1时，带来舍入误差的扩散。（小主元）2023/6/655第五十五页，共一百五十五页，编辑于2023年，星期二解

(方法1)高斯消去法(用4位有效数字)例4(小主元产生了大误差)2023/6/656第五十六页，共一百五十五页，编辑于2023年，星期二显然计算得到的解x1=0是一个很坏的结果，不能作为方程组的近似解.其原因是我们在消元计算时用了小主元0.00001。使得约化后的方程组元素数量级大大增长，经再舍入使得在计算(0.3),(1.2),(1.3)元素时发生了严重的相消情况，因此经消元后得到的三角形方程组就相当不准确了，所以产生了很大的误差，结果是不能用的.2023/6/657第五十七页，共一百五十五页，编辑于2023年，星期二

(方法2)列主元消元法，先交换行2023/6/658第五十八页，共一百五十五页，编辑于2023年，星期二

这个例子告诉我们，在采用高斯消去法解方程组时，小主元可能产生麻烦，故应避免绝对值小的主元素akk(k)，对一般矩阵来说，最好每一步选取系数矩阵(或消元后的低阶矩阵)中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数值稳定性，这就是全主元素消去法，在选主元时要花费较多机器时间，目前主要使用的是列主元素消去法，下面介绍列主元素消去法，并假定线性方程组(1.1)的系数矩阵A∈Rn×n为非奇异的.2023/6/659第五十九页，共一百五十五页，编辑于2023年，星期二设方程组(2.1)的增广矩阵为首先在A的第1列中选取绝对值最大的元素作为主元素，例如2023/6/660第六十页，共一百五十五页，编辑于2023年，星期二然后交换B的第1行与第i1行，经第1次消元计算得(A|b)→(A(2)|b(2))重复上述过程，设已完成第k-1步的选主元素，交换两行及消元计算，(A|b)约化为其中A(k)的元素仍记为aij，b(k)的元素仍记为bi.2023/6/661第六十一页，共一百五十五页，编辑于2023年，星期二第k步的选主元素(在A(k)右下角方阵的第1列内选)，即确定ik，使交换(A(k)|b(k))第k行与ik行的元素，再进行消元计算，最后将原方程组化为

(k=1,2,,n-1)2023/6/662第六十二页，共一百五十五页，编辑于2023年，星期二

回代求解

算法1(列主元素消去法)设Ax=b.本算法用A的具有行交换的列主元素消去法，消元结果冲掉A，乘数mij冲掉aij，计算解x冲掉常数项b，行列式存放在det中.

1.det←1

2.对k=1,2,,n-1

(1)按列选主元2023/6/663第六十三页，共一百五十五页，编辑于2023年，星期二

(2)如果aik,k=0，则计算停止(det(A)=0)

(4)消元计算对于i=k+1,,n

①aik←mik=aik/akk

②对于j=k+1,,n

aij←aij-mik·akj

(3)如果ik=k，则转(4)

换行akj↔aik,j(j=k,k+1,,n)

bk↔bik

det←-det

③bi←bi-mik·bk

(5)det←akk·

det2023/6/664第六十四页，共一百五十五页，编辑于2023年，星期二

3.如果ann=0，则计算停止(det(A)=0)

4.回代求解

(1)bn←bn/ann

(2)对于i=n-1,,2,1

5.det←ann·

det

例4的方法2用的就是列主元素消去法.2023/6/665第六十五页，共一百五十五页，编辑于2023年，星期二

下面用矩阵运算来描述解(1.1)的列主元素消去法.列主元素消去法为其中Lk的元素满足|mik|≤1

(k=1,2,,n-1)，是初等置换矩阵.

利用(1.15)得到(1.15)2023/6/666第六十六页，共一百五十五页，编辑于2023年，星期二简记为其中下面就n=4来考察一下矩阵其中可知亦为单位下三角矩阵，其元素的绝对值不超过1.记(1.16)2023/6/667第六十七页，共一百五十五页，编辑于2023年，星期二由(1.16)得到PA=LU.其中P为排列矩阵，L为单位下三角矩阵，U为上三角矩阵.这说明对(1.1)应用列主元消去法相当于对(A|b)先进行一系列行交换后对PAx=Pb再应用高斯消去法.在实际计算中我们只能在计算过程中做行的交换.总结以上的讨论我们有2023/6/668第六十八页，共一百五十五页，编辑于2023年，星期二

定理8(列主元素的三角分解定理)如果A为非奇异矩阵，则存在排列矩阵P使PA=LU，其中L为单位下三角矩阵，U为上三角矩阵.在编程过程中，L元素存放在数组A的下三角部分，U元素存放在数组A的上三角部分，由纪录主行的整型数组Ip(n)可知P的情况.2023/6/669第六十九页，共一百五十五页，编辑于2023年，星期二消元法的应用1.求行列式高斯消元法2.求逆矩阵(用高斯-若当消去法)2023/6/670第七十页，共一百五十五页，编辑于2023年，星期二例子分别用高斯消元法和列主元消元法求矩阵的行列式.解：高斯消元法大数“吃掉”了小数！2023/6/671第七十一页，共一百五十五页，编辑于2023年，星期二列主元消元法2023/6/672第七十二页，共一百五十五页，编辑于2023年，星期二

全主元素消去法介绍在消元计算过程中，若每次选主元素不局限在方框列中，而在整个主子矩阵中选取,便称为全主元高斯消去法,此时增加的步骤为:

①，确定ik

,jk,若=0，给出|Ａ|=0的信息，退出计算.2023/6/673第七十三页，共一百五十五页，编辑于2023年，星期二

行交换:(kjn)

列交换:(kin)②作如下行交换和列交换值得注意的是,在全主元的消去法中，由于进行了列交换，x各分量的顺序已被打乱.因此必须在每次列交换的同时，让机器“记住”作了一次怎样的交换，在回代解出后将x各分量换回原来相应的位置，这样增加了程序设计的复杂性.此外作①步比较大小时，全主元消去法将耗用更多的机时.2023/6/674第七十四页，共一百五十五页，编辑于2023年，星期二但全主元消去法比列主元消去法数值稳定性更好一些.实际应用中，这两种选主元技术都在使用.选主元素的高斯消元法是一种实用的算法，它可以应用于任意的方程组Ax=b，只要|A|≠0.2023/6/675第七十五页，共一百五十五页，编辑于2023年，星期二6.2.追赶法在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条插值函数，都会要求解系数矩阵为对角占优的三对角线方程组Ax=f，即:(2.1)2023/6/676第七十六页，共一百五十五页，编辑于2023年，星期二其中,当|i-j|>1时,aij=0,且满足如下的对角占优条件:①|b1|>|c1|>0;②|bi|≥|ai|+|ci|,ai,ci≠0,(i=2,3,…,n-1).③|bn|>|an|>0.我们利用矩阵的直接三角分解法来推导解三对角线方程组(2.1)的计算公式.由系数阵A的特点，可以将A分解为两个三角矩阵的乘积，即A=LU,其中取L下三角阵,取U为单位上三角阵,这样求解方程组Ax=f的方法称为追赶法.设2023/6/677第七十七页，共一百五十五页，编辑于2023年，星期二其中为待定系数，比较(2.2)两边即得(2.2)2023/6/678第七十八页，共一百五十五页，编辑于2023年，星期二下面我们用归纳法证明(2.4)对i=1是成立的.现设(2.4)对i-1成立，求证对i亦成立.(2.3)(2.3)(2.4)2023/6/679第七十九页，共一百五十五页，编辑于2023年，星期二由归纳法假设0<|i-1|<1，又由(2.4)及A的假设条件有(2.4)2023/6/680第八十页，共一百五十五页，编辑于2023年，星期二简记为Ax=f,A满足条件（1）（2）（3）小结：2023/6/681第八十一页，共一百五十五页，编辑于2023年，星期二用归纳法可以证明，满足上述条件的三对角线性方程组的系数矩阵A非奇异，所以可以利用矩阵的直接三角分解法来推导解三对角线性方程组的计算公式，用克洛特分解法，将A分解成两个三角阵的乘积，设A=LU：

按乘法展开

则可计算

可依次计算当，由上式可惟一确定L和U2023/6/682第八十二页，共一百五十五页，编辑于2023年，星期二例6.9用追赶法求解三对角方程组

解由Ly=f

解出y又由Ux=y解出x2023/6/683第八十三页，共一百五十五页，编辑于2023年，星期二求解Ax=f等价于求解两个三角形方程组.①Ly=f,求y;②Ux=y,求x.从而得到解三对角线方程组的追赶法公式.(1).计算{i},{i}的递推公式小结：2023/6/684第八十四页，共一百五十五页，编辑于2023年，星期二我们将计算系数及的过程称为追的过程,将计算方程组的解的过程称为赶的过程.合起来就是追赶法.追赶法公式实际上就是把高斯消去法用到求解三对角线方程组上去的结果.这时由于A特别简单,因此使得求解的计算公式非常简单,而且计算量仅为5n-4次乘除法,而另外增加一个方程组Ax=f2仅增加3n-2次乘除法运算,易见追赶法的计算量是比较小的.2023/6/685第八十五页，共一百五十五页，编辑于2023年，星期二总结上述讨论有下面定理

定理11设有三对角线方程组Ax=f,其中A满足条件(1),(2),(3),则A为非奇异矩阵且追赶法计算公式中的{i},{i}满足由定理11的(1),(2)说明追赶法计算公式中不会出现中间结果数量级的巨大增长和舍入误差的严重累积.2023/6/686第八十六页，共一百五十五页，编辑于2023年，星期二补充：矩阵三角分解法高斯消去法有很多变形，有的是高斯消去法的改进、改写，有的是用于某一类特殊矩阵的高斯消去法的简化.下面我们将介绍矩阵的直接三角分解法，解特殊方程组用的平方根法及追赶法.

定义

如果L为单位下三角阵，U为上三角阵，则称A=LU为杜里特尔(Doolittle)分解；如果L为下三角阵，U为单位上三角阵，则称A=LU为克劳特(Crout)分解.2023/6/687第八十七页，共一百五十五页，编辑于2023年，星期二1、直接三角分解法(LU分解)

在6.1已经通过高斯消去法得到一个将A分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，A=LU，其中并由定理7得到这种分解是唯一的.2023/6/688第八十八页，共一百五十五页，编辑于2023年，星期二

将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L、U元素的递推公式，而不需要任何中间步骤，这就是所谓直接三角分解法.一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求两个三角形方程组.①Ly=b，求y；②Ux=y，求x.

1.不选主元的三角分解法设A为非奇异矩阵，且有分解式A=LU，其中L为单位下三角矩阵，U为上三角矩阵，即2023/6/689第八十九页，共一百五十五页，编辑于2023年，星期二其中2023/6/690第九十页，共一百五十五页，编辑于2023年，星期二a11a12a1k

a1n

u11u12u1k

u1n

第1框a21a22a2k

a2n

l21u22u2k

u2n第2框

ak1ak2

akk

akn

lk1lk2

ukk

ukn第k框

an1an2ankann

ln1ln2

lnk

unn第n框比较式A=LU

两端的元素,按下图所示顺序逐框进行，先求ukj，后求lik.由第一框可得2023/6/691第九十一页，共一百五十五页，编辑于2023年，星期二得假设前k-1框元素已求出，则由2023/6/692第九十二页，共一百五十五页，编辑于2023年，星期二有了矩阵A

的LU分解计算公式,解线性方程组Ax=b

就转化为依次解下三角方程组

Ly=b

与上三角方程组Ux=y

其计算公式如下：

Ly=b

Ux=y2023/6/693第九十三页，共一百五十五页，编辑于2023年，星期二例6.10求矩阵解用紧凑形式计算的LU(Doolittle)分解.

u11=2

u12=1u13=4

2023/6/694第九十四页，共一百五十五页，编辑于2023年，星期二所以u11=2u12=1u13=4

得行列式2023/6/695第九十五页，共一百五十五页，编辑于2023年，星期二由于在计算机实现时当uri计算好后ari就不用了，因此计算好L,U的元素后就存放在A的相应位置.例如最后在存放A的数组中得到L,U的元素.由直接三角分解计算公式，需要计算形如∑aibi的式子，可采用“双精度累加”，以提高精度.2023/6/696第九十六页，共一百五十五页，编辑于2023年，星期二如果已经实现了A=LU的分解计算，且L,U保存在A的相应位置，则用直接三角分解法解具有相同系数的方程组Ax=(b1,b2,…,bm)是相当方便的，每解一个方程组Ax=bj仅需要增加次乘除法运算.矩阵A的分解公式(3.2),(3.3)又称为杜里特尔(Doolittle)分解.直接分解法大约需要n3/3次乘除法，和高斯消去法计算量基本相同.2023/6/697第九十七页，共一百五十五页，编辑于2023年，星期二

2.选主元的三角分解法从直接三角分解公式了看出当urr=0时计算将中断，或者当urr绝对值很小时，按分解公式直接计算可能引起舍入误差的积累.但如果A非奇异，我们可通过交换A的行实现矩阵PA的LU分解，因此可采用与列主元消去法类似的方法(可以证明下述方法与列主元消去法等价)，将直接三角分解法修改为(部分)选主元的三角分解法.2023/6/698第九十八页，共一百五十五页，编辑于2023年，星期二设第r-1步分解已完成，这时有设r步分解需用到(3.2)式及(3.3)式，为了避免用小的数urr作除数，引进量2023/6/699第九十九页，共一百五十五页，编辑于2023年，星期二于是有取交换A的r行与ir行元素，将调到(r,r)位置(将(i,j)位置的新元素仍记为lij及aij

),于是有|lir|1(i=r+1,…,n).由此再进行第r步分解计算.2023/6/6100第一百页，共一百五十五页，编辑于2023年，星期二算法2(选主元的三角分解法)设Ax=b，其中A为非奇异矩阵.本算法采用选主元的三角分解法，用PA的三角分解冲掉A，用整型数组Ip(n)记录主行，解x存放在b内.

1.对于r=1,2,…,n

(1)计算si

(2)选主元2023/6/6101第一百零一页，共一百五十五页，编辑于2023年，星期二上述计算过程完成后就实现了PA的LU分解，且U保存在A的上三角部分，L保存在A的下三角部分，排列矩阵P由Ip(n)最后记录可知.(这时有|lir|1)

(3)交换A的r行与ir行元素

(4)计算U的第r行元素，L的第r列元素2023/6/6102第一百零二页，共一百五十五页，编辑于2023年，星期二

求解方程组Ly=Pb及Ux=y.

2.对于i=1,2,,n-1

(1)t←Ip(i)

(2)如果i=t则转(3)

bi

←→bt

(3)(继续循环)

3.

4.2023/6/6103第一百零三页，共一百五十五页，编辑于2023年，星期二

利用算法2的结果(实现PA=LU三角分解)，则可以计算A的逆矩阵A-1=U-1L-1P.

利用PA的三角分解计算A-1步骤：

(1)计算上三角矩阵的逆阵U-1；

(2)计算U-1L-1；

(3)交换U-1L-1列(利用Ip(n)最后记录).

上述方法求A-1大约需要n3次乘除运算.2023/6/6104第一百零四页，共一百五十五页，编辑于2023年，星期二6.3平方根法实际问题中Ax=b，系数矩阵A大多是对称正定矩阵，所谓平方根法，就是利用对称正定矩阵的三角分解而得到求解对称正定方程组的一种有效方法，目前在计算机上广泛应用平方根法解此类方程组.

对称矩阵的LDLT

分解法当A为对称矩阵,且其顺序主子式均不为0时,由定理7可知A有唯一的LU分解为(3.1)式.为了利用A的对称性，将U再分解为2023/6/6105第一百零五页，共一百五十五页，编辑于2023年，星期二其中D为对角矩阵，U0为上三角矩阵，于是

A=LU=LDU0.(3.6)又A=AT=U0TDLT,由分解的唯一性即得U0T=L.代入(3.6)得到对称矩阵A的分解式A=LDLT.总结上述讨论有下面定理.2023/6/6106第一百零六页，共一百五十五页，编辑于2023年，星期二

定理9(对称矩阵的三角分解定理)设A为n阶对称矩阵，且A的各阶顺序主子式均不为零，则A可以唯一分解为A=LDLT,其中L是单位下三角矩阵，D是对角矩阵.如果A为对称正定矩阵，则A的分解式A=LDLT中D的对角元素di均为正数.事实上,由A的对称正定性,定理6的推论成立,即2023/6/6107第一百零七页，共一百五十五页，编辑于2023年，星期二于是有由定理6得到其中L1=LD1/2为下三角矩阵.2023/6/6108第一百零八页，共一百五十五页，编辑于2023年，星期二

定理10(对称正定矩阵的三角分解或乔雷斯基(Cholesky)分解)如果A为n阶对称正定矩阵，则存在一个实的非奇异下三角矩阵L使A=LLT，当限定L的对角元素为正时，这种分解是唯一的.下面我们用直接分解方法来确定计算L元素的递推公式，因为2023/6/6109第一百零九页，共一百五十五页，编辑于2023年，星期二其中lii>0(i=1,2,,n).由矩阵乘法及ljk=0(当j<k时),得于是得到解对称正定方程组Ax=b的平方根法计算公式对于j=1,2,,n2023/6/6110第一百一十页，共一百五十五页，编辑于2023年，星期二求解Ax=b，即求解两个三角形方程组①Ly=b，求y；②LTx=y，求x.由计算公式(3.7)知所以有2023/6/6111第一百一十一页，共一百五十五页，编辑于2023年，星期二于是上面分析说明，分析过程中元素ljk的数量级不会增长且对角元素ljj恒为正数.于是有结论不选主元素的平方根法是一个数值稳定的方法.当求出L的第j列元素时，LT的第j行元素亦算出.所以平方根约需n3/6次乘除法，大约为一般直接LU分解法计算量的一半.2023/6/6112第一百一十二页，共一百五十五页，编辑于2023年，星期二例6.11用平方根法求解对称正定方程组解首先对A进行Cholesky分解求解Ly=b，得y1=2,y2=3.5,y3=1.求解LTx=y，得x1=1,x2=1,x3=1.2023/6/6113第一百一十三页，共一百五十五页，编辑于2023年，星期二例6.12平方根法求解方程组

解:因方程组系数矩阵对称正定,设A=,即：由Ly=b解得由解得

由此例可以看出，平方根法解正定方程组的缺点是需要进行开方运算。为避免开方运算，我们改用单位三角阵作为分解阵，即把对称正定矩阵A分解成

的形式，其中

2023/6/6114第一百一十四页，共一百五十五页，编辑于2023年，星期二为对角阵，而

是单位下三角阵,这里分解公式为2023/6/6115第一百一十五页，共一百五十五页，编辑于2023年，星期二据此可逐行计算

运用这种矩阵分解方法,方程组Ax=b即可归结为求解两个上三角方程组

和其计算公式分别为

和求解方程组的上述算法称为改进的平方根法。这种方法总的计算量约为，即仅为高斯消去法计算量的一半。2023/6/6116第一百一十六页，共一百五十五页，编辑于2023年，星期二总之，用平方根法解对称正定方程组时，计算L的元素ljj需要用到开方运算.为了避免开方，采用定理10的分解式A=LDLT.即由矩阵乘法，并注意ljj=1,ljk=0(j<k)，得2023/6/6117第一百一十七页，共一百五十五页，编辑于2023年，星期二于是得到L的元素及D的对角元素公式：为了避免重复计算，我们引进tij=lijdj

.由(3.9)得到按行计算L,T(T=LD)元素的公式：d1=a11

对于i=1,2,…,n.2023/6/6118第一百一十八页，共一百五十五页，编辑于2023年，星期二求解Ly=b及DLTx=y的计算公式为计算公式(3.10),(3.11)称为改进的平方根法.对于i=1,2,…,n.2023/6/6119第一百一十九页，共一百五十五页，编辑于2023年，星期二6.4误差分析6.4.1矩阵的条件数由于A(或b)元素是测量得到的,或者是计算的结果,在第一种情况A(或b)常带有某些观测误差,在后一种情况A(或b)又包含有舍入误差.因此我们处理的实际矩阵是A+A(或b+b),下面我们来研究数据A(或b)的微小误差对解的影响.即考虑估计x-y,其中y是(A+A)y=b的解.考虑线性方程组Ax=b,其中设A为非奇异矩阵,x为方程组的精确解.2023/6/6120第一百二十页，共一百五十五页，编辑于2023年，星期二首先考察一个例子.

例8设有方程组记为Ax=b,它的精确解为x=(2,0)T

.现在考虑常数项的微小变化对方程组解的影响,即考察方程组2023/6/6121第一百二十一页，共一百五十五页，编辑于2023年，星期二也可表示为A(x+x)=b+b,其中b=(0,0.0001)T,y=x+x,x为(6.1)的解.显然方程组(5.2)的解为x+x=(1,1)T.我们看到(5.1)的常数项b的第2个分量只有1/1000的微小变化,方程组的解却变化很大.这样的方程组称为病态方程组.

定义7如果矩阵A或常数项b

的微小变化(小扰动),引起方程组Ax=b解的巨大变化,则称此方程组为“病态”方程组,其系数矩阵A称为“病态”矩阵(相对于方程组而言),否则称方程组为“良态”方程组,A称为“良态”矩阵.2023/6/6122第一百二十二页，共一百五十五页，编辑于2023年，星期二应该注意,矩阵的“病态”性质是矩阵本身的特性,下面我们希望找出刻画矩阵“病态”性质的量.设有方程组

Ax=b,(5.3)其中A为非奇异矩阵,x为(5.3)的精确解.以下我们研究方程组的系数矩阵A(或b)的微小误差(小扰动)时对解的影响.2023/6/6123第一百二十三页，共一百五十五页，编辑于2023年，星期二

(1)现设A是精确的,x为Ax=b的精确解，当方程组右端有误差b,受扰解为

x+x,则

A(x+x)=b+b,Ax=b,x=A-1b,

||x||≤||A-1||||b||.(5.4)由(5.3)有||b||≤||A||||x||.于是由(5.4)及(5.5),得2023/6/6124第一百二十四页，共一百五十五页，编辑于2023年，星期二定理21设A是非奇异矩阵,Ax=b≠0

,且

A(x+x)=b+b,则上式给出了解x的相对误差的上界,常数项b的相对误差在解中放大||A-1||||A||倍.

(2)现设b是精确的,x为Ax=b的精确解，当A有微小误差(小扰动)A,受扰解为

x+x,则(A+A)(x+x)=b,(A+A)x=-(A)x.(5.6)2023/6/6125第一百二十五页，共一百五十五页，编辑于2023年，星期二如果A不受限制的话,可能A+A奇异,而(A+A)=A(I+A-1A).由定理20知,||A-1A||<1时,(I+A-1A)-1存在.由(5.6)式x=-(I+A-1A)-1A-1(A)x.因此设||A-1||||A||<1,即得2023/6/6126第一百二十六页，共一百五十五页，编辑于2023年，星期二定理22设A是非奇异矩阵,Ax=b≠0

,且

(A+A)(x+x)=b.如果||A-1||||A||<1,则(5.7)式成立.如果A充分小,且在条件||A-1||||A||<1下,那么(5.7)式说明矩阵A的相对误差||A||/||A||在解中可能放大||A-1||||A||倍.总之,量||A-1||||A||越小,由A(或b)的相对误差引起的解的相对误差就越小;量||A-1||||A||越大,解的相对误差就可能越大.所以量||A-1||||A||事实上刻画了解对原始数据变化的灵敏程度,即刻画了方程组的“病态”程度,于是引进下述定义:2023/6/6127第一百二十七页，共一百五十五页，编辑于2023年，星期二

(3)现设x为Ax=b的精确解，当A有微小误差(小扰动)A,而b同时也有微小误差b(小扰动)时,受扰解为

x+x,则还可以推出相对误差估计式为

补充2023/6/6128第一百二十八页，共一百五十五页，编辑于2023年，星期二定义8设A是非奇异矩阵,称数cond(A)v=||A-1||v||A||v(v=1,2或)为矩阵A的条件数

.由此看出矩阵的条件数与范数有关.矩阵的条件数是一个十分重要的概念.由上面讨论知,当A的条件数相对的大,即cond(A)>>1时,则(5.3)是“病态”的(即A是“病态”矩阵,或者说A是坏条件的,相对于方程组),当A的条件数相对的小,则(5.3)是“良态”的(或者说A是好条件的).注意,方程组病态性质是方程组本身的特性.A的条件数越大,方程组的病态程度越严重,也就越难用一般的计算方法求得比较准确的解.2023/6/6129第一百二十九页，共一百五十五页，编辑于2023年，星期二例如对前面例8的矩阵作分析由于条件数cond(A)1很大,可见矩阵A的病态程度十分严重,故由此方程组的解误差非常大.计算其条件数2023/6/6130第一百三十页，共一百五十五页，编辑于2023年，星期二通常使用的条件数,有(1)cond(A)∞

=||A-1||||A||;(2)A的谱条件数;当A为对称矩阵时其中1,n为A的绝对值最大和绝对值最小的特征值.2023/6/6131第一百三十一页，共一百五十五页，编辑于2023年，星期二