线性代数方程组的数值解法ne.ppt

1、第3章 线性方程组的数值解法,3.1 引言与预备知识 3.2 高斯消去法 3.3矩阵三角分解法 3.4向量和矩阵的范数误差分析 3.5迭代方法,这一章讨论线性方程组,3.1 引言与预备知识,的数值解法.,在自然科学和工程技术中,很多问题归结为解线性方程组.例如电学中的网络问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性方程组，而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵(例如，阶数不超过150). 另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多).,3.1.1 引

2、言,有的问题的数学模型中虽不直接表现为含线性方程组,但它的数值解法中将问题“离散化”或“线性化”为线性方程组.因此线性方程组的求解是数值分析课程中最基本的内容之一.,关于线性方程组的解法一般有两大类：,但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解. 本章将阐述这类算法中最基本的和具有代表性的算法就是高斯消元法,以及它的某些变形和应用.这类方法是解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组(例如，大型带状方程组)的有效方法.,1. 直接法,经过有限次的算术运算,可以求得方程组的精确解(假定计算过程没有舍入误差).如线性代数课程中提到的克莱姆算法就是一种直接法.但该法

3、对高阶方程组计算量太大,不是一种实用的算法.,就是用某种极限过程去逐步逼近方程组精确解的方法. 迭代法具有计算机的存储单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛条件和收敛速度问题.迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)的重要方法.,为了讨论线性方程组数值解法，需复习一些基本的矩阵代数知识.,2. 迭代法,3.1.2 向量和矩阵,基本概念：,用Rmn表示全部mn实矩阵的向量空间；,用Cmn表示全部mn复矩阵的向量空间.,(由数排成的矩阵表，称为m行n列矩阵).,其中aj为A的第j列的m维列向量. 同理,(m维列向量).,其中biT为

4、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) 矩阵的行列式,设ARnn，则A的行列式可按任一行(列)展开,其中Aij为aij的代数余子式，Aij=(-1)i+jMij，Mij为元素aij的余子式.,行列式性质：,定理1 设ARnn，则下述命题等价：,(1) 对任何bRn，方程组Ax=b有唯一解.,(

5、2) 齐次方程组Ax=0只有唯一解零解x=0.,(3) det(A)0.,(4) A-1存在.,3.2 高斯消去法,本节介绍高斯消去法(逐次消去法)及消去法和矩阵三角分解之间的关系. 虽然高斯消去法是一种古老的求解线性方程组的方法(早在公元前250年我国就掌握了解方程组的消去法)，但由它改进、变形得到的选主元素消去法、三角分解法仍然是目前计算机上常用的有效方法.我们在中学学过消去法,高斯消去法就是它的标准化的、适合在计算机上自动计算的一种方法.,3.2.1 高斯消去法,设有线性方程组,或写成矩阵形式,简记为Ax=b.,如: 上三角矩阵所对应的线性方程组,解为,当m=n时，对三角形方程组的解非常

6、容易，如,例如:,下三角矩阵所对应的线性方程组,计算量（乘除法的主要部分）都为 n2/2.,解为,因此，我们将一般的线性方程组化成等价的三角形方程组来求解.,首先举一个简单的例子来说明消去法的基本思想.,例1 用消去法解方程组,解 第1步，将方程(2.2)乘上-2加到方程(2.4)上去，消去(2.4)中的未知数x1，得到,第2步，将方程(2.3) 加到方程(2.5)上去，消去(2.5)中的未知数x2，得到与原方程组等价的三角形方程组,显然，方程组是(2.6)是容易求解的，解为,上述过程相当于对方程的增广阵做初等行变换,其中ri用表示矩阵的第i行.,由此看出，用消去法解方程组的基本思想是用逐次消

7、去未知数的方法把原方程组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个

