《计算方法》课件1第2章

第二章线性代数方程组求解方法2.1向量与矩阵基本知识2.2高斯消去法2.3矩阵的三角分解2.4矩阵的条件数与方程组的性态2.5解线性代数方程组的迭代法2.6基本迭代法2.7迭代法的收敛性习题2

2.1向量与矩阵基本知识

2.1.1引言

在自然科学和工程技术中,很多问题的解决常常归结为求解线性代数方程组:例如电学中的网络问题、船体数学放样中建立三次样条函数问题、用最小二乘法求实验数据的曲线拟合问题、解非线性方程组以及用差分法或有限元法求解偏微分方程等。这些线性代数方程组的系数矩阵大致分为两种：一种是低阶稠密矩阵；另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多)。

线性代数方程组的解法从理论上可分为两类：直接法和迭代法。

1.直接法

直接法就是不考虑计算过程中的舍入误差，经过有限次的运算得到方程组精确解的方法。本章将阐述这类算法中最基本的高斯顺序消去法及其某些变形，这类方法是解低阶稠密线性代数方程组及某些大型稀疏线性代数方程组(例如大型带状线性代数方程组)的有效方法。

2.迭代法

迭代法就是采用某种极限过程，用线性代数方程组的近似解逐步逼近精确解的方法。迭代法具有需要计算机的存储单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛性及收敛速度问题。迭代法是求解大型稀疏矩阵对应的线性代数方程组(尤其是由微分方程离散化后得到的大型稀疏线性代数方程组)的重要方法。应当注意的是，这两种方法的使用并没有严格的界限，有时在解某一个问题时会将两种方法混合使用。用直接法求解线性代数方程组时，由于在模型建立过程及求解过程中存在误差的影响，得到的解也不精确，因此在某些精度要求较高的问题中，经常把用直接法得到的解再用迭代法进行若干步迭代，以达到更高精度的要求。有关解线性代数方程组的迭代方法，我们将在第三章介绍。2.1.2向量和矩阵

下面给出一些本书要用到的定义和基本运算。

用Rm×n表示全体m×n阶实矩阵在实数域上所构成的向量空间，Cm×n表示全体m×n阶复矩阵在复数域上所构成的向量空间。

其中,ai为A的第i(i=1,2,…，n)列，bTj为A

的第j(j=1，2，

…，m)行。

(1)矩阵加法：

(2)矩阵与标量的乘法：

(3)矩阵与矩阵乘法:

(4)矩阵的转置:

(5)单位矩阵:

其中,

(6)非奇异矩阵:设A,B∈Rn×n,若AB=BA=I，则称B是A的逆矩阵，记为A-1。若A-1存在，则称A为可逆矩阵，也称A为非奇异矩阵。

注意，(A-1)T=(AT)-1，且如果AB均为非奇异矩阵，则有(AB)-1=B-1A-1。

(7)方阵的行列式：设A∈Rn×n，则A的行列式可按任一行(列)展开，即

(i=1，2，…，n)

其中，Aij为矩阵元素aij的代数余子式，Aij=(-1)i+jMij，而Mij为aij

的余子式。矩阵行列式的性质：2.1.3特殊矩阵

设A=(aij)∈Rn×n，i，j=1，2，…，n，有下列特殊矩阵：

(1)若i≠j时，有aij=0，则A为对角矩阵。

(2)若当|i-j|>1时，有aij=0，则A为三对角矩阵。

(3)若当j>i(i>j)时，有aij=0，则A为上(下)三角矩阵。

(4)若当i>j+1时，有aij=0，则A为上海森伯格(Hessenberg)矩阵。

(5)若AT=A，则A为对称矩阵。

(6)设A∈Cn×n，若AH=A(即A的元素取共轭后再转置等于A)，则A为埃尔米(Hermit)矩阵或H－矩阵。

(7)若①AT=A；②对任意非零向量x∈Rn，有xTAx>0，则A为实对称正定矩阵。

(8)若A-1=AT，则A为正交矩阵。

(9)设A∈Cn×n，若A-1=AH，则A为酉矩阵。

(10)对单位矩阵I交换其第i行与第j(j≥i)行(或交换第i列与第j(j≥i)列)，得到的矩阵记为Iij，它称为初等置换矩阵，显然:

①IijA=A1,A1为交换A中第i行与第j行得到的矩阵；

②AIij=A2,A2为交换A中第i列与第j列得到的矩阵。

(11)由有限个初等置换矩阵的乘积得到的矩阵称为置换矩阵。

定理2.1设A∈Rn×n，则下述命题等价：

(1)对任何b∈Rn×n，方程组Ax=b若有解，则只有唯一解；

(2)齐次线性方程组Ax=0只有唯一零解x=0；

(3)det(A)≠0；

(4)A-1存在；

(5)A为满秩，即rank(A)=n。

定理2.2设A∈Rn×n为实对称正定矩阵，则:

(1)A为非奇异矩阵，且A-1亦是对称正定矩阵；

(2)记Ak为A的顺序主子矩阵，则Ak也是对称正定矩阵，其中，

(3)A的特征值λi>0(i=1,2,…，n)；

(4)A的顺序主子式都大于零，即det(Ak)>0(k=1,2,…，n)。

定理2.3设A∈Rn×n为对称矩阵，若det(Ak)>0(k=1,2,…，n)，或A的特征值λi>0(i=1,2,…，n)，则A为对称正定矩阵。

定理2.4(Jordan(若当)标准型)设A为n

阶矩阵，若存在一个非奇异矩阵P，使得其中,λ1，λ2,…，λr为A的互不相同的特征值，为若当块，ni≥1(i=1,2,…，r)，且，这就是矩阵A的若当标准型。

