数值计算方法第三章讲义

第3章线性方程组的数值解法三、矩阵的三角分解及其应用四、线性代数方程组的性态与误差分析五、迭代法一、线性方程组二、高斯(Gauss)消元法六、共轭梯度法求解上一页下一页返回

一、线性方程组实际问题中的线性方程组分类：按系数矩阵中零元素的个数：稠密线性方程组稀疏线性方程组按未知量的个数：高阶线性方程组低阶线性方程组(如1000)(80%)按系数矩阵的形状对称正定方程组三角形方程组三对角占优方程组

实际中，存在大量的解线性方程组的问题。很多数值方法到最后也会涉及到线性方程组的求解问题：如样条插值的M和m关系式，曲线拟合的法方程，方程组的Newton迭代等问题。解线性方程组的两类数值方法：直接法:经过有限次四则运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法。（一般有限步内得不到精确解）特点：准确，可靠，理论上得到的解是精确的。适合于中小规模的方程组。特点：速度快，但有误差。适合于大规模的方程组。Cramer法则是一种直接法，但无法实用

直接法概述直接法是将原方程组化为一个或若干个三角形方程组的方法．对于线性方程组根据Cramer(克莱姆)法则,若若用初等变换法求解,则对其增广矩阵作行初等变换:经过n-1次同解即以上求解线性方程组的方法称为Gauss消元法则都是三角形方程组上述方法称为直接三角形分解法不论是Gauss消元法还是直接三角形分解法,最终都归结为解三角形方程组一、三角形线性代数方程组的直接法若记下三角形线性方程组上三角形线性方程组Gauss消去法二、直接法——即按前推过程求解其解为其解为:按回代过程求解对方程组，作如下的变换，解不变：①交换两个方程的次序②一个方程的两边同时乘以一个非0的数③一个方程的两边同时乘以一个非0数，加到另一个方程二、Gauss消元法因此，对增广矩阵(A，b)作如下的变换，解不变：①交换矩阵的两行②某一行乘以一个非0的数③某一行乘以一个非0数，加到另一行消元法就是对增广矩阵作上述行的变换，变为三角形线性方程组，而后求根。

Gauss消元法的基本思想：用矩阵表示：=首先将A化为上三角阵，再回代求解。例3.1.1解化为同解方程组化为上三角方程组按回代过程求解上述求解三元线性代数方程组的方法即为Gauss消元法。

Gauss消元法计算过程矩阵表示记为：对其增广矩阵施行行初等变换(1)定义行乘数相当于第i个方程-第一个方程×行乘数→新的第i方程—同解！第一个方程不动！

且(1)定义行乘数相当于第i个方程-第二个方程×行乘数→新的第i方程—同解！第一、二个方程不动！

(1)定义行乘数相当于第i个方程-第k个方程×行乘数→新的第i方程—同解！第一至第k个方程不动！

系数矩阵与常数项：

