数值分析ppt第5章-解线性方程组的直接方法

第5章解线性方程组的直接方法5.1引言与预备知识5.2高斯消去法5.3高斯主元素消去法5.4矩阵三角分解法5.5向量和矩阵的范数5.6误差分析5.7矩阵的正交三角化及应用这一章讨论线性方程组5.1引言与预备知识的数值解法.在自然科学和工程技术中,很多问题归结为解线性方程组.例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵(例如，阶数不超过150).另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多).5.1.1引言有的问题的数学模型中虽不直接表现为含线性方程组,但它的数值解法中将问题“离散化”或“线性化”为线性方程组.因此线性方程组的求解是数值分析课程中最基本的内容之一.

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

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

2.迭代法5.1.2向量和矩阵

基本概念：用Rm×n表示全部m×n实矩阵的向量空间；用Cm×n表示全部m×n复矩阵的向量空间.(由数排成的矩阵表，称为m行n列矩阵).其中aj为A的第j列的m维列向量.同理(m维列向量).其中biT为A的第i行的n维行向量.

矩阵的基本运算：

(1)矩阵加法

(2)矩阵与标量的乘法

(3)矩阵与矩阵的乘法

(4)转置矩阵

(5)单位矩阵其中

(6)非奇异矩阵则称B是A的逆矩阵，记为A-1，且(A-1)T=(AT)-1.如果A-1存在，则A称为非奇异矩阵.如果A、B均为非奇异矩阵，则有(AB)-1=B-1A-1.

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

行列式性质：5.1.3特殊矩阵设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的共轭转置矩阵).

(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)置换阵由初等置换阵的乘积得到的矩阵.

定理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.

定理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).

定理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使得为若当(Jordan)块.其中

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

(2)如果A的特征值各不相同，则其若当标准型必为对角矩阵diag(λ1,λ2,,λn).5.2高斯消去法

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

例1

用消去法解方程组

解第1步，将方程(2.2)乘上-2加到方程(2.4)上去，消去(2.4)中的未知数x1，得到第2步，将方程(2.3)加到方程(2.5)上去，消去(2.5)中的未知数x2，得到与原方程组等价的三角形方程组显然，方程组是(2.6)是容易求解的，解为上述过程相当于对方程的增广阵做初等行变换其中ri用表示矩阵的第i行.由此看出，用消去法解方程组的基本思想是用逐次消去未知数的方法把原方程组Ax=b化为与其等价的三角形方程组，而求解三角形方程组可用回代的方法求解.换句话说，上述过程就是用初等行变换将原方程组系数矩阵化为简单形式(上三角矩阵)，从而求解原方程组(2.1)的问题转化为求解简单方程组的问题.或者说，对系数矩阵A施行一些行变换(用一些简单矩阵左乘A)将其约化为上三角矩阵.这就是高斯消去法.下面讨论求解一般线性方程组的高斯消去法.由将(2.1)记为A(1)x=b(1)，其中

(1)第1步(k=1).

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

mi1=ai1(1)/a11(1)(i=2,3,…,m).用-mi1乘(2.1)的第一个方程，加到第i个(i=2,3,,m)方程上，消去(2.1)的从第二个方程到第m个方程中的未知数x1，得到与(2.1)等价的方程组(2.7)简记为A(2)x=b(2)，其中A(2),b(2)的元素计算公式为

(2)第k次消元(k=1,2,,s,s=min(m-1,n)).

设上述第1步，

，第k-1步消元过程计算已经完成，即已计算好与(2.1)等价的方程组，简记为A(k)x=b(k)，其中

设akk(k)0，计算乘数

mik=aik(k)/akk(k)

(i=k+1,,m).用-mik乘(2.8)的第k个方程加到第i个(i=k+1,,m)方程上，消去从第k+1个方程到第m个方程中的未知数xk，得到与(2.1)等价的方程组A(k+1)x=b(k+1)，其中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)为上阶梯.

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

如果A是非奇异矩阵，且akk(k)0(k=1,2,,n-1)，求解三角形方程组(2.10)，得到求解公式(2.10)的求解过程(2.11)称为回代过程.

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