(1)当A的若当标准型中所有若当块Ji均为一阶块时，此标准型变为对角型矩阵；

(2)若A的特征值各不相同，则若当标准型必为对角阵diag(λ1，λ2，…，λn)。2.1.4向量与矩阵的范数

首先给出向量范数定义。

定义2.1对任意向量x∈Rn，若对应一个非负实值函数‖x‖，满足下列3条时，则称‖x‖为向量x的范数：

(1)非负性：‖x‖≥0，‖x‖=0当且仅当x=0；

(2)齐次性：对任意实数α，有‖αx‖=|α|‖x‖；

(3)三角不等式性：对任意x,y∈Rn，有

‖x+y‖≤‖x‖+‖y‖。常用的范数有：

此处0<p<∞，显然,当p=1,2及p→+∞时，由p－范数可得到上面的1－范数、2－范数和∞－范数。

例2.1计算向量x=(1，-2，3)T的各种范数。

解

向量范数有如下一些性质。

定理2.5(连续性定理)设‖x‖是x∈Rn的某种向量范数，则‖x‖是分量x1，x2，…，xn的连续函数。

证明令εi=(0，…，0，1，0，…，0)T，它是Rn空间中第i(i=1，2，…，n)个单位向量，则对任何x=(x1，x2，

…，xn)T和y=(y1，y2，…，yn)T，有：根据向量范数的齐次性和三角不等式性质，有

所以当|xi-yi|→0(i=1，2，…，n)时，有‖x‖→‖y‖。

定理2.6(等价性定理)设‖·‖p以及‖·‖q是Rn上两种向量范数，则存在正常数c1，c2使得

c1‖x‖p≤‖x‖q≤c2‖x‖p

对任何x∈Rn成立。

证明当x=0时结论显然成立。下证x≠0时结论也成立。记

它是Rn中在范数‖·‖∞意义下的单位球面。由定理2.5可知，范数‖·‖p是有界闭集Ω上的连续函数，故它在Ω上有最大值M1和最小值m1。

而在Ω内，于是有：

故类似地，‖x‖q也是有界闭集Ω上的连续函数，它在Ω上也有最大值M2和最小值m2，即有：

因此或

令，，则范数的等价性说明，一种范数可由另一种范数所控制，因而一般地有，Rn上所有范数是等价的。

下面给出矩阵范数定义。

设Rn×n是n阶实矩阵集合按照实数域上矩阵的线性运算构成的线性空间。

定义2.2对任意A∈Rn×n，若对应一个非负实数函数

‖A‖，满足以下4条，则称‖A‖为矩阵A的范数：

(1)非负性，即‖A‖≥0，‖A‖=0当且仅当A=0；

(2)齐次性，即对任意实数α∈R，有

‖αA‖=|α|‖A‖；

(3)三角不等式性，即对任意A，B∈Rn×n，有

‖A+B‖≤‖A‖+‖B‖；

(4)对任意A，B∈Rn×n，有‖AB‖≤‖A‖‖B‖。进一步，若对给定的矩阵范数‖·‖M，它与某个向量范数‖·‖V满足条件(5)，则称矩阵范数‖·‖M与向量范数‖·‖V相容。

(5)对任意A∈Rn×n，x∈Rn，有

‖Ax‖V≤‖A‖M‖x‖V成立。

设A=(aij)n×n∈Rn×n，常用的矩阵范数有：

在矩阵范数中还有一种由向量范数导出的矩阵范数。

定义2.3设x∈Rn，A∈Rn×n，‖x‖p为向量x的某种向量范数，定义矩阵A的非负函数为

可以验证它满足矩阵范数的条件，称它为从属于向量范数‖·‖p的矩阵范数，简称为从属范数，

有时也称为算子范数。

显然，由向量范数‖x‖p所导出的矩阵范数‖A‖p与该向量范数‖x‖p是相容的。

关于从属范数，有下述结论。

定理2.7矩阵范数‖A‖1，‖A‖∞，‖A‖2分别是向量范数‖x‖1，‖x‖∞，‖x‖2的从属范数。

证明

(1)先证明是‖x‖∞的

从属范数。

设x=(x1,x2,…，xn)T≠0，A≠0,令则于是对任意非零向量x∈Rn，有

下面说明至少有一个x0∈Rn,x0≠0使设

取向量x0=(β1,β2,…，βn)T，

(j=1，2，

…，n)，由符号函数sign定义知‖x0‖∞=1，而

即

(2)再证明是‖x‖1的从属范数。

设则有：

所以对任意非零向量x∈Rn，都有:设

取向量x0=(β1,β2,…，βn)T，其中

,

则显然有‖x0‖1=1,而

所以

(3)最后证明是‖x‖2的从属范数。

显然对任意x∈Rn，有xTATAx≥0，从而ATA为对称半正定矩阵。设ATA的特征值从大到小依次为λ1≥λ2≥…≥

λn≥0，特征值相应的标准正交特征向量为u1,u2,…，un，

对任意非零向量x∈Rn，有:

所以另一方面，若取x=u1，则有:

故2.2高斯消去法

高斯顺序消去法是一个古老的求解线性代数方程组的方法(早在公元前250年，我国就掌握了解线性代数方程组的消去法)，对这个方法进行改进，变形得到的主元消去法、三角分解法等仍是目前计算机上常用的方法。2.2.1高斯顺序消去法

设有线性代数方程组

(2.1)或写为矩阵形式

简记为Ax=b。首先举一个简单的例子来说明顺序消去法的基本思想。

例2.2用顺序消去法解线性代数方程组