8、(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个方程中的

9、未知数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

10、)的求解过程(2.11) 称为回代过程.,注意：设Ax=b，其中ARnn为非奇异矩阵，如果a11(1)=0，由于A为非奇异矩阵，所以A的第1列一定有元素不等于零，例如al10，于是可交换两行元素(即r1rl)，将al1 调到(1,1)位置，然后进行消元计算，这时A(2)右下角矩阵为n-1阶非奇异矩阵，继续这过程，高斯消去法照样可进行计算.,总结上述讨论即有,定理5 设Ax=b，其中ARnn.,(1) 如果akk(k)0 (k=1,2,n)，则可通过高斯消去法将Ax=b约化为等价的三角形方程组(2.10)，且计算公式为,(a) 消元计算(k=1,2, ,n-1),(b) 回代计算,(2) 如果A

11、为非奇异矩阵，则可通过高斯消去法(及交换两行的初等变换)将方程组Ax=b约化为等价的三角形方程组(2.10).,例2 解方程组,解：消元,回代得,习题:求解下列方程组,，,(1),(2),3.2.2 高斯主元素消元法,由高斯消去法知道, 在消元过程中可能有akk(k)0的情况，这时消去法将无法进行；即使主元素akk(k)0但很小时，用其作除数，会导致其它元素数量级的严重增长和舍入误差的扩散，最后也使得计算解不可靠.,先看一个例子,高斯消去法也称主元素消去法 (条件det A0),即当akk(k)=0 时，高斯消元法无法进行；或 | akk(k) |1时，带来舍入误差的扩散。（小主元）,解 (方

12、法1)高斯消元法(用4位有效数字),例3,(小主元产生了大误差),(方法2)列主元消元法，先交换行,这个例子告诉我们，在采用高斯消去法解方程组时，小主元可能产生麻烦，故应避免绝对值小的主元素akk(k)，对一般矩阵来说，最好每一步选取系数矩阵(或消元后的低阶矩阵)中绝对值最大的元素作为主元素，以使高斯消去法具有较好的数值稳定性，这就是全主元素消去法，在选主元时要花费较多机器时间，目前主要使用的是列主元素消去法,上例2,其中系数矩阵,3.2.3 矩阵三角分解,将上三角矩阵A(n)记为U，由(2.14)得到,其中,为单位下三角矩阵.,这就是说，高斯消去法实质上产生了一个将A分解为两个三角形矩阵相乘

13、的因式分解，于是我们得到如下重要定理，它在解方程组的直接法中起着重要作用.,定理 (矩阵的LU分解) 设A为n阶矩阵，如果A的顺序主子式Di0 (i=1,2,n-1)，则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，即 A=LU, 且这种分解是唯一的.,对于例1，系数矩阵,由高斯消去法 m21=0，m31=2，m32=-1，故,从而得到求矩阵行列式的计算公式为,消元法的应用,1. 求行列式,高斯消元法,2. 求逆矩阵(用高斯-若当消去法),例子 分别用高斯消元法和列主元消元法求矩阵的行列式.,解： 高 斯 消 元 法,大数“吃掉”了小数！,列 主 元 消 元 法,例4 用高斯-若当方

14、法求下列矩阵的逆矩阵.,解 由分块矩阵,所以得,3.3 矩阵三角分解法,高斯消去法有很多变形，有的是高斯消去法的改进，改写，有的是用于某一类特殊矩阵的高斯消去法的简化. 下面我们将介绍矩阵的直接三角分解法，解特殊方程组用的平方根法及追赶法.,定义 如果 L为单位下三角阵，U为上三角阵，则称A=LU为杜里特尔(Doolittle)分解；如果 L为下三角阵，U为单位上三角阵，则称A=LU为克劳特(Crout)分解.,3.3.1 直接三角分解法(LU分解),在3.2.3已经通过高斯消去法得到一个将A分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，A=LU，其中,并由定理7得到这种分解是唯一的.,

15、将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L,U元素的递推公式，而不需要任何中间步骤，这就是所谓直接三角分解法. 一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求两个三角形方程组., Ly=b，求y；, Ux=y，求x.,1. 不选主元的三角分解法,设A为非奇异矩阵，且有分解式A=LU，其中L为单位下三角矩阵，U为上三角矩阵，即,其中,a11 a12 a1k a1n u11 u12 u1k u1n 第1框 a21 a22 a2k a2n l21 u22 u2k u2n 第2框 ak1 ak2 akk akn lk1 lk2 ukk ukn 第k框 an1 an2 a