调到(1,1)位置，然后进行消元计算，这时A(2)右下角矩阵为n-1阶非奇异矩阵，继续这过程，高斯消去法照样可进行计算.总结上述讨论即有

定理5

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

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

(a)消元计算(k=1,2,,n-1)

(b)回代计算

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

算法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

.

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

数akk(k)在高斯消去法中有着突出的作用，称为约化的主元素(简称主元).

算法2(回代算法)

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

对于j=n,n-1,,1

(1)xi←bi

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

xi←xi-uij·xj

(3)xi←xi/uii

这个算法需要n(n-1)/2次乘除法运算.例子

解方程组解：消元回代得

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

由此，需要对算法1进行修改，首先研究原来矩阵A在什么条件下才能保证akk(k)0(k=1,2,,s).下面的定理给出了这个条件.

定理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)，即利用行列式的性质，我们有

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

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

定理6得证.

推论

如果A的顺序主子式Di0(i=1,2,,n-1),则5.2.2矩阵三角分解

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

设(2.1)的系数矩阵A∈Rn×n的各顺序主子式均不为零，对于A施行初等行变换相当于用初等矩阵左乘A，于是对(2.1)施行第1次消元后化为(2.7)，这时A(1)化为A(2)，b(1)化为b(2)，即L1A(1)=A(2),L1b(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.14)

将上三角矩阵A(n)记为U，由(2.14)得到其中为单位下三角矩阵.

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

定理7(矩阵的LU分解)

设A为n阶矩阵，如果A的顺序主子式Di0(i=1,2,,n-1)，则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，即A=LU,且这种分解是唯一的.

证明

根据以上高斯消去法的矩阵分析，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.证毕.

例2

对于例1，系数矩阵由高斯消去法m21=0，m31=2，m32=-1，故从而得到求矩阵行列式的计算公式为5.3高斯主元素消元法

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

高斯消去法也称主元素消去法

(条件det

A0)

即当akk(k)=0

时，高斯消元法无法进行；或

|akk(k)|<<1时，带来舍入误差的扩散。（小主元）解

(方法1)高斯消元法(用4位有效数字)例3(小主元产生了大误差)

(方法2)列主元消元法，先交换行

这个例子告诉我们，在采用高斯消去法解方程组时，小主元可能产生麻烦，故应避免绝对值小的主元素akk(k)，对一般矩阵来说，最好每一步选取系数矩阵(或消元后的低阶矩阵)中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数值稳定性，这就是全主元素消去法，在选主元时要花费较多机器时间，目前主要使用的是列主元素消去法，本节主要介绍列主元素消去法，并假定(2.1)的A∈Rn×n为非奇异的矩阵.5.3.1列主元素消去法设方程组(2.1)的增广矩阵为首先在A的第1列中选取绝对值最大的元素作为主元素，例如然后交换B的第1行与第i1行，经第1次消元计算得(A|b)→(A(2)|b(2))

重复上述过程，设已完成第k-1步的选主元素，交换两行及消元计算，(A|b)约化为其中A(k)的元素仍记为aij，b(k)的元素仍记为bi.第k步的选主元素(在A(k)右下角方阵的第1列内选)，即确定ik，使交换(A(k)|b(k))第k行与ik行的元素，再进行消元计算，最后将原方程组化为

(k=1,2,,n-1)

回代求解

算法3(列主元素消去法)

设Ax=b.本算法用A的具有行交换的列主元素消去法，消元结果冲掉A，乘数mij冲掉aij

，计算解x冲掉常数项b，行列式存放在det中.

1.det←1

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

(1)按列选主元

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

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

(a)aik←mik=aik/akk

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

aij←aij-mik·akj

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

换行

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

bk↔bik

det←-det

(c)bi←bi-mik·bk

(5)det←akk·

det

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

4.回代求解

(1)bn←bn/ann

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

5.det←ann·

det

例3的(方法2)用的就是列主元素消去法.

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

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

利用(3.1)得到简记为其中下面就n=4来考察一下矩阵其中由习题3知亦为单位下三角矩阵，其元素的绝对值不超过1.记由(3.2)得到PA=LU.其中P为排列矩阵，L为单位下三角矩阵，U为上三角矩阵.这说明对(2.1)应用列主元消去法相当于对(A|b)先进行一系列行交换后对PAx=Pb再应用高斯消去法.在实际计算中我们只能在计算过程中做行的交换.总结以上的讨论我们有