解第一步，将此方程组第一个方程乘以－2，加到方程组的第三个方程上，消去第三个方程中的未知数x1，得到与原方程组等价的方程组第二步，将此方程组第二个方程加到方程组第三个方程上，消去第三个方程中的未知数x2，得到与原方程组等价的三角形方程组

对此方程组经过回代，容易求得其解为x*=(1，2，3)T。上述过程用矩阵表述如下:由此看出，高斯顺序消去法的基本思想是：对线性代数方程组(2.1)式所对应的增广矩阵(A

b)进行一系列“把某一行的非零常数倍加到另一行上”的初等行变换，使得

(A

b)中A的对角线以下的元素全变为0，从而使(2.1)式等价地转化为容易求解的上三角形线性代数方程组，再通过回代得到上三角形线性代数方程组的解，即得到线性代数方程组(2.1)式的解。为方便叙述，我们记

(2.2)

即

1.消去过程

第一步，设，令,把(2.2)式中第1行的li1倍依次加到(2.2)式中的第i行(i=2，3，…，n)，则(2.2)式中第i(i≥2)行第j列(2≤j≤n)位置处的元素为因此，(2.2)式化为

(2.3)第二步，设，令，把(2.3)式中第2行的li2倍依次加到(2.3)式中的第i行(i=3，4，…，n)，则(2.3)式中第i(i≥3)行第j列(3≤j≤n)位置处的元素为因此，(2.3)式化为

(2.4)依次做下去，一直做到第n－1步，即有

(2.5)

至此,A(1)x=b(1)转化为A(n)x=b(n)。

2.回代过程

容易看出，线性代数方程组(2.1)式经过上述n-1步消去过程以后变为(2.5)式所对应的上三角形线性代数方程组。若

以及上述消去过程中假设,

…,

，则可得xn,xn-1,…，x2，x1为

(k=n－1，n－2，…，2，1)

(2.6)

此即为(2.1)式的解。注意：在Ax=b中，由于A∈Rn×n，为非奇异矩阵，如果a11=0，则A的第一列一定有元素不等于0，比如ai1≠0，于是交换第一行和第i行，将ai1调到(1，1)的位置，然后进行消元计算，这时A(2)右下角矩阵为n－1阶非奇异矩阵。继续此过程，高斯消去法照样可继续进行。总结上述讨论即有定理2.8。

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

(1)如果(k=1,2,…，n)，则可通过高斯顺

序消去法将Ax=b化为等价的上三角形线性代数方程组(2.5)式，且计算公式为

(a)消元计算(k=1，2，…，n－1)：

(b)回代计算：

(2)如果A为非奇异矩阵，则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为(2.5)式。

(3)如果A的各阶顺序主子式di≠0(i=1,2,…，n)，则可通过高斯顺序消去法将方程组Ax=b约化为(2.5)式，即

(2.7)下面计算高斯顺序消去法中乘除法运算的工作量。

(1)消去过程中的计算量。在第k步的消去过程中，对矩阵A(k)的元素作n－k次除法运算和(n－k)2次乘法运算，对b(k)作n－k次乘法运算，于是消去过程中总的乘除运算量为

(2)回代过程中的计算量。计算xk需作(n－k)次乘法运算和一次除法运算，于是回代过程中的运算量为

因此，用高斯顺序消去法计算线性代数方程组(2.1)式解的乘除法运算次数为在浮点机上完成一次乘除运算比完成一次加减运算所耗的时间要多得多，故当一个算法中的加减运算次数与乘除运算次数相差不多时，仅考虑乘除运算次数就可以体现该算法的计算量。

如果在消去的过程中考虑进行行交换，那么求线性代数方程组(2.1)式的解所进行的乘除法运算量不超过对于一个含有20个未知量的线性代数方程组，用高斯消去法求解，总的乘除法运算量不超过3440次，用每秒1亿次运算速度的计算机来完成则是一瞬间的事情，而如果采用线性代数课程中学过的克莱姆法则求解，则其乘除法的运算量大约为9.7×1020，也就是说，即使用每秒1亿次运算速度的计算机计算，也需要30多万年。因此，算法的好坏对计算能力的提高至关重要。

例2.3采用四位有效数字，求下列线性代数方程组的解：

解方法1采用不进行行交换的高斯消去法求解。把的第一行乘以－1764加到第二行,有:

回代可得x2=1.001，x1=－10.00。

方法2采用进行行交换的高斯消去法求解。

交换增广矩阵第一、二行，有：把的第一行乘以－0.0005670加到第二行,有:

回代可得x2=1.000，x1=10.00。从原线性代数方程组可看出精确解为x1=10.00，x2=1.000，即为方法2所求的解。那么方法1得出的解为何会有如此大的误差呢？这主要是因为在计算x2时，产生的小误差0.001导致了在求解x1时产生了更大的误差，即误差扩大了20000倍。

从以上例子可以看出，高斯顺序消去法是有缺陷的。它并没有考虑同列对角线以下元素绝对值的大小关系，有时会导致用绝对值较小的数字做分母而出现较大的数字与较小数字进行相加减运算而被“吞掉”的现象，以至于产生某一未知元较小的误差却引起其它元更大的误差。为了改变这一缺陷，实际应用时较多采用高斯主元消去法。2.2.2高斯主元消去法

根据主元选取范围不同，主元消去法又分为列主元消去法和全主元消去法。

1.列主元消去法

设线性代数方程组(2.1)式的增广矩阵为首先在的第一列中选取绝对值最大的元素

作为第一列的主元，即

然后交换中的第一行与第i行，经一次消元计算得

此消去过程与高斯顺序消去法完全相同。重复上述过程，设已完成第k－1步的选主元素，交换两行及消元过程后(A

b)已约化为第k步选主元素，在右下角方阵的第一列内选取绝对值最大的元素作为这一列的主元，即