16、nk ann ln1 ln2 lnk unn 第n框,比较式 A=LU 两端的元素, 按下图所示顺序逐框进行，先求ukj，后求lik. 由第一框可得,得,假设前k -1框元素已求出，则由,有了矩阵 A 的LU分解计算公式 , 解线性方程组 Ax=b 就转化为依次解下三角方程组 Ly=b 与上三角方程组 Ux=y 其计算公式如下：, Ly=b, Ux=y,例5 求矩阵,解 用紧凑格式,的LU (Doolittle)分解.,所以,得行列式,习题:已知 b=(1,2,0)T,1.求A的LU分解，其中L为单位下三角矩阵，U为上三角矩阵 2.利用分解结果求Ax=b的解,解得:x=(3,1,-1),习题:

17、已知方程组,求A的LU分解，其中 L为单位下三角矩阵， U为上三角矩阵,3.3.2 平方根法,实际问题中Ax=b，系数矩阵A大多是对称正定矩阵，所谓平方根法，就是利用对称正定矩阵的三角分解而得到求解对称正定方程组的一种有效方法，目前在计算机上广泛应用平方根法解此类方程组.,对称矩阵的LDLT 分解法,当A为对称矩阵, 且其顺序主子式均不为0时, 由定理可知A有唯一的LU分解.,为了利用A的对称性，将U再分解为,其中D为对角矩阵，U0为上三角矩阵，于是,A=LU=LDU0. (4.6),又 A=AT=U0T(DLT),由分解的唯一性即得 U0T=L.,代入(4.6)得到对称矩阵A的分解式A=LD

18、LT. 总结上述讨论有,定理10(对称矩阵的三角分解定理) 设A为n阶对称矩阵，且A的各阶顺序主子式均不为零，则A可以唯一分解为 A=LDLT, 其中L 是单位下三角矩阵，D是对角矩阵.,如果A为对称正定矩阵，则A的分解式A=LDLT中D的对角元素di均为正数.,于是有,由定理6得到,其中L1=LD1/2为下三角矩阵.,定理11(对称正定矩阵的三角分解或乔雷斯基(Cholesky)分解) 如果A为n阶对称正定矩阵，则存在一个实的非奇异下三角矩阵L使A=LLT，当限定L的对角元素为正时，这种分解是唯一的.,下面我们用直接分解方法来确定计算L元素的递推公式，因为,其中lii0(i=1,2,n).

19、由矩阵乘法及ljk=0(当jk时),得,于是得到解对称正定方程组Ax=b的平方根法计算公式,对于 j=1,2,n,求解Ax=b，即求解两个三角形方程组,(1) Ly=b，求y；(2) LTx=y，求x.,由计算公式知道,所以有,例1.,用平方根法解对称正定方程组,解:,即,所以原方程组的解为,本例中出现了大量的根式运算,原因为,考虑改变分解方式,例题 用平方根法求解对称正定方程组,解 首先对A进行Cholesky分解,求解Ly=b，得 y1=2, y2=3.5, y3=1.,求解LTx=y，得 x1=1, x2=1, x3=1.,由公式(4.7)看出，用平方根法解对称正定方程组时，计算L的元素

20、ljj需要用到开方运算. 为了避免开方，我们下面用定理10的分解式 A=LDLT.,即,称为改进的平方根法.,求解Ly=b及DLTx=y的计算公式为,3.3.3 追赶法,在一些实际问题中，例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条插值函数，都会要求解系数矩阵为对角占优的三对角线方程组 Ax=f，即:,其中, 当|i-j|1时, aij=0, 且满足如下的对角占优条件:,(a) |b1|c1|0; (b) |bi|ai|+|ci|, ai,ci0, (i=2,3,n-1). (c) |bn|an|0.,我们利用矩阵的直接三角分解法来推导解三对角线方程组(4.12)的计算