定理8(列主元素的三角分解定理)

如果A为非奇异矩阵，则存在排列矩阵P使PA=LU，其中L为单位下三角矩阵，U为上三角矩阵.在编程过程中，L元素存放在数组A的下三角部分，U元素存放在数组A的上三角部分，由纪录主行的整型数组Ip(n)可知P的情况.消元法的应用1.求行列式高斯消元法2.求逆矩阵(用高斯-若当消去法)例子

分别用高斯消元法和列主元消元法求矩阵的行列式.解：高斯消元法大数“吃掉”了小数！列主元消元法全主元素消去法介绍在(5.8)中，若每次选主元素不局限在方框列中，而在整个主子矩阵中选取,便称为全主元高斯消去法,此时增加的步骤为:

①

，确定ik

,jk,若=0，给出|Ａ|=0的信息，退出计算.

行交换:(kjn)

列交换:(kin)②作如下行交换和列交换值得注意的是,在全主元的消去法中，由于进行了列交换，x各分量的顺序已被打乱.因此必须在每次列交换的同时，让机器“记住”作了一次怎样的交换，在回代解出后将x各分量换回原来相应的位置，这样增加了程序设计的复杂性.此外作①步比较大小时，全主元消去法将耗用更多的机时.但全主元消去法比列主元消去法数值稳定性更好一些.实际应用中，这两种选主元技术都在使用.选主元素的高斯消元法是一种实用的算法，它可以应用于任意的方程组Ax=b，只要|A|≠0.5.3.2高斯-若当消去法高斯-若当消去法是高斯消去法的一种变形.高斯消去法是消去对角元素下方的元素.如果同时消去对角元素上方和下方的元素，而且将对角元素化为1，这种方法就称为高斯-若当(Gauss-Jordan)消去法.设高斯-若当消去法已完成k-1步，于是Ax=b化为等价方程组A(k)x=b(k)，其中增广阵在第k步计算时(k=1,2,,n),考虑将第k行上、下的第k列元素都进行消元计算化为零，且akk化为1.

1.按列选主元素，即确定ik使

2.换行(当ik≠k时)，交换第k行与第ik行元素(j=k,k+1,,n),

3.计算乘数mik=-aik/akk(i=1,2,,n,i≠k)

mkk=1/akk.

(mik可保存在存放aik的单元中).

4.消元计算

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

bi←bi+mikbk

(i=1,2,,n,且i≠k)

5.计算主行

akj←akjmkk

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

bk←bkmkk.上述过程结束后，即当k=n时有显然xi=bi,i=1,2,,n就是Ax=b的解.说明用高斯-若当消去法虽比高斯消去法略复杂些，但将A约化为单位矩阵，计算解就在常数项位置得到，因此用不着回代求解，故也称为无回代的高斯消元法.它的计算量大约需要n3/2次乘除法，要比高斯消去法大，但用高斯-若当消去法求一个矩阵的逆矩阵还是比较合适的.

定理9(高斯-若当法求逆矩阵)

设A为非奇异矩阵，方程组AX=In的增广矩阵为C=(A|In)，如果对C应用高斯-若当方法化为(In|T)

，则A-1=T.事实上,求A的逆矩阵A-1,即求n阶矩阵X使AX=In,其中In为单位矩阵，将X按列分块X=(x1,x2,,xn)，I=(e1,e2,,en)，于是求解AX=In等价于求解n个方程组Axj=ej(j=1,2,,n).我们可用高斯-若当方法求解AX=In.

例4用高斯-若当方法求下列矩阵的逆矩阵.

解由分块矩阵所以得参看书上p181-例4.将高斯-若当消元法算法略加改造便成了求逆的算法.求解方程组

AX=I

对k=1,2,,n

1.按列选主元,确定ik

2.换行,交换第k行和ik行(j=k,k+1,,n)(j=1,2,,n)

3.计算乘数mik=-aik／akk

,(i=1,2,,i≠k)

mkk=1／akk,

5.计算主行

akj←akjmkk

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

xkj←xkjmkk