然后交换的第i行与第k行，再进行消元计算。如此重复上述过程，直到最后将原线性代数方程组约化为回代求解

从上面分析可以看出，列主元消去法除了每步需要按列选出主元，然后进行对换外，其消去过程与高斯顺序消去法是相同的。

2.全主元消去法

全主元消去法是对于一般的矩阵，每一步都在要处理的矩阵(或消元后的低阶矩阵)中选取绝对值最大的元素作为主元，从而使高斯消去法具有更好的数值稳定性的方法。虽然全主元消去法的求解结果更加可靠，但由于全选主元每步耗时更多，而且要进行列交换，那么所求未知量x1,x2,…,xn的次序就会被打乱，因此，实际应用中一般使用列主元消去法。

2.3矩阵的三角分解

由线性代数知道，对矩阵进行一次初等变换，就相当于用相应的初等矩阵去左乘或右乘原来的矩阵。因此，我们把上述求解线性代数方程组的高斯消去法用矩阵乘法表示出来，即可得到求解线性代数方程组的另一种直接方法——矩阵的三角分解法。设(2.1)式的系数矩阵A∈Rn×n的各阶顺序主子式均不为0。由于对A施行初等行变换相当于用相应的初等矩阵左乘矩阵A，因此对(2.1)式施行第一次消元后，化为

，其中A(1)化为A(2)

，b(1)化为b(2)，即其中，一般第k步消元后，化为，其中A(k)化为A(k+1)，b(k)化为b(k+1)，相当于

LkA(k)=A(k+1)，Lkb(k)=b(k+1)

其中，重复这个过程，最后得到

(2.8)

其中,将上三角形矩阵A(n)记为U，由(2.8)式得

其中,

为单位下三角矩阵。由以上讨论可知，高斯顺序消去法实质上产生了一个将A分解为两个三角形矩阵相乘的因式。为明确起见，对矩阵的三角分解有如下定义。

定义2.4设A为n(n≥2)阶方阵，称A=LU为矩阵A的三角分解，其中L是下三角矩阵，U是上三角矩阵。

定义2.5若L是单位下三角矩阵，U是上三角矩阵，则三角分解A=LU称为杜利脱尔(Doolittle)分解；若L是下三角矩阵，U是单位上三角矩阵，则称A=LU为克罗脱(Grout)分解。应该注意的是，矩阵A=LU的三角分解是不唯一的，如

其中，L1=LD，U1=D－1U。对角矩阵D为可逆的非单位矩阵，显然L1≠L，U1≠U。为了保证A的三角分解的唯一性，我们有如下定理。

定理2.9如果非奇异矩阵A满足下列三个条件之一，

则矩阵A有唯一的杜利脱尔分解A=LU和唯一的克罗脱分解

。

(1)A的各阶顺序主子式det(Ak)>0(k=1,2,…，n)；

(2)

(3)A为对称正定矩阵。

证明仅给出条件(1)下杜利脱尔分解唯一性的证明，其余情形读者可练习证明。

根据以上高斯消去法的矩阵分析，当det(Ak)>0(k=1,2,…,n)时，杜利脱尔分解过程可以进行到底，即A=LU的存在性已证。以下证明分解的唯一性。

设A为非奇异矩阵，且

A=LU=L1U1

其中,L、L1

为单位下三角矩阵，U、U1为上三角矩阵。由于U１－1存在，因此

容易证明上(下)三角矩阵的逆矩阵仍然为上(下)三角矩阵，因而上式右边为上三角矩阵，左边为单位下三角矩阵，故上式要成立，两边都必须等于单位矩阵，从而U=U1，L=L1。证毕。2.3.1直接三角分解法

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

Ly=b，求y

Ux=y，求x

的问题，而这两个线性代数方程组只要回代，就可以求出其解。

1.不选主元的三角分解法

设A为非奇异矩阵，且有分解式A=LU，其中L为单位下三角矩阵，U为上三角矩阵，即

(2.9)第一步，用L的第一行分别乘U的第j(j=1，2，…，n)列，比较两边可得

a1j=u1j

(j=1，2，…，n)

分别用L的第i(i=2，3，…，n)行乘U的第一列，比较两边可得

ai1=li1u11

(i=2，3，…，n)

即得

(i=2，3，…，n)

这样就求出了L的第一列和U的第一行的所有元素。第二步，用L的第二行分别乘U的第j(j=2，3，…，n)列，比较两边可得

a2j=l21u1j+u2j

(j=2，3，…，n)

即得

u2j=a2j－l21u1j

(j=2，3，…，n)

分别用L的第i(i=3，4，…，n)行乘U的第二列，比较两边可得

ai2=li1u12+li2u22

即得

这样就求出了L的第二列和U的第二行的所有元素。依次进行下去，一直到第k－1步，即已求出L的前k－1列和U的前k－1行的所有元素。

第k步，用L的第k行分别乘U的第j(j=k，k+1，…，n)列，比较两边可得

akj=lk1u1j+…+lk，k－1uk－1，j+ukj

(j=k，k+1，…，n)

即得

ukj=akj－(lk1u1j+…+lk，k－1uk－1，j)

(j=k，k+1，…，n)分别用L的第i(i=k+1，…，n)行乘U的第k列，比较两边可得

aik=li1uik+…+li，k－1uk－1，k+likukk

(i=k+1，…，n)

即得总结上述讨论，得到用直接三角分解法求解Ax=b(要求A的所有顺序主子式都不为零)的计算公式：及求解Ly=b，Ux=y的计算公式：

例2.4用杜利脱尔分解法求线性代数方程组的解：

解设第一步有：

u11=1,