21、公式. 由系数阵A的特点，可以将A分解为两个三角矩阵的乘积，即 A=LU, 其中取L下三角阵, 取U为单位上三角阵, 这样求解方程组Ax=f的方法称为追赶法. 设,其中 为待定系数，比较(4.13)两边即得,求解Ax=f等价于求解两个三角形方程组.,(1) Ly=f, 求y; (2) Ux=y, 求x.,从而得到解三对角线方程组的追赶法.,我们将计算系数 及 的过程称为追的过程, 将计算方程组的解 的过程称为赶的过程. 合起来就是追赶法.,追赶法公式实际上就是把高斯消去法用到求解三对角线方程组上去的结果. 这时由于A特别简单, 因此使得求解的计算公式非常简单,追赶法的计算量比较小的.,3.4方

22、程组的性态、条件数,在分析方程组的解的误差及下章中迭代法的收敛性时，常产生一个问题,即如何判断向量 x 的“大小”，对矩阵也有类似的问题. 本节介绍n维向量范数和nn矩阵的范数. 向量范数是三维欧氏空间中向量长度概念的推广，在数值分析中起着重要作用.,3.4.1 向量与矩阵的范数,首先将向量长度概念推广到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

23、时成立;,我们还可以用其他办法来度量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(向

24、量的范数) 如果向量xRn(或Cn)的某个实值函数N(x)=|x|，满足条件:,由(3)可推出不等式.,(4) | |x|-|y| | |x-y|.,下面给出几种常用的向量范数, 设x=(x1,x2,xn)T.,容易证明前三种范数是的p-范数特殊情况, 其中,向量的-范数(最大范数),向量的1-范数,向量的p-范数(0p+),向量的欧氏范数,例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 , 其中|为向量的任一种范

25、数.,下面我们将向量范数概念推广到矩阵上去. 可以将Rnn中的矩阵A=(aij)nn当作n2维向量, 则由向量的2-范数可以得到Rnn中矩阵的一种范数,称为A的Frobenius范数. |A|F显然满足正定性、齐次性及三角不等式.,下面我们给出矩阵范数的一般定义.,定义4(矩阵的范数) 如果矩阵ARnn的某个实值函数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| (相容性)；,

26、则称 N(A) 是Rnn上的一个矩阵范数(或模).,上面我们定义的F(A)=|A|F就是Rnn上的一个矩阵范数.,上述条件(即不等式)就称为矩阵范数与向量范数的相容性条件.,由于在大多数与估计有关的问题中, 矩阵和向量会同时参与讨论, 所以希望引进一种矩阵的范数, 它是和向量范数相联系而且和向量范数相容的, 即对任何向量xRn及矩阵ARnn都成立,|Ax|A| |x|.,为此我们再引进一种矩阵的范数.,定义5(矩阵的算子范数) 设xRn, ARnn, 给出一种向量范数|x|v(如v=1,2或), 相应地定义一个矩阵的非负函数,可验证|A|v满足定义4(见下面定理), 所以|A|v是Rnn上矩阵

27、的一个范数, 称为A的算子范数.,显然这种矩阵的范数|A|v依赖于向量范数|x|v的具体含义. 也就是说, 当给出一种具体的向量范数|x|v时, 相应地就得到了一种矩阵范数|A|v.,定理18 设xRn, A=(aij)Rnn, 则有常用的算子范数,（称为A的行范数）,（称为A的列范数）,（称为A的2-范数）.,其中max(ATA)表示ATA的最大特征值, 由于矩阵2-范数与ATA 的特征值有关，所以又称为A的谱范数.,例7 设,则,求相应的各种范数.,解,因为,令,得ATA的特征值,故,可见,定义6 设ARnn的特征值为i (i=1,2,n) , 称,为A的谱半径.,例子 设, 求A的谱半径