(j=1,2,…,n)应该注意到，在上述算法中都加进了选列主元的步骤，它可以保证在方阵A非奇异的情况下都能求得A-1.求逆算法中不宜使用全面选主元，否则会改变逆阵元素的排列，增加程序设计的复杂性.

4.消元

aij←aij+mik

akj,(i=1,2,…,n,且i≠k,j=k+1,…,n)

xij←xij+mik

xkj,(i=1,2,…,n,且i≠k,j=1,2,…,n)得A-1=(xij)n×n.5.4矩阵三角分解法高斯消去法有很多变形，有的是高斯消去法的改进，改写，有的是用于某一类特殊矩阵的高斯消去法的简化.下面我们将介绍矩阵的直接三角分解法，解特殊方程组用的平方根法及追赶法.

定义

如果

L为单位下三角阵，U为上三角阵，则称A=LU为杜里特尔(Doolittle)分解；如果

L为下三角阵，U为单位上三角阵，则称A=LU为克劳特(Crout)分解.5.4.1直接三角分解法(LU分解)

在5.2.2已经通过高斯消去法得到一个将A分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，A=LU，其中并由定理7得到这种分解是唯一的.

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

①Ly=b，求y；

②Ux=y，求x.

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

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

a1n

u11u12u1k

u1n

第1框a21a22a2k

a2n

l21u22u2k

u2n第2框

ak1ak2

akk

akn

lk1lk2

ukk

ukn

第k框

an1an2ank

ann

ln1ln2

lnk

unn

第n框比较式A=LU

两端的元素,按下图所示顺序逐框进行，先求ukj，后求lik.由第一框可得得假设前k-1框元素已求出，则由有了矩阵A

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

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

Ly=b

与上三角方程组Ux=y

其计算公式如下：

Ly=b

Ux=y例5

求矩阵解用紧凑格式的LU(Doolittle)分解.

u11=2

u12=1u13=4

所以u11=2u12=1u13=4

参看书上p185-例5.得行列式

2.选主元的三角分解法

见书p186讲解，自己看为主.5.4.2平方根法实际问题中Ax=b，系数矩阵A大多是对称正定矩阵，所谓平方根法，就是利用对称正定矩阵的三角分解而得到求解对称正定方程组的一种有效方法，目前在计算机上广泛应用平方根法解此类方程组.

对称矩阵的LDLT

分解法当A为对称矩阵,且其顺序主子式均不为0时,由定理7可知A有唯一的LU分解为(4.1)式.为了利用A的对称性，将U再分解为其中D为对角矩阵，U0为上三角矩阵，于是

A=LU=LDU0.(4.6)又A=AT=U0T(DLT),由分解的唯一性即得U0T=L.代入(4.6)得到对称矩阵A的分解式A=LDLT.总结上述讨论有

定理10(对称矩阵的三角分解定理)

设A为n阶对称矩阵，且A的各阶顺序主子式均不为零，则A可以唯一分解为A=LDLT,其中L是单位下三角矩阵，D是对角矩阵.如果A为对称正定矩阵，则A的分解式A=LDLT中D的对角元素di均为正数.事实上,由A的对称正定性,5.2节中推论成立,即于是有由定理6得到其中L1=LD1/2为下三角矩阵.

定理11(对称正定矩阵的三角分解或乔雷斯基(Cholesky)分解)

如果A为n阶对称正定矩阵，则存在一个实的非奇异下三角矩阵L使A=LLT，当限定L的对角元素为正时，这种分解是唯一的.下面我们用直接分解方法来确定计算L元素的递推公式，因为其中lii>0(i=1,2,,n).由矩阵乘法及ljk=0(当j<k时),得于是得到解对称正定方程组Ax=b的平方根法计算公式对于j=1,2,,n求解Ax=b，即求解两个三角形方程组