u12=2,

u13=3

l21=1,

l31=1

第二步有：

u22=1,

u23=2

l32=1

第三步有：

u33=1

l33=1于是有

由Ly=b得y=(2,1,1)T；由Ux=y得x=(1，－1，1)T。

直接三角分解法大约需要次乘除法运算，与高斯顺序消去法的计算量基本相同。

2.选主元的三角分解法

从直接三角分解公式可以看出，当ukk=0时，或者当ukk绝对值很小时，按分解公式计算可能会引起舍入误差的积累，但如果A为非奇异矩阵，则可通过对A进行行交换后再分解的方法实现LU分解。也就是说，可采用与列主元消去法类似的方法，将直接三角分解法修改为(部分)选主元的三角分解法。

定理2.10对非奇异矩阵A，存在n阶置换阵

(ik≥k)，使得PA有杜利脱尔分解，即PA=LU。

证明由于A为非奇异矩阵，因而有

Ln－1Pn－1…L2P2L1P1A=U

成立，此式等价于用带行交换的高斯消元法将A化为上三角形矩阵。其中Lk为单位下三角矩阵；Pk为初等置换矩阵，表示进行第k行与它下面的某行进行交换，且PkPk=I；U

为上三角形矩阵。从而有其中，容易验证为单位下三角矩阵，因此

例2.5用带行交换的杜利脱尔分解求线性代数方程组的解。

解

于是于是

即令P=P2，L－1=P2L1P2，则

这样我们得到A的杜利脱尔分解为PA=LU。因此，由Ax=b有PAx=Pb，即LUx=Pb，故

回代计算可得y=(5，7，4)T，x=(2，3，4)T。2.3.2平方根法

许多问题的求解所归纳出来的线性代数方程组Ax=b，其系数矩阵常常具有对称正定性。由上节讨论可知，这样的系数矩阵A一定有唯一的杜利脱尔分解，因而完全可用前述方法分解并求解。然而利用矩阵A的对称正定性，我们可以建立更好的三角分解法，这就是下面的平分根法。

定义2.6设A是n(n≥2)阶实对称正定矩阵，L是非奇异的下三角矩阵，则称A=LLT为矩阵A的乔列斯基(Cholesky)分解。

乔列斯基分解法计算公式如下所述。设A为实对称正定矩阵，令A=LLT，即第一步，分别用L的第i(i=1，2，…，n)行乘LT的第一列，比较两边元素可得

a11=l11l11，ai1=l11li1

(i=2，3，…，n)

即得

第二步，分别用L的第i(i=2，3，…，n)行乘以LT第二列，比较两边元素可得

即得依次进行下去，直到第n－1步。

第n步，用L的第n行乘LT的第n列，比较两边元素可得

即得综上所述，可以得到解对称正定方程组Ax=b的平方根法计算公式，对于i=1，2，…，n，有：而求解线性代数方程组Ax=b相当于求解两个三角形方程组：

Ly=b

(求y)

LTx=y

(求x)

从而有：由计算公式①知：

所以这说明在乔列斯基分解过程中，元素lis的平方不会超过A的最大对角元素，因而舍入误差的放大受到了限制。因此，用平方根方法求解对称正定矩阵所对应的线性代数方程组时可以不考虑选主元问题。平方根法约需次乘除法，大约为一般直接LU分解法计算量的一半。

在上述计算lis的公式中，由于每步都有开方运算，因此也称乔列斯基分解法为平方根法。在用平方根法解对称正定矩阵所对应的线性代数方程组时，不可避免地要用到开方运

算，为了避免开方运算，我们可以对乔列斯基分解法作以改进，从而得到以下定理。

定理2.11

n(n≥2)阶对称正定矩阵A必有如下分解：

A=LDLT

其中，L是单位下三角矩阵；D是对角元素全为正的对角矩阵，并且这种分解是唯一的。

证明因为A为对称正定矩阵，A的顺序主子式全大于0，从而A有唯一的杜利脱尔分解，因为

所以故有将等式右端三个矩阵分别记为LDR，因为AT=A，所以有：

A=LDR=RTDLT=AT

或

A=L(DR)=RT(DLT)

由杜利脱尔分解的唯一性可知，L=RT，LT=R，从而

A=LDLT

即上式分解的唯一性成立。应用定理2.11的分解式A=LDLT解线性代数方程组Ax=b，其分解计算量与LLT分解计算量差不多，但LDLT分解不需要开方计算。

例2.6用改进的平方根法解方程组：

解

设由矩阵乘法得即得于是有

从而由Ly=b得y=(2,3,1)T；由DLTx=y得x=(1,1,1)T。

关于乔列斯基分解的唯一性有如下定理。

定理2.12

n(n≥2)阶对称正定矩阵A一定有乔列斯基分解A=LLT，当限定L对角线上元素全为正时，乔列斯基分解是唯一的。

证明由定理2.11可知，存在单位下三角矩阵L1，对角线元素全为正的对角矩阵D=diag(u11,u22，…，unn)，使得

记

则有

再令，显然L是非奇异的下三角矩阵，故A有乔列斯基分解A=LLT，当L的对角线元素全为正时，由杜利脱尔分解的唯一性可推知乔列斯基分解的唯一性。2.3.3解三对角方程组的追赶法

有些实际问题，如二阶常微分方程边值问题、热传导问题等导出的线性代数方程组的系数矩阵是三对角矩阵，即Ax=d，其中，并且满足条件

即三对角矩阵A是严格对角占优矩阵，此时A的各阶顺序主子式不为零，因此A有唯一的杜利脱尔分解。但根据A是三对角矩阵的特点，易知三对角矩阵有如下更特殊的杜利脱尔三角分解式：