Gauss消元法的运算量乘法次数：除法次数：全部回代过程需作乘除法的总次数为于是Gauss消元法的乘除法运算总的次数为数级类似地，可得Gauss消元法的加减法运算总的次数为由于计算机完成一次乘除运算所耗时间要远远多于完成一次加减运算的时间，而且按统计规律，在一个算法中，加减运算和乘除运算次数大体相当，故在衡量一个算法的运算量时只需统计乘除的运算次数。Gauss消元法乘除法约为2700次而如果用Cramer法则的乘除法运算次数约为用行列式定义回代求解从Gauss消元法的过程自然想到：如果在消元过程中把对角线上方未知量的系数也化为零，使方程组化成每个方程只含一个未知量的方程组，则无须回代就可得到解。这种消元法称为Gauss－Jordan消元法。但经计算，它需要的乘除法次数为显然，其运算量超过Gauss消元法。所以用Gauss－Jordan消元法求解线性代数方程组是不可取的。不过，用它来求逆矩阵却有方便之处。三、选主元Gauss消元法例3.1.3用Gauss消元法解线性方程组(用3位十进制浮点数计算）解本方程组的精度较高的解为用Gauss消元法求解(用3位十进制浮点数计算)回代后得到主元－9999与精确解相比,该结果相当糟糕！究其原因,在求行乘数时用了很小的数0.0001作除数!

上述例题表明：Gauss消元法是不稳定的算法，其中小主元是不稳定的根源。因而在Gauss消元法中要避免小主元的出现。这就需要采用“选主元”的技术。所谓选主元，简单地说就是选取绝对值最大的元素作主元。

全选主元Gauss消元法在此范围内选主元需进行元素比较

列选主元Gauss消元法在此范围内选主元不必选主元的情况：当系数矩阵A对称正定或严格对角占优或不可约对角占优时，可不必选主元。现用列选主元消元法重解例3.1.3：将方程组的1,2行交换，即得回代后得到这是一个相当不错的结果！0.9999例3.1.4解线性方程组(用8位十进制尾数的浮点数计算)解这个方程组和例3.1.3一样，若用Gauss消元法计算会有小数作除数的现象，若采用换行的技巧（即列选主元消元法），则可避免。绝对值最大不需换行经过回代后可得事实上,方程组的准确解为

列选主元Gauss消元法的计算步骤?

列选主元Gauss消元法的流程图开始输出无解信息消元换行停机回代求解本算法主要由四个模块组成一些特殊情况,主元就在对角线上,不需选主元.元素满足如下条件的矩阵即对角线上每一元素的绝对值均大于同行其他各元素绝对值之和,这样的矩阵称为按行严格对角占优矩阵,简称严格对角占优矩阵.例：性质：严格对角占优矩阵必定非奇异.上一页下一页

返回

补充一、LU分解法

/*LU

Factorization.*/就n=3的情况分析顺序消去法的消元过程.上一页下一页

返回

矩阵的三角(LU)分解法三、直接法—可见,消元过程相当于下述矩阵乘法运算：由分块乘法可得：直接计算可得由(*)式可得上一页下一页

返回

Step1:记M(1)=，则Stepn

1:其中M(k)=

以上分析推广到n阶线性方程组可得上一页下一页

返回

记为L单位下三角阵/*unitarylower-triangularmatrix*/记

U=A

的

LU

分解/*LUfactorization*/也称A

的Doolittle分解上一页下一页

返回

道立特分解法/*DoolittleFactorization*/：——LU

分解的紧凑格式/*compactform*/根据矩阵乘法法则,先比较等式两边第1行和第1列元素有：上一页下一页

返回

设已定出U的第1行到第k-1行的元素L的第1列到第k-1列的元素上一页下一页

返回

(1),(2)两式就是LU分解的一般计算公式,其结果与高斯消去法所得结果完全一样.避免了中间过程的计算,故称为A的直接分解公式.上一页下一页

返回

按上图逐框求出矩阵A的LU分解,这种计算方法称为紧凑格式法。上一页下一页

返回

定理若矩阵A非奇异,则A能分解为LU的充分必要条件是A的所有顺序主子式均不为0.定理若非奇异矩阵A有LU分解,则此分解是唯一的.上一页下一页

返回

矩阵的Doolittle分解法的Matlab程序function[1,u]=lu_Doolittle1(A)%A为可逆方阵，l为返回下三角矩阵，u为上三角阵n=length(A);u=zeros(n);l=eye(n);u(1,:)=A(1,:);l(2:n,1)=A(2:n,1)/u(1,1);fork=2:nu(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);end例：利用系数矩阵的LU分解,求解方程组解：LU分解的紧凑格式为上一页下一页

返回

在Matlab命令窗口执行>>A=[1111;1213;4321;61137];>>[l,u]=lu_Doolittle1(A)则求解原方程组等价于求解下面两个方程组Ly=b及Ux=y上一页下一页

返回

在Matlab命令窗口执行>>A=[1111;1213;4321;61137];>>b=[6121445]’;[y,x]=LU_s(A,b)A=LU若U为单位上三角矩阵，则称为Crout分解.Crout分解相应的求解公式可类似得到。上一页下一页

返回

二、平方根法——对称

正定矩阵的分解法将对称

正定阵

A做LU

分解记为定理设n阶对称正定矩阵A，则存在唯一的单位下三角阵L及对角阵D使得。称为矩阵A的乔累斯基分解上一页下一页

返回

定理设矩阵A对称正定，则存在唯一的对角元为正的下三角阵L，使得。称为对称正定矩阵A的乔累斯基分解利用乔累斯基（Cholesky）分解式来求解Ax=b的方法也称Cholesky方法或平方根法三、三对角方程组追赶法若A满足Gauss消去法可行的条件，则可用LU分解法求解其中：上一页下一页

返回

解方程组Ax=d分为两步，即求解Ly=d和Ux=y，计算公式如下：上述方法为求解三对角方程组的追赶法，也称Thomas算法.上一页下一页

返回

追赶法公式简单，计算量和存储量都小，整个求解过程只需要5n-4次乘除运算。在分析方程组的解的误差及下章中迭代法的收敛性时，常产生一个问题,即如何判断向量x

的“大小”，对矩阵也有类似的问题.本节介绍n维向量范数和n×n矩阵的范数.向量范数是三维欧氏空间中向量长度概念的推广，在数值分析中起着重要作用.四、线性代数方程组的性态与误差分析首先将向量长度概念推广到Rn(或Cn)中.称为向量x,y的数量积(内积).将非负实数称为向量x的欧氏范数.定义3.1设x=(x1,x2,

,xn)T,y=(y1,y2,

,yn)T

Rn

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

定理设x,y

Rn

(或Cn),则1.(x,x)=0,当且仅当x=0时成立;我们还可以用其他办法来度量Rn中向量的“大小”.例如,对于x=(x1,x2)T

Rn,可以用一个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||(三角不等式).则称f(x)=||x||是Rn(或Cn)上的一个向量范数(或模).定义3.2(向量的范数)如果向量x

Rn(或Cn)的某个实值函数f(x)=||x||，满足条件:由(3)可推出不等式.(4)|||x||-||y|||

||x-y||.下面给出几种常用的向量范数,设x=(x1,x2,

,xn)T.容易证明前三种范数是的p-范数特殊情况,其中向量的

-范数(最大范数)向量的1-范数向量的p-范数(0<p<+)向量的欧氏范数,2范数显然并且由于

例计算向量x=(1,-2,3)T的各种范数.

解定义3.3设{x(k)}为Rn中一向量序列,x*

Rn,记x(k)=(x1(k),

,xn(k))T,x*=(x1*,

,xn*)T.如果则称x(k)收敛于x*.记为定理3.1(向量范数的连续性)设函数f(x)=||x||为Rn上任一向量范数,则f(x)是x的分量x1,x2,

,xn的连续函数.

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

Rn有

证明只要就||x||s=||x||∞证明上式成立即可,即证明存在常数c1,c2>0,使得考虑n元函数记S={x|||x||∞=1,x

Rn},则S是一个有界闭集.由于f(x)为S上的连续函数,所以f(x)于S上达到最大值和最小值,即存在x′,x″

S使得设x

Rn且x≠0,则x/||x||∞

S,从而有即可以证明常用范数满足下述关系以及

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

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

Rn×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上的一个矩阵范数.上述条件(即不等式)就称为矩阵范数与向量范数的相容性条件.由于在大多数与估计有关的问题中,矩阵和向量会同时参与讨论,所以希望引进一种矩阵的范数,它是和向量范数相联系而且和向量范数相容的,即对任何向量x

Rn及矩阵A

Rn×n都成立||Ax||≤||A||||x||.为此我们再引进一种矩阵的范数.定义3.5(矩阵的算子范数)设x

Rn,A

Rn×n,给出一种向量范数||x||v(如v=1,2或

),相应地定义一个矩阵的非负函数可验证||A||v满足定义4(见下面定理),所以||A||v是Rn×n上矩阵的一个范数,称为A的算子范数.定理3.4设||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.定理3.5设x

Rn,A=(aij)

Rn×n,则有常用的算子范数（称为A的行范数）,（称为A的列范数）,（称为A的2-范数）.其中

max(ATA)表示ATA的最大特征值,由于矩阵2-范数与ATA的特征值有关，所以又称为A的谱范数.定义3.6设A

Rn×n的特征值为

i(i=1,2,

,n),称为A的谱半径.例子设,求A的谱半径.解得A的特征值故得A的谱半径为由定理3.5看出,计算一个矩阵的||A||∞,||x||1还是比较容易的,而矩阵的2-范数在||x||2计算上不方便,但是矩阵的2-范数具有许多好的性质,它在理论上是非常有用的.我们指出,对于复矩阵(即A

Cn×n)定理3.5中1,2显然也成立,对于3应改为例7设则求相应的各种范数.解因为令得ATA的特征值故可见定理3.7(特征值上界)设A

Rn×n,则ρ(A)≤||A||,即A的谱半径不超过A的任何一种算子范数(对||A||F亦成立).证明设

是A的任一特征值,x为相应的特征向量,则Ax=

x,由(5.7)得注意到||x||≠0,即得||≤||A||.定理3.6,其中||·||为矩阵的任一种范数.定理3.8定理3.9如果||A||<1,则I±A为非奇异矩阵,且其中||·||是指矩阵的算子范数.证明用反证法.若det(I-A)=0,则(I-A)x=0有非零解,即存在x0≠0使Ax0=x0,||Ax0||/||x0||=1,故||A||≥1,这与假设矛盾.又由(I-A)(I-A)-1=I,有(I-A)-1=I+A(I-A)-1,从而类似对加法有(I+A)(I+A)-1=I,(I+A)-1=I-A(I+A)-1,得求‖A‖1,‖A‖2,‖A‖∞,ρ(A).例

设,解:显然‖A‖1=4,‖A‖∞=4这里,我们指出,对于实对称矩阵A,有ρ(A)=||A||2.例.求矩阵A的各种常用范数解:由于特征方程为容易计算计算较复杂对矩阵元素的变化比较敏感不是从属范数较少使用使用最广泛性质较好注：

Frobenius范数不是算子范数。

我们只关心有相容性的范数，算子范数总是相容的。

即使A中元素全为实数，其特征根和相应特征向量/*eigenvector*/仍可能是复数。将上述定义中绝对值换成复数模均成立。若不然，则必存在某个向量范数||·||v使得对任意A成立。Counterexample?

谱半径/*spectralradius*/定义矩阵A的谱半径记为

(A)=，其中

i

为

A的特征根。ReIm

(A)定理对任意算子范数||·||有证明：由算子范数的相容性，得到将任意一个特征根

所对应的特征向量代入定理若A对称，则有证明：A对称若

是A的一个特征根，则

2必是A2的特征根。又：对称矩阵的特征根为实数，即

2(A)为非负实数，故得证。对某个

A的特征根

成立所以2-范数亦称为谱范数。矩阵的条件数由于A(或b)元素是测量得到的,或者是计算的结果,在第一种情况A(或b)常带有某些观测误差,在后一种情况A(或b)又包含有舍入误差.因此我们处理的实际矩阵是A+A(或b+b),下面我们来研究数据A(或b)的微小误差对解的影响.即考虑估计x-y,其中y是(A+A)y=b的解.考虑线性方程组Ax=b,其中设A为非奇异矩阵,x为方程组的精确解.二、方程组的性态条件数与摄动理论例，考查以下三个方程组及其准确解

其准确解其准确解其准确解可以看到，后两个方程组与第一个方程组相比，系数矩阵或右端向量仅有0.0005以下的误差，但准确解却相差很大。对这样的方程组，无论用多么稳定的算法求解，一旦计算中产生误差就使解面目全非，所以该方程组的性态很差。上一页下一页

返回

定义：若方程组Ax=b的系数矩阵A与右端向量b的微小变化（小扰动）,将引起解向量x产生巨大变化，则称此方程组为病态方程组，其系数矩阵A称为病态矩阵，否则称Ax=b为良态方程组，称A为良态矩阵.方程组的病态程度与Ax=b对A和b的扰动的敏感程度有关。求解时，A和的误差对解有何影响？

设A精确，有误差，得到的解为，即绝对误差放大因子又相对误差放大因子上一页下一页

返回

设精确，A有误差，得到的解为，即(只要A充分小，使得是关键的误差放大因子，称为A的条件数，记为cond(A),越则A越病态，难得准确解。大上一页下一页

返回

定义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的条件数越大,方程组的病态程度越严重,也就越难用一般的计算方法求得比较准确的解.注:

cond(A)的具体大小与||·||的取法有关，但相对大小一致。

cond(A)取决于A，与解题方法无关。常用条件数有：cond

1(A)cond

(A)cond2(A)条件数的性质：

A可逆，则cond

p

(A)1；

A可逆，

R

则cond(

A)=cond(A)；

A正交，则cond

2

(A)

=1；

A可逆，R正交，则cond

2

(RA)

=cond

2

(AR)=cond(A)2。上一页下一页

返回

通常使用的条件数,有(1)cond(A)∞

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

;(2)A的谱条件数;当A为对称矩阵时其中λ1,λn为A的绝对值最大和绝对值最小的特征值.

条件数的性质:1.对任何非奇异矩阵,都有cond(A)v≥1.事实上

cond(cA)v=cond(A)v.2.设A为非奇异矩阵且c≠0(常数),则3.如果A为正交矩阵,则cond(A)2=1;如果A为非奇异矩阵,R为正交矩阵,则

cond(RA)2=cond(AR)2=cond(A)2.精确解为例计算cond

2(A)

。A1=解：考察A的特征根39206>>1

测试病态程度：给一个扰动，其相对误差为此时精确解为2.0102>200%上一页下一页

返回

例：Hilbert阵cond

(H2)

=27cond

(H3)

748cond

(H6)

=2.9106cond

(Hn)

asn

注：一般判断矩阵是否病态，并不计算A1，而由经验得出。

行列式的值很大或很小（如某些行、列近似相关）；

元素间的数量级相差大，且无规则；

主元消去过程中出现小主元；

特征值相差大数量级。上一页下一页

返回

直接法:

经过有限次运算后可求得方程组精确解的方法(不计舍入误差!)直接法比较适用于中小型方程组。直接法得到的解理论上是准确的，但它们的计算量都是n3数量级，存储量为n2数量级，这在n比较小的时候还比较合适（n<400）。对于现在的很多实际问题，往往要我们求解n很大的高阶方程组，而且这些矩阵往往是稀疏的。对于这类矩阵，用直接法时在运算中很难保持稀疏性，程序设计复杂，且会耗费大量的时间和存储单元。因此我们有必要引入一类新的方法：迭代法。四、迭代法

迭代法能保持矩阵的稀疏性，具有计算简单，存储量小，编制程序容易等优点，并在许多情况下收敛较快。故迭代法能有效地解一些高阶方程组。

迭代法的基本思想是构造一串收敛到解的序列，即建立一种从已有近似解计算新的近似解的规则。由不同的计算规则得到不同的迭代法，本节介绍的几种迭代法都各有其限定的应用范围，这是因为对于一个给定的方程组，某些迭代法收敛得很快，而有些迭代法可能不收敛，或收敛得很慢，以至于无实用价值。迭代法:从解的某个近似值出发,通过构造一个无穷序列去逼近精确解的方法。(一般有限步内得不到精确解)——迭代法的一般形式及其收敛性四、迭代法——几种古典迭代算法（一）雅可比(Jacobi)迭代法五、迭代法建立迭代格式：记

Jacobi迭代的矩阵形式易知，Jacobi迭代有

Jacobi迭代法的计算过程Jacobi迭代法分量形式的Matlab程序function[x,k]=jacobif(A,b,x0,ep,Nmax)n=length(A);k=0;ifnargin<5Nmax=500;endifnargin<5ep=1e-5;endifnargin<3x0=zeros(n,1);y=zeros(n,1);endx=x0;x0=x+2*ep;whilenorm(x0-x,inf)>ep&k<Nmax,k=k+1;x0=x;fori=1:ny(i)=b(i);