(1)Ly=b，求y；(2)LTx=y，求x.由计算公式知道所以有于是上面分析说明，分析过程中元素ljk的数量级不会增长且对角元素ljj恒为正数.于是有结论不选主元素的平方根法是一个数值稳定的方法.当求出L的第j列元素时，LT的第j行元素亦算出.所以平方根约需n3/6次乘除法，大约为一般直接LU分解法计算量的一半.例题用平方根法求解对称正定方程组解首先对A进行Cholesky分解求解Ly=b，得y1=2,y2=3.5,y3=1.求解LTx=y，得x1=1,x2=1,x3=1.由公式(4.7)看出，用平方根法解对称正定方程组时，计算L的元素ljj需要用到开方运算.为了避免开方，我们下面用定理10的分解式A=LDLT.即由矩阵乘法，并注意ljj=1,ljk=0(j<k)，得于是得到L的元素及D的对角元素公式：为了避免重复计算，我们引进tij=lij

dj

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

求解Ly=b及DLTx=y的计算公式为计算公式(4.10),(4.11)称为改进的平方根法.5.4.3追赶法在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条插值函数，都会要求解系数矩阵为对角占优的三对角线方程组Ax=f，即:其中,当|i-j|>1时,aij=0,且满足如下的对角占优条件:(a)|b1|>|c1|>0;(b)|bi|≥|ai|+|ci|,ai,ci≠0,(i=2,3,…,n-1).(c)|bn|>|an|>0.我们利用矩阵的直接三角分解法来推导解三对角线方程组(4.12)的计算公式.由系数阵A的特点，可以将A分解为两个三角矩阵的乘积，即A=LU,其中取L下三角阵,取U为单位上三角阵,这样求解方程组Ax=f的方法称为追赶法.设其中为待定系数，比较(4.13)两边即得求解Ax=f等价于求解两个三角形方程组.

(1)Ly=f,求y;(2)Ux=y,求x.从而得到解三对角线方程组的追赶法.我们将计算系数及的过程称为追的过程,将计算方程组的解的过程称为赶的过程.合起来就是追赶法.追赶法公式实际上就是把高斯消去法用到求解三对角线方程组上去的结果.这时由于A特别简单,因此使得求解的计算公式非常简单,而且计算量仅为5n-4次乘除法,而另外增加一个方程组Ax=f2仅增加3n-2次乘除法运算,易见追赶法的计算量是比较小的.总结上述讨论有

定理12

设有三对角线方程组Ax=f,其中A满足条件(a),(b),(c),

则A为非奇异矩阵且追赶法计算公式中的满足由定理12的10,20说明追赶法计算公式中不会出现中间结果数量级的巨大增长和舍入误差的严重累积.5.5向量与矩阵的范数在分析方程组的解的误差及下章中迭代法的收敛性时，常产生一个问题,即如何判断向量x

的“大小”，对矩阵也有类似的问题.本节介绍n维向量范数和n×n矩阵的范数.向量范数是三维欧氏空间中向量长度概念的推广，在数值分析中起着重要作用.首先将向量长度概念推广到Rn(或Cn)中.称为向量x,y的数量积(内积).将非负实数称为向量x的欧氏范数.

定义1设x=(x1,x2,,xn)T,y=(y1,y2,,yn)TRn

,将实数下面定理可在线性代数书中找到.

定理13设x,yRn

(或Cn),则

1.(x,x)=0,当且仅当x=0时成立;我们还可以用其他办法来度量Rn中向量的“大小”.例如,对于x=(x1,x2)TRn,可以用一个x的函数N(x)=max|xi|来度量x的“大小”,而且这种度量x的“大小”的方法计算起来比欧氏范数方便.在许多应用中,对度量x的“大小”的函数N(x)都要求是正定的、齐次的且满足三角不等式.下面我们给出向量范数的一般定义.

(1)||x||0(||x||=0当且仅当x=0)(正定性),(2)||αx||=|α|

||x||,

对任何αR(或αC)(齐次性),(3)||x+y||||x||+||y||(三角不等式).则称N(x)=||x||是Rn(或Cn)上的一个向量范数(或模).定义2(向量的范数)如果向量xRn(或Cn)的某个实值函数N(x)=||x||，满足条件:由(3)可推出不等式.

(4)|||x||-||y|||||x-y||.下面给出几种常用的向量范数,设x=(x1,x2,,xn)T.容易证明前三种范数是的p-范数特殊情况,其中向量的-范数(最大范数)向量的1-范数向量的p-范数(0<p<+)向量的欧氏范数

例6

计算向量

x=(1,-2,3)T的各种范数.