A=LU且由矩阵乘法运算，比较上式两边对应元素可得如下计算公式：从而将求解方程组Ax=d化为求解两个二对角方程组Ly=d和Ux=y，计算公式如下：

如果把A进行克罗脱分解，则有：

即有由矩阵乘法运算，比较上式两边对应元素可得如下计算公式：从而将求解方程组Ax=d化为求解两个方程组和，计算公式为追赶法公式实际上就是把高斯消去法用到求解三对角线方程组上去的结果。这时由于A特别简单，因此使得求解的计算公式非常简单，而且计算量仅为5n－4次乘法，而另外增加一个方程组仅增加3n－2次乘法运算，可见追赶法的计算量是比较小的。

例2.7用追赶法求解三对角线性代数方程组：

解令A=LU，即利用矩阵乘法运算，再比较两边对应元素可得从而从而由Ly=d及Ux=y可算出

2.4矩阵的条件数与方程组的性态

考虑线性代数方程组Ax=b，其中A为非奇异矩阵，x为方程组的精确解。由于A(或b)元素是测量得到的，或者是计算的结果，在第一种情况下，A(或b)常带有某种观测误差，而在后一种情况下A(或b)又包含有舍入误差，因此参与线性代数方程组运算的实际矩阵是A+δA(或b+δb)。本节讨论线性代数方程组的解对系数矩阵A与常数项b的扰动(误差)的敏感性问题。先来观察两个例子。

例2.8线性代数方程组：

的精确解为(1,1)T，若A及b作微小扰动，扰动后的线性代数方程组为则其精确解为(－2，7)T。由此可见，在本例中，A与b的微小变化引起了x的很大的变化，这表明x对A与b的扰动是敏感的。此例中，detA=－10－4，方程组由两条几乎平行的直线组成，因而当其中一条直线有很小的变化时，新交点与原交点相差较远。

例2.9线性代数方程组：的精确解为x=(9.2,－12.6,4.5,－1.1)T，如果对方程组的系数矩阵A作微小扰动，则变为

从而其精确解为x=(－8.1,137,－34,22)T。同样，解也有很大的误差。下面来详细分析线性代数方程组Ax=b中向量b和系数矩阵A的微小扰动对解x的影响。首先来看一些相关基本概念。

如例2.8和例2.9所示，当方程组的系数矩阵A或b有微小变化时，方程组的解可能会产生很大的变化，这样的方程组称为病态方程组。

定义2.7如果矩阵A或常数项b的微小变化，引起方程组Ax=b解的巨大变化，则称此方程组为“病态方程组”，相对于方程组而言，矩阵A称为“病态”矩阵，否则称方程组

为“良态”方程组，A为“良态”矩阵。

应该注意，矩阵的“病态”性质是矩阵本身的特性，下面希望找出刻画矩阵“病态”性质的量。设有方程组：

Ax=b

(2.10)其中，A为非奇异矩阵，x为(2.10)式的精确解，以下我们来研究方程组的系数矩阵A(或b)的微小误差(扰动)对解的影响问题。

1.常数项b有扰动

设A是精确的，b有误差δb，解为x+δx，即

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

则

δx=A－1δb

‖δx‖≤‖A－1‖‖δb‖

(2.11)

由(2.10)式得‖b‖≤‖A‖‖x‖，即

(2.12)

于是由(2.11)式及(2.12)式得定理2.13。

定理2.13设A是非奇异矩阵，Ax=b≠0,且A(x+δx)=b+δb,则

(2.13)

上式给出了解的相对误差的上界，常数项b的相对误差在解中可能放大‖A－1‖‖A‖倍。

2.系数矩阵A有扰动

设(2.10)式中b是精确的，A有微小误差(扰动)δA，解为x+δx，即

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

则

(A+δA)δx=－(δA)x

而

A+δA=A(I+A－1δA)当‖A－1δA‖<1时，(I+A－1δA)－1存在，则有

δx=－(I+A－1δA)－1A－1(δA)x

因此设

‖A－1‖‖δA‖<1

即得

(2.14)

3.系数矩阵A和常数项b都有扰动

设(2.10)式中A和b都有小扰动时，解为x+δx，即

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

则

(A+δA)δx=δb－δAx

得

δx+A－1δAδx=A－1δb－A－1δAx从而

‖δx‖－‖A－1‖‖δA‖‖δx‖≤‖A－1‖‖δb‖

+‖A－1‖‖δA‖‖x‖

两边同除‖x‖，并由(2.12)式得当‖A－1‖‖δA‖<1时，可得

(2.15)

分析(2.13)~(2.15)式可知，无论方程组(2.10)式中右端常数项b有扰动，还是系数矩阵A有扰动，或者两者都有扰动，总之，量‖A－1‖‖A‖愈小，由A(或b)的相对误差

引起的解的相对误差就愈小，量‖A－1‖‖A‖愈大，解的相对误差就愈大，所以量‖A－1‖‖A‖实际上刻画了解对原始数据变化的灵敏程度，即刻画了方程组的“病态”程度。

定义2.8设A为n阶非奇异矩阵，称数

Cond(A)v=‖A‖v‖A－1‖v

为矩阵A的条件数。其中，‖·‖v为Rn×n中的某种矩阵范数。

由此看出，矩阵的条件数与矩阵范数有关。矩阵的条件数是一个十分重要的概念，由上面讨论知，当A的条件数相对比较大，即Cond(A)>>1时，方程组Ax=b是“病态”的(即A为“病态”或者说相对于解方程组来说A是坏条件的)；当A的条件数相对较小时，方程组Ax=b是“良态“的(或者说A是好条件的)。注意，方程组病态性质是方程组本身的特性，A的条件数愈大，方程组的病态程度愈严重，也就愈难用一般的计算方法求比较精确的解。常用的条件数有：条件数的性质有：

