数值分析知识内容 (9).pdf

《数值分析知识内容 (9).pdf》由会员分享，可在线阅读，更多相关《数值分析知识内容 (9).pdf（14页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 3.5 方程组的误差分析 为了对线性方程组的直接法作出误差分析，为了讨论方程组迭代法的收敛性，需要向量空间nR中向量（或nnR中矩阵）及向量序列极限的大小，而向量的大小需要引进范数的概念来度量.nR中向量范数是3R中向量长度概念的推广，在数值分析中起着重要作用.3.5.1 向量序列的极限 【定义 1】(向量序列的极限)在n维向量空间nR（或nC）中，设)(kx为向量序列，记为 nTknkkkRxxxx),()()(2)(1)(及nTnRxxxx),(*2*1*.如果n个数列极限存在，且),2,1(lim*)(nixxikik，则称)(kx收敛于*x，且记为*)(limxxkk.例 10 向量

2、空间3R，设向量序列3)()21,)11(,1(RkkxTkkk，则由于k时，01k，ekk)11(，021k，因此有向量序列的极限*)()0,0(xexTk.例 11 设33R中，kkkkkkkkkkkx221/10/12/11/sin0/sin11)(.注意到1221lim211lim11limkkkkkkk和kkkkk1lim0sinlim，所以 Ixkk)(lim（单位矩阵）.3.5.2 向量的范数 【定义 2】（向量的内积）设Tnxxxx),(21，nTnRyyyy),(21（或nC）.将实数yxyxT),(niiiyx1（或复数niiiHyxyxyx1),(）称为向量x与y的内积（

3、或数量积）.非负实数21122)(),(niixxxx称为向量x的长度，即向量x的欧氏范数.将向量长度概念推广，得到向量范数定义.【定义 3】（向量的范数）如果对于任意的向量nRx（或nC），都有一个实值函数xxN)(与之对应且满足：（1）非负性：0 x，0 x当且仅当0 x；（2）齐次性：xx，R（或C）；（3）三角不等式：yxyx，对任意向量nRyx,（或nC）.则称xxN)(为nR（或nC）上的一个向量范数.下面是一些常用的向量范数.（1）向量的 1-范数：njjxx11；（2）向量的 2-范数：21122)(),(njjxxxx；（3）向量的-范数（最大范数）：jnjxx1max；（4

4、）向量的p-范数：pxxpnjpjp1,)(11.例 12 验证1x符合向量范数定义.解 1）非负性：由定义 3，显然01x，而01x，即),2,1(0njxj，也即0 x；2）齐次性：对于R，nRx，总有 1111xxxxnijnjj；3）三角不等式：111111yxyxyxyxnjjnjjnjjj.因此，由范数定义，1x是nR中的向量范数.由已知范数可以构造新的范数.例 13 设a是nR的一种范数，A是n阶非奇异矩阵，定义)(nabRxAxx，则b仍是nR的一种范数.证明 验证b符合范数定义.1）非负性：对非零向量nRx，则0Ax，从而0abAxx；2）齐次性：baabxkAxkkxAkx

5、)(，nRk；3）三角不等式：当nRyx,时，bbaaabyxAyAxyxAyx)(.因此，b是nR上的范数.一个向量的不同范数一般是不相等的.例如，nTRx)1,1,1(时，则有1,21xnxnx.这就给我们提出问题：范数是用来度量逼近程度的尺度，而范数的计算又不唯一，那么哪一种范数才能反映出真正的逼近性态呢？反映不同范数之间的联系定理如下.【定理 4】（范数的连续性）设x是nR中向量Tnxxxx),(21 的范数，则它是nxxx,21的n元连续函数.【定理 5】（范数的等价性）nR上定义的任何两种范数x与x是等价的，即存在正数210kk 使对一切nRx，成立不等式 xkxxka21.（3.

6、17）范数的等价关系（3.17）说明了：任何两种范数作为逼近度量的尺度，逼近性态是一样的.即，如果0 x（或），则0 x（或）.可以证得，常用的向量范数等价关系如下：xnxx1，xnxx2，1211xxxn.3.5.3 矩阵的范数 【定义 4】（矩阵范数）nnR是n阶方阵全体的集合，如果nnRA，A的某个非负的实值函数AAN)(满足：（1）非负性：0A，而00AA；（2）齐次性：RkAkkA,（或Ck）；（3）三角不等式：BABA，对任意的nnRBA,；（4）相容性：BABA，对任意的nnRBA,.则称AAN)(是nnR上的一个矩阵范数.Frobenius 范数：如果把nnR中的方阵理解为2n

7、维向量，则由向量 2-范数的定义，可以得到nnR中矩阵的一种范数 21112)(ninjijFaA （3.18）称为A的 Frobenius 范数（简称F范数）.F范数显然满足矩阵范数定义.由于在大多数与估计有关的问题中，矩阵和向量会同时参与讨论，所以希望引进一种矩阵的算子范数，它是和向量范数相联系而且和向量范数相容的，即对任何向量nRx及nnRA，总有 ppxAAx.（3.19）下面构造矩阵的算子范数.利用公式（3.19）ppxAAx，从而AxxApp.令pxxy，则y是单位向量，在闭球：1py内pAy是变量的连续函数（定理 4），因此，它能达到最大值A，即 AAypyp1max.【定义 5

8、】(矩阵的算子范数)是nR上的任一向量范数，则 AxAyxAxAnnnRxxRyyRxx110maxmaxmax （3.20）为nnR上的一个矩阵范数，称为矩阵的算子范数.可以验证公式（3.20）符合矩阵范数的定义.常用的算子范数是从属于向量),2,1(pp范数的算子范数：（1）列范数：njijjaA11max；（3.21）（2）2-范数：AATA2，（3.22）其中AAT是方阵AAT的最大特征值；（3）行范数：njijiaA1max.（3.23）证明（1）列范数：设1x=1，则njTjnjnjjjnjjjxaxaxaAx11211),(，njjniijniinjijninjjijxaxaxa

9、Ax1111111)(niijjnjjniijjaxa111max)max(.设kj 时niijjniikaa11max，则取Tx)0,0,1,0,0,0(0（第k个分量为 1）时，niniikijjaaAx110max.（2）2-范数：因为2122max AxAx，但是),(),(22AxAxAxAxAxT，而AAT显然是对称、正定的.设021n为AAT的特征值，而n,21为 对应的标准正交向量组.0 x为单位向量，则有 nncccx22110，且 1),(1200220niicxxx.由于 1121121100220),(),(niiniiiniiiiniiiTccccAxAxAx，当10

10、 x时，111111221),(),(AAAT，所以 AAxTAxA12122max.例 14 设4321A，求范数FAAAA,21.解 642,31max1A.743,21maxA.4772.53016941FA.由于20141410AAT，所以2211522115,22115maxAAT，因此4650.5221152A.向量和矩阵范数的 MATLAB 函数 在 MATLAB 中，norm()函数求向量与矩阵的范数，其命令格式为 norm(X,p).当 X 为向量或矩阵时，norm(X,p)表示向量或矩阵 X 的 p 范数，例如，norm(X,1)表示 X的 1 范数，norm(X,2)表示

11、 X 的 2 范数，norm(X,inf)表示 X 的范数.对于矩阵 X，norm(X,fro)表示 X 的 F 范数.对于向量 X，p 可以取任意的数值和 inf，但对于矩阵 X，p 只能取上述四种值.缺省情况表示 2 范数，即 norm(X)=norm(X,2).对于例 14 中矩阵 A，调用 MATLAB 的 norm()函数，输入 A=1,-2;-3,4 A1=norm(a,1),A2=norm(a),A_inf=norm(a,inf),A_fro=norm(a,fro)计算得到 A1=6 A2=5.4650 A_inf=7 A_fro=5.4772 又如，输入向量 x=1,-2,0,

12、4，则得 norm(x,1)=7,norm(x)=4.5826,norm(x,inf)=4.3.6 矩阵的条件数与直接法的误差分析 解线性方程组的直接法产生误差的主要原因：1）不同的算法及舍入误差的影响；2）方程组本身固有的问题（病态或良态）.前面我们分析了方程组直接法的不同算法，本节我们将分析方程组的状态并估计算法的误差，即原始数据扰动对解的影响.考虑n阶线性方程组bAx，其中A为非奇异矩阵.由于A（或b）的数值是测量得到的，或者是计算的结果，在第一种情况下A（或b）常带有某些观测误差，在后一种情况A（或b）又包含有舍入误差.因此我们处理的实际矩阵是AA（或bb），下面我们来研究数据A（或b

13、）的微小误差对解的影响.首先考虑一个例子.例 15 设方程组bAx，即00001.8800001.626221xx，它的精确解为T)1,1(.现在考虑系数矩阵和右端项的微小变化对方程组解的影响，即考察方程组 00002.8899999.526221xx，其解变为T)2,10(.扰动后方程组的解面目全非了，真所谓“差之毫厘，失之千里”，这种现象的出现是完全由方程组的性态决定的.利用 MATLAB 函数解线性方程组，输入 a=2 6;2 6.00001;b=8,8.00001;x=ab 得到解：x=1 1 输入 a=2 6;2 5.99999;b=8,8.00002;y=ab 得到解：y=10 -

14、2【定义 6】如果矩阵A或常数项b的微小变化，引起方程组bAx 解的巨大变化，则称此方程组为病态方程组，矩阵A称为病态矩阵（相对于方程组而言），否则称方程组为良态方程组，矩阵A称为良态矩阵.我们需要一种能刻划矩阵和方程组“病态”程度的标准.3.6.1 线性方程组的误差分析 设线性方程组为 bAx，（3.24）其中nnRA，nRbx,，且A非奇异.*x为精确解，x为解的误差，记xxx.设A为A的误差，b为b的误差.下面分别讨论x与A，b的关系.一、b 有误差而 A 无误差情形 将带有误差的右端项和带误差的解向量代入方程组，则 bbxxA)(.由于bAx，而得到bAx1，从而bAx1.另一方面，由

15、(3.24)式取范数，有xAb，即)0(1bbAx.因此可得解的相对误差估计式（3.25）.【定理 6】设A是非奇异矩阵，0 bAx，且bbxxA)(，则有误差估计式 bbAcondxx)(,（3.25）其中AAAcond1)(称为方阵 A 的条件数.【注】（1）解的相对误差是右端项b的相对误差的 cond(A)倍；（2）如果条件数越大，则解的相对误差就可能越大；（3）条件数成了刻划矩阵的病态程度和方程组解对A或b扰动的敏感程度.【定义 7】称条件数很大的矩阵为“病态”矩阵；称病态矩阵对应的方程组为病态方程组.反之，则称矩阵为良态矩阵，对应的方程组为良态方程组.二、A 及 b 都有误差的情形【

16、定理 7】设方程组0 bAx，A为非奇异矩阵，A及b都有误差，若A的误差A非常小，使11AA，则有误差估计式 bbAAAAAcondAcondxx)(1)(*.(3.26)证明 带有误差的方程组为 bbxxAA)(.由于bAx，代入上式整理得)()(1*11xAAxAAbAx.将上式两端取范数，利用向量范数的三角不等式及矩阵和向量范数的相容条件，则有 xAAxAAbAx1*11，整理可得)()1(*11xAbAxAA.若A足够小，使得11AA，则*111xAbAAAx.从而由bAx*1，有 bbAAAAAcondAcondxx)(1)(*.【注】仅A或b有误差是(3.26)式中0b或0A的特例

17、.例 16 已知方程组bAx 中 210603760471213,615141514131413121bA，若Tb)001.0,001.0,001.0(时，估计解的相对误差.解 由于,001.0,108.3333 1012132bb由式（3.25）有误差估计 0186.0101213001.020152xx，xx比右端项的相对误差6102308.9bb扩大了 2015 倍.例 17 设在例 15 方程组中，b有扰动b=(0，0.00001)T，试计算)(Acond，并说明b对解向量x的影响.解 由A求得1000001000003000005.3000001A，则 61108.400001.85

18、.600000)(AAAcond，%6006800001.0108.4)(6bbAcondxx.这说明右端项向量b其分量的万分之一的变化，可能引起解向量x有 600%的变化，如果我们事先不作分析，其解就难以置信了.因此，A是严重病态矩阵，相应的方程组是病态方程组.3.6.2 矩阵的条件数及其性质 常用的条件数有 AAAcond1)(，1111)(AAAcond，)()()(minmax2212AAAAAAAcondTT，分别称为矩阵A的-条件数、1-条件数和 2-条件数.条件数的性质（1）当TAA 时，)()()(minmax2AAAcond；（2）当A对称正定时，)()()(minmax2A

19、AAcond，其中，)(maxA和)(minA分别是矩阵A的按模最大和按模最小的特征值.（3）若1A存在，则有),0)()(),()(,1)(1RccAcondcAcondAcondAcondAcond；（4）若U为正交矩阵，即IUUT，则 1)(2Ucond，)()()(222AUcondUAcondAcond.判别一个矩阵是否病态是件极其重要的事情.MATLAB 提供了 cond(X,p)函数求矩阵 X 的p 条件数，例如，cond(X,1)：X 的 1 条件数；cond(X,2)：X 的 2 条件数；cond(X,inf)：X 的条件数；cond(X,fro)：X 的 Frobenius

20、 条件数.p 缺省情况表示 2 条件数，即 cond(X)=cond(X,2).例 18 下列希尔伯特（Hilbert）矩阵是一族著名的病态矩阵.)12/(1)1/(1/1)1/(13/12/1/12/11nnnnnHn.用 MATLAB 函数计算条件数nH.输入 for n=3:8 cond(hilb(n)end 得到 3 至 8 阶希尔伯特矩阵的条件数分别为：524.0568 1.5514e+004 4.7661e+005 1.4951e+007 4.7537e+008 1.5258e+010 由此可见，随着n的增加，nH的病态可能越严重.nH常出现在数据拟合和函数逼近中.例 19 用 M

21、ATLAB 求解线性方程组bxHn，其中eHbn，Te)1,1,1(.显然，这个方程组的精确解为Tex)1,1,1(*.下面用 MATLAB 求解，讨论n为不同值的情况.（1）当5n时，输入 n=5;H=hilb(n);b=H*ones(n,1);x=Hb;x 得到 x=1.0000 1.0000 1.0000 1.0000 1.0000 其解没有什么问题.（2）当10n时，输入 n=10;H=hilb(n);b=H*ones(n,1);x=Hb;x 得到 x=1.0000 1.0000 1.0000 1.0000 0.9999 1.0002 0.9996 1.0004 0.9998 1.00

22、00 其解有误差，但误差不大，可以接受.（3）当20n时，输入 n=20;H=hilb(n);b=H*ones(n,1);x=Hb;x 得到 x=1.0000 1.0000 0.9997 1.0039 0.9794 1.0126 1.3920 -1.1443 6.6138 -7.8690 11.9311 -15.2355 26.1941 -24.9072 16.5064 -10.4564 19.5516 -18.4902 10.8570 -0.9390 此时方程组的解已经面目全非了，基本上看不出解的各个分量为 1.方程组的解与精确值之间的误差如表 3-1 所示.表 3-1 Hilbert 系数

23、矩阵方程组的解的误差 n 5 10 20 60 cond(hilb(n)4.7661105 1.60251013 1.90841018 2.31911019 2*xx 2.673310-12 6.143110-4 106.1591 6.3902103 表 3-1 列出的误差大得超过我们的想像.原因在于 Hilbert 系数矩阵当20n时，条件数已达到 1018，尽管计算机计算的舍入误差很小，但由于巨大的条件数，还会产生很大的计算误差.对于病态方程组，通常的方法无法得到它的准确解，需要采用一些特殊的处理方法.3.6.3 病态方程组的处理 对于病态方程组，可采用高精度的算术运算，如双精度或扩充精度

24、，以改善或减轻方程组的病态程度.如果用无限精度运算即不存在舍入误差的话，即使条件数很大，也没有病态可言.我们也可对病态方程组作预处理，使改善方程组系数矩阵的条件数.例 20 设方程组bAx，即 666610/32000110/210/310/4100020003000921x，试验证其为病态方程组，且对其作预处理PbPAx，使)()(AcondPAcond.解（1）用 MATLAB 函数计算系数矩阵的条件数，得10109201.1)(Acond，显然方程组为病态方程组.（2）令)10,10,1(63 diagP，使PbPAx，其中 321,234123921PbPA，则有)(354.87)(A

25、condPAcond.显然，经过预处理后的方程组PbPAx 是良态的.奇异值分解法解病态方程组.奇异值分解（Singular-Value Decomposition）简称 SVD，对于n阶矩阵A，必存在正交矩阵U，V和对角阵S，使得A有奇异值分解 TUSVA.（3.27）有关奇异值分解的理论知识省略.在 MATLAB 中，函数 svd()表示矩阵的奇异值分解，其命令格式为 U,S,V=svd(A)其中，U，V为正交矩阵，S为对角阵.例如，求 4 阶 Hilbert 矩阵的奇异值分解，输入 U,S,V=svd(hilb(4)得到 U=-0.7926 0.5821 -0.1792 -0.0292

26、-0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146 S=1.5002 0 0 0 0 0.1691 0 0 0 0 0.0067 0 0 0 0 0.0001 V=-0.7926 0.5821 -0.1792 -0.0292 -0.4519 -0.3705 0.7419 0.3287 -0.3224 -0.5096 -0.1002 -0.7914 -0.2522 -0.5140 -0.6383 0.5146 矩阵S对角线上的元素称为奇异值，由大到小排列.由公式（

27、3.27）得到TUVSA11，因此对于线性方程组bAx 的解，有 niiiTiTvbubUVSbAx111，（3.28）其中，i是矩阵S对角线上的元素，iu是正交矩阵U的第i列，iv是正交矩阵V的第i列.由式（3.28）可知，当i接近 0 时，会对计算解向量产生较大的误差.如何克服它的影响呢？由公式（3.27）得到 niTiiivuA1，当i接近 0 时，后面的项对A的贡献非常地小.取一正的阈值（不妨设为），当i时，其相应的项就不再计算了，即 iTiiivuA，由此得到方程组解的近似公式 iiiTivbux.(3.29)由于在式（3.29）中的i，不会产生较大的计算误差.奇异值分解求解线性方程

28、组的MATLAB程序：svd_equation.m function x=svd_equation(A,b)%奇异值分解方法求解线性方程组，其中，%A为方程组的系数矩阵；%b为方程组的右端项；%ep为精度；%x为方程组的解.ep=1e-10;n=length(A);x=zeros(n,1);U,S,V=svd(A);sigma=diag(S);for i=1:n if abs(sigma(i)=ep x=x+(U(:,i)*b)/sigma(i)*V(:,i);end end 用函数svd_equation（）求解Hilbert系数矩阵的方程组，会得到较好的计算结果.对于例19，当20n时，输

29、入 n=20;H=hilb(n);b=H*ones(n,1);x=svd_equation(H,b);x 得到 x=1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 显然，得到的解向量没有问题.评注 对于良态问题，高斯消去法也可能给出很坏的结果，即说明这个算法是不稳定的，在高斯消去法中引进选主元的技巧，就得到了解方程组的完全主元素消去法和列主元素消去法，选主元素技巧的

30、根本作用是为了对舍入误差的增长加以控制，由此，在非病态情况下，完全选主元素消去法及列主元消去法是数值稳定的算法，这两种方法都是计算机上解线性方程组 的有效方法，但通常用列主元消去法即可.列主元消去法所用的 CPU 时间比不选主元的高斯消去法略多一点，完全主元素消去法所花费的 CPU 时间是列主元消去法的一倍，如果要求计算精度较高时，则选用完全主元素消去法.从代数上看，直接分解法和高斯消去法本质上一样，但如果我们采用“双精度累加”计算iiba，那么直接三解分解法的精度要比高斯消去法高.对于对称正定方程组，采用不选主元素的平方根法（或改进的平方根法）求解适宜.理论分析指出，在非病态情况下，解对称正定方程组的平方根法是一个稳定的算法，在工程计算中使用比较广泛.追赶法是解对角元占优势的三对角线方程组的有效方法，它具有计算量少、方法简单、算法稳定等优点.对于病态问题，最好扩大运算字长，例如采用双精度或扩充精度，通常，用双精度就能比较好地计算很多病态问题.也可对病态方程组作预处理，改善方程组系数矩阵的条件数.奇异值分解法也会得到较好的计算结果.