解定义3设{x(k)}为Rn中一向量序列,x*Rn,记x(k)=(x1(k),,xn(k))T,x*=(x1*,,xn*)T.如果则称x(k)收敛于x*.记为定理14(向量范数的连续性)设非负函数N(x)=||x||为Rn上任一向量范数,则N(x)是x的分量x1,x2,,xn的连续函数.

证明设其中只须证明x→y时N(x)→N(y)即成.事实上即定理15(向量范数的等价性)设||x||s,||x||t为Rn上任意的两种范数,则存在常数c1,c2>0,使得对一切xRn有

证明只要就||x||s=||x||∞证明上式成立即可,即证明存在常数c1,c2>0,使得考虑n元函数记S={x|||x||∞=1,xRn},则S是一个有界闭集.由于f(x)为S上的连续函数,所以f(x)于S上达到最大值和最小值,即存在x′,x″S使得设xRn且x≠0,则x/||x||∞S,从而有即可以证明常用范数满足下述关系以及

注意定理15不能推广到无穷维空间.由定理15可得结论:如果在一种范数意义下向量序列收敛时,则在任何一种范数意义下该向量序列均收敛.定理16,其中||·||为向量的任一种范数.

证明显然,而对于Rn上任一种范数||·||,由定理15,存在常数c1,c2>0,使于是又有下面我们将向量范数概念推广到矩阵上去.可以将Rn×n中的矩阵A=(aij)nn当作n2维向量,则由向量的2-范数可以得到Rn×n中矩阵的一种范数称为A的Frobenius范数.||A||F显然满足正定性、齐次性及三角不等式.下面我们给出矩阵范数的一般定义.定义4(矩阵的范数)如果矩阵ARn×n的某个实值函数N(A)=||A||，满足条件:

(1)||A||

0(||A||=0当且仅当A=0)(正定性)；

(2)||cA||=|c|||A||,c为实数(齐次性)；

(3)对任意A,B有||A+B||||A||+||B||(三角不等式)

(4)对任意A,B有||AB||||A||||B||(相容性)；则称N(A)

是Rn×n上的一个矩阵范数(或模).上面我们定义的F(A)=||A||F就是Rn×n上的一个矩阵范数.上述条件(即不等式)就称为矩阵范数与向量范数的相容性条件.由于在大多数与估计有关的问题中,矩阵和向量会同时参与讨论,所以希望引进一种矩阵的范数,它是和向量范数相联系而且和向量范数相容的,即对任何向量xRn及矩阵ARn×n都成立||Ax||≤||A||||x||.为此我们再引进一种矩阵的范数.定义5(矩阵的算子范数)设xRn,ARn×n,给出一种向量范数||x||v(如v=1,2或),相应地定义一个矩阵的非负函数可验证||A||v满足定义4(见下面定理),所以||A||v是Rn×n上矩阵的一个范数,称为A的算子范数.定理17设||x||v是Rn上的一个向量范数,则||A||v是Rn×n上矩阵的范数,且满足相容条件证明由(5.6)相容条件(5.7)是显然的.现只验证定义4中条件(4).由(5.7),有当x≠0,有故显然这种矩阵的范数||A||v依赖于向量范数||x||v的具体含义.也就是说,当给出一种具体的向量范数||x||v时,相应地就得到了一种矩阵范数||A||v.定理18设xRn,A=(aij)Rn×n,则有常用的算子范数（称为A的行范数）,（称为A的列范数）,（称为A的2-范数）.其中max(ATA)表示ATA的最大特征值,由于矩阵2-范数与ATA

的特征值有关，所以又称为A的谱范数.证明见书p202.由定理18看出,计算一个矩阵的||A||∞,||x||1还是比较容易的,而矩阵的2-范数在||x||2计算上不方便,但是矩阵的2-范数具有许多好的性质,它在理论上是非常有用的.我们指出,对于复矩阵(即ACn×n)定理18中1,2显然也成立,对于3应改为例7设则求相应的各种范数.解因为令得ATA的特征值故可见参见书p203-例7.

定义6

设ARn×n的特征值为i(i=1,2,,n),称为A的谱半径.例子