(1)对任意非奇异矩阵A∈Rn×n，有Cond(A)≥1；

(2)对任意非奇异矩阵A∈Rn×n及c∈R，c≠0，有Cond(cA)=Cond(A)；

(3)对于任意正交矩阵A，有Cond(A)2=1;

(4)如果P为正交矩阵，A为非奇异矩阵，则有Cond(PA)2=Cond(AP)2=Cond(A)2。

例2.10已知希尔伯特(Hilbert)矩阵:

计算H3的条件数。

解

(1)计算H3的条件数Cond(H3)10：

同理可计算得Cond(H6)∞=2.9×107，Cond(H7)∞=

9.85×108。可见，随着n的增大，Cond(Hn)2急剧增加，因此，以Hn为系数矩阵的线性代数方程组Hnx=b是严重病态的。

(2) 考虑方程组

设H3及b有微小误差(取3位有效数字)，有：

(2.16)简写为(H3+δH3)(x+δx)=b+δb，方程组H3x=b与(2.16)式的精确解分别为

x=(1，1，1)T

x+δx=(1.089512538,0.487967062,1.491002798)T于是δx=(0.089512538,－0.512032938,

0.491002798)T，可以算得：

这就是说,H3与b的相对误差不超过0.2%，而引起的解的相对误差则超过了50%。由上述例题可以看出，对于给定的线性代数方程组Ax=b，要判断它是否病态并不容易，因为用条件Cond(A)=‖A‖‖A－1‖判断时要求出A－1，而计算A－1是比较费劲的，那么在实际计算中如何发现病态情况呢？

(1)如果在系数矩阵A的三角分解时(尤其是用主元素消去法解(2.10)式)出现小主元，则对大多数矩阵来说，A可能是病态的。

(2)系数矩阵A的行列式值相对很小，或系数矩阵某些行近似线性相关，则A可能是病态的。

(3)系数矩阵A的元素间数量级相差很大，并且无一定规则，则A可能是病态的。一般病态线性代数方程组的求解是比较困难的，线性代数方程组给定时，其系数矩阵的条件数也是随之确定的，因而线性代数方程组的性态是方程组本身固有的性质，与求解线性方程组的算法无关。

在计算机上求解线性方程组，都是所给方程组的扰动方程组，这是因为将系数矩阵和常数项输入计算机后，机器要对数作十进制和二进制转化，由于字长的限制总有误差。对于良态方程组，只要求解方法稳定，即可得到比较满意的计算结果；但对于病态方程组，即使采用稳定性好的算法也未必能得到理想解。对于病态方程组Ax=b，实际应用中使用下述两种方法进行处理：

(1)采用高精度的算术运算。如采用双精度运算，使由于舍入误差的放大而损失若干有效数位后，还能保留一些有效数位，改善或减轻“病态”影响。

(2)对方程组进行预处理。如用可逆对角矩阵对方程组进行矩阵平衡，降低系数矩阵的条件数。具体作法是寻找可逆对角矩阵D，使方程组Ax=b转化为等价方程组DAx=Db，并且使DAx=Db的条件数相对较小，如计算(i=1，2，…)，取D=diag(S－11，S－12，

…，S－1n)，那么Cond(DA)的值比Cond(A)要小得多。

例2.11对线性代数方程组：

(2.17)

计算条件数Cond(A)∞。

解

可见原方程组为病态方程组。现对方程组进行预处理。

取，

得

设DA=A′，Db=b′,则有同解方程组:

A′x=b′

(2.18)而

于是当用列主元消去法解(2.17)式时(计算到三位数字)：

得到很坏的结果：x2=1，x1=0。

现用列主元消去法解(2.18)式，得到

从而得到较好的解：x1=1，x2=1。2.5解线性代数方程组的迭代法

考虑线性代数方程组:

Ax=b

(2.19)其中,A为非奇异矩阵，当A为低阶稠密矩阵时，选主元消去法是解(2.19)式的有效方法，但是，对由工程技术中产生的大型稀疏矩阵方程组(A的阶数n

很大，零元素很多，例如求某些偏微分方程数值解所产生的线性方程组，n≥104)，利用迭代法求解(2.19)式是合适的。在计算机的内存和运算方面，迭代法通常都可利用A中有大量零元素的特点。本节将介绍迭代法的一些基本理论及雅可比(Jacobi)迭代(简称J－迭代)、高斯赛德尔(Gauss-Seidel)迭代(简称GS-迭代)、逐次超松弛(SuccesiveOver-Relaxation)迭代(简称SOR-迭代)和对称逐次超松弛(SymmetricSuccessiveOverRelaxation)迭代(简称SSOR迭代)。其中，超松弛迭代法应用很广泛。

先看下面简例，以便了解迭代法的基本思想。

例2.12求解线性代数方程组：

(2.20)

解记为Ax=b，其中，线性代数方程组的精确解x*=(3，2，1)T。现将(2.20)式改写为

(2.21)

也可简写为

x=Bx+f其中，任取初始值，比如取x(0)=(0，0，0)T，代入(2.21)式右边(若(2.21)式为等式即求得线性代数方程组的解，但一般不满足)，得到新的值x(1)=(x(1)1，x(1)2，x(1)3)T=(2.5,3,3)T，再将x(1)代入(2.21)式右边得到x(2)。反复利用这个计算程序，得到一个向量序列和一般的计算公式(迭代公式)：

(2.22)简写为

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

其中，k表示迭代次数(k=0，1，2，…)。

迭代到第10次时有：

x(10)=(3.000032,1.999838,0.9998813)T