28、.,解,得A的特征值,故得A的谱半径为,定理19(特征值上界) 设ARnn, 则(A)|A|,即A的谱半径不超过A的任何一种算子范数(对|A|F亦成立).,定理20 如果ARnn为对称矩阵, 则|A|2=(A).,求A1, A2 , A , (A).,例子 设 ,解: 显然A1 =4,A=4,这里, 我们指出, 对于实对称矩阵A, 有(A)=|A|2.,习题:,= ,3.4.2 误差分析,5.6.1 矩阵的条件数,由于A(或b)元素是测量得到的, 或者是计算的结果, 在第一种情况A(或b)常带有某些观测误差, 在后一种情况A(或b)又包含有舍入误差. 因此我们处理的实际矩阵是A+A(或b+b)

29、, 下面我们来研究数据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的微小变化, 方程组的解却变化很大. 这样的方程组称为病态方程

30、组.,定义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

31、+b, Ax=b, x=A-1b,|x|A-1| |b|. (6.4),由(6.3)有,|b|A| |x|.,于是由(6.4)及(6.5), 得,定理22 设A是非奇异矩阵, Ax=b0 , 且,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

32、-1A)-1存在. 由(6.6)式,x= -(I+A-1A)-1A-1(A)x.,因此,设|A-1| |A|1, 即得,定理23 设A是非奇异矩阵, Ax=b0 , 且,(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|事实上刻画了解对原始数据变化的灵敏程度, 即刻画了

33、方程组的“病态”程度, 于是引进下述定义:,(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是好条件的).

34、 注意, 方程组病态性质是方程组本身的特性. A的条件数越大, 方程组的病态程度越严重, 也就越难用一般的计算方法求得比较准确的解.,例如对前面例8的矩阵作分析,由于条件数cond(A)1很大, 可见矩阵A的病态程度十分严重, 故由此方程组的解误差非常大.,计算其条件数,通常使用的条件数, 有,(1) cond(A) =|A-1|A| ;,(2) A的谱条件数;,当A为对称矩阵时,其中1, n为A的绝对值最大和绝对值最小的特征值.,习题,条件数的性质:,1. 对任何非奇异矩阵, 都有cond(A)v1. 事实上,cond(cA)v=cond(A)v.,2. 设A为非奇异矩阵且c0(常数), 则

35、,3. 如果A为正交矩阵, 则cond(A)2=1; 如果A为非奇异矩阵, R为正交矩阵, 则,cond(RA)2=cond(AR)2=cond(A)2.,由性质1知, 1cond(A), 即cond(A)总是大于等于1的数. 条件数反映了方程组的“病态程度”.条件数越小, 方程组的状态越好, 条件数很大时, 称方程组为病态方程组. 但多大的条件数才算病态则要视具体问题而定, 病态的说法只是相对而言.,条件数的计算是困难的, 这首先在于要算A-1, 而求A-1比解Ax=b的工作量还大, 当A确实病态时, A-1也求不准确；其次要求范数, 特别是求|A|2, |A-1|2又十分困难, 因此实际工

36、作中一般不先去判断方程组的病态.但是必须明白, 在解决实际问题的全过程中, 发现结果有问题, 同时数学模型中有线性方程组出现, 则方程组的病态可能是出问题的环节之一.,病态方程组无论选用什么方法去解,都不能根本解决原始误差的扩大,即使采用全主元消去法也不行.可以试用加大计算机字长,比如用双精度字长计算,或可使问题相对得到解决.如仍不行,则最好考虑修改数学模型,避开病态方程组.,如拟合问题中出现的正规方程组,就往往呈现病态,此时解决问题的方法之一是避开正规方程组,采用正交多项式拟合的方法,尽管后者比前者在理论上和实际计算中都复杂得多.,例如, n阶希尔伯特矩阵, 当n较大时, 是有名的病态矩阵.

37、 在函数逼近时, 有时得到方程组Hnx=b, 当n稍大, 用任何方法都难以解出理想的结果. 考查它的条件数表可见Hn在n稍大时, 即呈严重病态. 下面我们通过例题来考察Hn的病态情况.,例9 已知希尔伯特(Hilbert)矩阵,求条件数cond (H2)2, cond (H2)1和cond (H2).,解 因为,下面再分析H3的情况,(2) 考虑方程组,(1) 计算H3的条件数cond(H3),|H3|=11/6, |H3-1|=408, 所以cond(H3)=748.,同样可计算cond(H6)=2.9107, cond(H7)=9.85108.,由此看出当n越大时, Hn矩阵病态越严重.,

38、H3 x=(11/6, 13/12, 47/60)T=b,设H3及b有有微小误差(取3位有效数字)有,简记为(H3+H3)(x+x)=b+b. 方程组H3x=b与(6.8)的精确解分别为x=(1, 1, 1)T,x+x=(1.089512538, 0.48967062, 1.491002798)T.,于是 x=(0.0895, -0.5120, 0.4910)T,这就是说H3与b相对误差不超过0.3%, 而引起解的相对误差超过50%.,由上面的讨论, 要判别一个矩阵是否病态需要计算条件数cond(A) =|A-1|A|,而计算A-1是比较费劲的, 那么在实际计算中如何发现病态情况呢?,(1)

39、如果在A的三角约化时(尤其是用主元素消去法解(6.3)时)出现小主元, 对大多数矩阵来说, A是病态矩阵, 例如用选主元的直接三角分解法解方程组(6.8)(结果舍入为3位浮点数), 则有,(2) 系数矩阵A的行列式值相对说很小, 或系数矩阵某些行近似线性相关, 这时A可能病态.,(3) 系数矩阵A的元素间数量级相差很大, 并且无一定规则, 这时A可能病态.,用选主元素的消去法不能解决病态问题, 对于病态方程组可采用高精度的算术运算(采用双倍字长进行运算)或者采用预处理方法, 即将求解Ax=b转化为一等价方程组,选择非奇异矩阵P, Q使,cond(PAQ) cond(A).,一般选择P, Q为对

40、角阵或者三角矩阵.,当矩阵A的元素大小不均时, 对A的行(或列)引进适当的比例因子(使矩阵A的所有行或列按-范数大体上有相同的长度, 使A的系数均衡), 对A的条件数是有影响的. 这种方法不能保证A的条件一定得到改善.,其中A为非奇异矩阵, 当A为低阶稠密矩阵时, 第5章讨论的选主元消去法是有效的. 但对于大型稀疏矩阵方程组(A的阶数n很大，但零元素较多), 利用迭代法求解是合适的.,本章将介绍迭代法的一些基本理论及雅可比迭代法，高斯-赛德尔迭代法，超松弛迭代法，而超松弛迭代法应用很广泛。,下面举简例，以便了解迭代法的思想.,对线性方程组,Ax=b， (1.1),3.5 解线性方程组的迭代方法

41、,例1 求解方程组,记为Ax=b，其中,方程组的精确解是x*=(3,2,1)T. 现将改写为,或写为x=B0 x+f，其中,我们任取初始值，例如取x(0)=(0, 0, 0)T. 将这些值代入(1.3)式右边(若(1.3)式为等式即求得方程组的解，但一般不满足)，得到新的值x(1)=(x1(1), x2(1), x3(1)T=(3.5, 3, 3)T ，再将x(1)分量代入(1.3)式右边得到 x(2)，反复利用这个计算程序，得到一向量序列和一般的计算公式(迭代公式),简写为 x(k+1)=B0 x(k) +f，,其中k表示迭代次数(k=0,1,2,).,迭代到第10次有,从此例看出，由迭代法

42、产生的向量序列x(k)逐步逼近方程组的精确解是x*=(3,2,1)T. 即有,对于任何一个方程组x=Bx+f(由Ax=b变形得到的等价方程组)，由迭代法产生的向量序列x(k)是否一定逐步逼近方程组的解x*呢？回答是不一定. 请同学们考虑用迭代法解下述方程组,但 x(k)并不是所有的都收敛到解x*！,对于给定方程组x=Bx+f，设有唯一解x*，则,x*=Bx*+f .(1.5),又设x(0)为任取的初始向量, 按下述公式构造向量序列,x(k+1)=Bx(k)+f , k=0,1,2,.(1.6),其中k表示迭代次数.,定义1 (1)对于给定的方程组x=Bx+f , 用公式(1.6)逐步代入求近似

43、解的方法称为迭代法(或称为一阶定常迭代法，这里B与k无关). B称为迭代矩阵.,(2) 如果limx(k) (k)存在(记为x*), 称此迭代法收敛, 显然x*就是方程组的解, 否则称此迭代法发散.,3.5.1 迭代法的收敛性,6.3.1 一阶定常迭代法的基本定理,其中，A=(aij)Rnn为非奇异矩阵，记x*为(3.1)精确解，且设有等价的方程组,设线性方程组,Ax=b， (3.1),于是,设有解Ax=b的一阶定常迭代法,有意义的问题是：迭代矩阵B满足什么条件时，由迭代法产生的向量序列x(k)收敛到x*.,引进误差向量,由(3.3)式减(3.2)得到误差向量的递推公式,由6.1节可知，研究迭

44、代法(3.3)收敛性问题就是要研究迭代矩阵B满足什么条件时，有.,由上述讨论，需要研究x(k)的收敛性. 引进误差向量,由(1.6)减去(1.5)式，得(k+1)=B(k)(k=0,1,2,)，递推得,要考察x(k)的收敛性，就要研究B在什么条件下有lim(k)=0 (k)，亦即要研究B满足什么条件时有Bk0(零向量) (k) .,定义2 设有矩阵序列 Ak=(aij(k)Rnn 及A=(aij)Rnn，如果n2个数列极限存在且有,则Ak称收敛于A，记为lim (k).,例4 设有矩阵序列Ak, 其中Ak=Bk,而,且设|1，考查矩阵序列极限.,解 显然, 当|1时, 则有,矩阵序列极限概念可

45、以用矩阵算子范数来描述.,定理1 其中|为矩阵的任意一种算子范数.,定理2,定理3 设B=(bij)Rnn，则limBk=0 (k)(零矩阵)的充分必要条件是矩阵B的谱半径(B)1.,定理4(迭代法基本定理) 设有方程组,x=Bx+f. (3.4),及一阶定常迭代法,x(k+1)=Bx(k)+f. (3.5),对任意选择初始向量x(0)，迭代法(3.5)收敛的充要条件是矩阵B的谱半径(B)1.,定理4是一阶定常迭代法的基本理论.,定理3和定理4的结论和起来即为,(1) 迭代法x(k+1)=Bx(k)+f 收敛limBk=O；,(2) 迭代法x(k+1)=Bx(k)+f 收敛(B)1.,3.5.

46、1 基本迭代法,其中，A=(aij)Rnn为非奇异矩阵，下面研究任何建立Ax=b的各种迭代法.,设线性方程组,Ax=b， (2.1),其中，M为可选择的非奇异矩阵，且使Mx=d容易求解，一般选择A的某种近似，称M为分裂矩阵.,将A分裂为,A=M-N. (2.2),于是，求解Ax=b转化为求解Mx=Nx+b ，即求解,可构造一阶定常迭代法,其中 B=M-1N=M-1(M-A)=I-M-1A ， f=M-1b. 称 B=I-M-1A为迭代法的迭代矩阵，选取M矩阵，就得到解Ax=b的各种迭代法.,设aii0 (i=1,2,n)，并将A写成三部分,即A=D-L-U,3.5.2 雅可比(Jacobi)迭代法,设aii0 (i=1,2,n)，选取M为A的对角元素部分,即选取M=D(对角阵)，A=D-N，由(2.3)式得到解方程组Ax=b的雅可比(Jacobi)迭代法. 又称简单迭代法.,其中B=I-D-1A=D-1(L+U)=J, f=D-1b. 称J为解Ax=b的雅可比迭代法的迭代矩阵.,于是雅可比迭代法可写为矩阵形式,其Jacobi迭代矩阵为,下面给出雅可比迭代