设,求A的谱半径.解得A的特征值故得A的谱半径为定理19(特征值上界)设ARn×n,则ρ(A)≤||A||,即A的谱半径不超过A的任何一种算子范数(对||A||F亦成立).证明设是A的任一特征值,x为相应的特征向量,则Ax=x,由(5.7)得注意到||x||≠0,即得||≤||A||.定理20如果ARn×n为对称矩阵,则||A||2=ρ(A).

(证明留作练习题或查有关参考书.)定理21如果||B||<1,则I±B为非奇异矩阵,且其中||·||是指矩阵的算子范数.证明用反证法.若det(I-B)=0,则(I-B)x=0有非零解,即存在x0≠0使Bx0=x0,||Bx0||/||x0||=1,故||B||≥1,这与假设矛盾.又由(I-B)(I-B)-1=I,有(I-B)-1=I+B(I-B)-1,从而类似对加法有(I+B)(I+B)-1=I,(I+B)-1=I-B(I+B)-1,得求‖A‖1,‖A‖2,‖A‖∞,ρ(A).例子

设,解:

显然‖A‖1=4,‖A‖∞=4

这里,我们指出,对于实对称矩阵A,有ρ(A)=||A||2.5.6误差分析5.6.1矩阵的条件数由于A(或b)元素是测量得到的,或者是计算的结果,在第一种情况A(或b)常带有某些观测误差,在后一种情况A(或b)又包含有舍入误差.因此我们处理的实际矩阵是A+A(或b+b),下面我们来研究数据A(或b)的微小误差对解的影响.即考虑估计x-y,其中y是(A+A)y=b的解.考虑线性方程组Ax=b,其中设A为非奇异矩阵,x为方程组的精确解.首先考察一个例子.

例8

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

.现在考虑常数项的微小变化对方程组解的影响,即考察方程组也可表示为A(x+x)=b+b,其中b=(0,0.0001)T,y=x+x,x为(6.1)的解.显然方程组(6.2)的解为x+x=(1,1)T.我们看到(6.1)的常数项b的第2个分量只有1/1000的微小变化,方程组的解却变化很大.这样的方程组称为病态方程组.

定义7

如果矩阵A或常数项b

的微小变化(小扰动),引起方程组Ax=b解的巨大变化,则称此方程组为“病态”方程组,其系数矩阵A

称为“病态”矩阵(相对于方程组而言),否则称方程组为“良态”方程组,A称为“良态”矩阵.应该注意,矩阵的“病态”性质是矩阵本身的特性,下面我们希望找出刻画矩阵“病态”性质的量.设有方程组

Ax=b,(6.3)其中A为非奇异矩阵,x为(6.3)的精确解.以下我们研究方程组的系数矩阵A(或b)的微小误差(小扰动)时对解的影响.

(1)

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

x+x,则

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

||x||≤||A-1||||b||.(6.4)由(6.3)有||b||≤||A||||x||.于是由(6.4)及(6.5),得定理22设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.(6.6)如果A不受限制的话,可能A+A奇异,而(A+A)=A(I+A-1A).由定理21知,||A-1A||<1时,(I+A-1A)-1存在.由(6.6)式x=-(I+A-1A)-1A-1(A)x.因此设||A-1||||A||<1,即得定理23设A是非奇异矩阵,Ax=b≠0

,且

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

(3)

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

x+x,则还可以推出相对误差估计式为定义8设A是非奇异矩阵,称数cond(A)v=||A-1||v||A||v(v=1,2或∞)为矩阵A的条件数

.由此看出矩阵的条件数与范数有关.矩阵的条件数是一个十分重要的概念.由上面讨论知,当A的条件数相对的大,即cond(A)>>1时,则(6.3)是“病态”的(即A是“病态”矩阵,或者说A是坏条件的,相对于方程组),当A的条件数相对的小,则(6.3)是“良态”的(或者说A是好条件的).注意,方程组病态性质是方程组本身的特性.A的条件数越大,方程组的病态程度越严重,也就越难用一般的计算方法求得比较准确的解.

例如对前面例8的矩阵作分析由于条件数cond(A)1很大,可见矩阵A的病态程度十分严重,故由此方程组的解误差非常大.计算其条件数

通常使用的条件数,有

(1)cond(A)∞

=||A-1||∞||A||