‖ε(10)‖∞=0.000187

(ε(10)=x(10)－x*)

从此例看出，由迭代法产生的向量序列x(k)逐次逼近线性代数方程组的精确解x*。对于任何一个线性代数方程组x=Bx+f(由Ax=b变形得到的等价方程组)，由迭代法产生的向量序列x(k)是否一定逐步逼近线性代数方程组的解x*呢？回答是不一定。请读者考虑用迭代法解下述方程组：对于给定方程组x=Bx+f，设有唯一解x*，则

x*=Bx*+f

(2.23)

又设x(0)为任取的初始向量，按下述公式构造向量序列

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

(k=0，1，2，…)(2.24)其中，k表示迭代次数。

定义2.9

(1)对于给定的线性代数方程组x=Bx+f，用(2.24)式逐步代入求近似解的方法称为迭代法(或称为一阶定常迭代法，这里B与k无关)。

(2)如果存在(记为x*)，则称此迭代法收敛，显然，x*就是线性代数方程组的解，否则称迭代法发散。由上述讨论，需要研究{x(k)}的收敛性。引入误差向量：

ε(k+1)=x(k+1)－x*

由(2.24)式减去(2.23)式得

ε(k+1)=Bε(k)

(k=0,1,2,…)

递推得

ε(k)=Bε(k－1)=…=Bkε(0)

可见要考察{x(k)}的收敛性，就要研究B在什么条件下有

，即要研究B满足什么条件时有。2.6基本迭代法

对线性代数方程组:

Ax=b(2.25)

其中，A=(aij)∈Rn×n为非奇异矩阵，下面研究如何建立解Ax=b的各种迭代方法。

将A分裂为

A=M－N

(2.26)其中，M为非奇异矩阵，且要求线性代数方程组Mx=d容易求解，一般选择为A的某一部分元素构成的矩阵，称M为A的分裂矩阵。于是，求解Ax=b转化为求解Mx=Nx+b，

即

由此可构造一个迭代法

(2.27)

其中，B=M－1N=M－1(M－A)=I－M－1A，f=M－1b，称

B=I－M－1A为迭代法的迭代矩阵。通过选取不同的M矩阵，就可得到解Ax=b的各种迭代方法。设aii≠0(i=1,2,…，n)，将A分为三部分：

(2.28)2.6.1雅克比迭代法(J-迭代法)

由aii≠0(i=1,2,…，n)，选取M为A的对角元素组成的矩阵，即选取M=D(对角阵)，因此A=D－N。由(2.27)式得到解Ax=b的雅克比迭代法：

(2.29)

其中，BJ=I－D－1A=D－1(L+U)，f=D－1b，称BJ为求解Ax=b的雅克比迭代法的迭代矩阵。以下给出雅克比迭代的分量计算公式。令x(k)=(x(k)1，x(k)2，…，x(k)n)T，由雅克比迭代公式(2.29)有：

即有于是，解Ax=b的雅克比迭代法的计算公式为

(2.30)

由(2.30)式可知，雅克比迭代法的计算公式简单，每迭代一次只需计算矩阵和向量的乘法一次，且计算过程中原始矩阵A的元素始终不变。

例2.13用J迭代法计算线性代数方程组：

的近似解x(k+1)，要求‖x(k+1)－x(k)‖∞<10－5(精确解为x1=3,x2=4,x3=－5)。

解应用J-迭代法，取x(0)=(1，1，1)T，得到迭代公式对k=0，1，…，计算得

x(1)=(5.250000,7.000000,-5.750000)T

x(2)=(0.750000,2.125000,-4.250000)T

x(58)=(2.999996,3.999996,-4.999999)T

且满足‖x(58)－x(57)‖∞<10－5。2.6.2高斯-赛德尔迭代法(GS-迭代法)

选取矩阵M为A的下三角部分，即M=D－L(下三角矩阵)，因此A=M－N，由(2.27)式得到解Ax=b的高斯赛德尔迭代法：

(2.31)

其中，BG=I－(D－L)－1A=(D－L)－1U，f=(D－L)－1b，

称BG=(D－L)－1U为解Ax=b的高斯-赛德尔迭代法的迭代矩

阵。以下给出高斯-赛德尔迭代法的分量计算公式。记x(k)=(x(k)1，x(k)2，…，x(k)n)T，由(2.31)式有：

(D－L)x(k+1)=Ux(k)+b

即

得于是，解Ax=b的高斯赛德尔迭代法的计算公式为

(2.32)雅克比迭代法不使用变量的最新信息计算xi(k+1)，而由高斯-赛德尔迭代公式(2.32)可知，计算x(k+1)的第i个分量xi(k+1)时,利用了已经计算出的最新分量xj(k+1)(j=1，2，…，i－1),高斯-赛德尔迭代法可看做是雅克比迭代法的一种改进。由(2.32)式可知，高斯-赛德尔迭代法每迭代一次只需计算一次矩阵与向量的乘法。

例2.14用高斯-赛德尔迭代法解线性方程组(2.20)。

解应用高斯-赛德尔迭代公式，取x(0)=(0，0，0)T，得对k=0，1，…，计算

x(1)=(2.500000,2.090909,1.227273)T

x(7)=(3.000002,1.9999987,0.9999932)T

此时‖x*－x(7)‖∞<2.02×10－6。由此例可知，用高斯-赛德尔迭代法、雅克比迭代法解同一线性代数方程组(2.20)(且取x(0)=(0,0,0)T)时均收敛，而高斯-赛德尔迭代法比雅克比迭代法收敛快，(即取相同初始值x(0)，在同样精度要求下，高斯-赛德尔迭代法所需迭代次数较少)，