第五章 求解线性方程组的迭代法.ppt

1、本章研究的对象是n阶线性代数方程组,在自然科学和工程技术中，许多问题的解决常常归结为线性代数方程组的求解。,例如：电学中的网络问题；船体放样中三次样条求解问题；求解非线性方程组问题；常微分方程及偏微分方程数值求解问题等，都可归结为求解线性代数方程组。,有关线性代数方程组解的存在唯一性及解的结构等理论问题，线性代数已作了详细讨论。,本书只介绍线性代数方程组的两类求解方法：直接法和迭代法。直接法是指经过有限步计算后求得方程组精确解的方法。,第五章 求解线性方程组的直接法,(1.1),若用矩阵和向量的记号来表示,若A是非奇异矩阵，则方程组有唯一解。,记D = detA，应用Cramer法则可得,Di

2、是用b代替中D第i列而得到的相应行列式,在实际中，用Cramer法则求解(1.1)并不可行,P98,因此，用Cramer法则求解 (1.1) 的乘除次数为,n 阶行列式共有 n！项，每项 n 个因子，所以计算一个 n 阶行列式共需 (n -1) n! 次乘法.,应用Cramer法则求解方程组需要计算 (n +1) 个行列式，另外还需n次除法。,当 n= 10 时， N = 359251210 ； 当 n = 20 时，N = 9.70731020，,P98,2.1 高斯消去法,2.1.1 顺序高斯消去法,高斯消去法是一个古老的方法，但实践证明它仍是目前计算机上一种常用有效的方法，其基本思想是通

3、过初等变换将方程组 (1.1) 转化为一个等价的三角形方程组,(1.3),这个过程称为消元，得到 (1.3) 后，就可逐个求出 这个过程称为回代.,从 (1.3) 中的最后一个方程直接可求出,代入倒数第二个方程便有,一般有,具有过程如下：,2.1.1.1 顺序高斯消去法的计算过程,P100,为了统一符号，把方程组（1.1）改写成下面形式,用矩阵表示即为,其中,则（1.4）化为,将（1.4）第1行乘 后加到第 i 行(i =2, 3, ,n),其中,用矩阵符号表示为,P101,其中,将(1.6)中第2个方程乘 后加到第 i 行(i =3, 4, ,n),类似地，若,则(1.6)化为,其中,用矩阵

4、符号表示为,其中,P102,若 上述过程可以继续进行下去，重复n-1次后，可得到等价的三角形方程组,用矩阵符号表示为,其中,P102,消元过程中元素的计算公式可以归纳为,从上面消元过程我们可以看到，要顺序高斯消去法能进行下去的充分必要条件是A矩阵的所有顺序主子式不为零。,P102,回代过程算法,P102,P102,2.1.2 顺序高斯消去法求解算法,上述算法中，为了不增加存储量，消元过程直接在系数矩阵A中进行。 的引用主要是为了控制溢出。,2.1.3 高斯消去法的计算量,(1) 消元过程的计算量,第1次消元，计算 mi1 ( i = 2, 3, , n ) 需要n-1次除法，,计算,（i ,

5、j = 2, 3, , n ）,需要 次乘法,计算,（i = 2, 3, , n ）,需要 次乘法。,因此第一次消元的乘除次数为,加减法次数为,P103,第2次消元等价于一个 阶方程组的一次消元，,第2次消元后乘除次数为,加减法次数为,因此消元过程的计算量为：,乘除法次数,加减法次数,P103,（2）回代过程的计算量,从而顺序高斯消去法的计算量为：,乘除法次数,加减法次数,P103,在上述消元过程中，用到了矩阵初等变换：即将某一行减去另一行的常数倍。在后面的列主元和全主元高斯消去法中我们还会用到矩阵的行交换与列交换。今对相应初等变换的矩阵表示说明如下：,2.1.4 顺序高斯消去法的矩阵解释(矩

6、阵表示),利用单位坐标向量 ei ( i = 1, 2, , n) ，,交换单位矩阵I 的 i , j两行所得矩阵记为 Pij,单位矩阵I可以表示成,(1)交换矩阵的两行、两列的变换,P104,用Pij 左乘矩阵A,其作用就是交换A的第 i 行和第 j行;用Pij 右乘矩阵A,其作用就是交换A的第 i 列和第 j列。,（2）矩阵的第 i ( i = 2, 3, , n ) 行减去 mi1乘第1行的变换,这种变换相当于用矩阵,左乘矩阵A。,高斯消去法的消元过程，从矩阵运算的角度来看，从（1.5）化为（1.7）相当于系数矩阵和左端项左乘一个初等变换矩阵M1 ，即,类似地，记,P105,则,P105

7、,从而,注意到,均为下三角矩阵且对角线元素为1，因此，,为对角线元素为1的下三角矩阵，,记,则,即系数矩阵A的所有顺序主子式不为零，,其中L为对角线元素为1的单位下三角矩阵，,只要,由顺序高斯消去法我们看到：,U为上三角矩阵,P111,定理1.1,如果n阶矩阵A的所有顺序主子式不为零，则A有的LU分解。(第一部分),则A一定可以分解成,1.1.2 列主元高斯消去法,主元高斯消去法根据主元选取范围的不同，又分为列主元消去法和全主元消去法.,P105,在顺序高斯消去法的消元过程中，若出现,则消元无法进行，,即使,就会导致其它元素量级的巨大增长和舍入误差的扩散，引起解的失真。,用它作为除数，,先看一

8、个例子：,主元高斯消去法是为了控制舍入误差而提出来的一种算法。,例1 方程组,该方程组精确解为,解 利用顺序高斯消去法取4位有效数进行运算得,将第1个方程的 1764倍加到第2个方程，我们得,P105,通过回代我们得,x 2 = 1.001 , x 1 = 10.00,与准确解相差很大。这是由于用0.003000作除数,将第1个方程的 0.0005760倍加到第2个方程，我们得,我们交换一下方程的位置,再利用顺序高斯消去法求解。,回代，我们得,x 1 = 10.00 , x 2 = 1.000,因此，为使这种不稳定现象发生的可能性减至最小，必须在每次消元前选择一个相对较大的系数作为主元。,1.

9、1.2.1 列主元高斯消去法计算过程,第一步,在系数矩阵 的第1列元素中选出绝对值最大的元素 称之为第1列的主元，即,P106,小主元是不稳定的根源,如果i11 , 则交换第1行和第i1行元素及相应的右端项，交换完后再进行消元得到：,（1.12）,第二步，在系数矩阵 第2列 以下 n-1 个元素中选主元 ，即,如果i22，则对调（1.12）中的第2行与第i2行元素及相应的右端项，交换完后再进行消元得：,一般地，第k步，求 I k 使,如果i k k，则交换第 k行与第I k行及相应的右端项，经消元后得：,列主元消去法除了每步需要按列选主元并可能进行矩阵行交换外，其消元过程与顺序高斯消去法的过程

10、一样，归纳以上过程有以下算法。,P106,2.2.2 列主元高斯消去算法,（1）输入A、b及 ,P106,（3）回代,（4）打印 （i=1,2,n）,2.2.3 列主元高斯消去法的计算量,列主元高斯消去法乘除及加减次数与顺序高斯消去法相同，增加了选主元时 n(n-1) 次比较以及交换方程次序所需时间。,1.1.2.4 列主元高斯消去法的矩阵解释,记,则说明第1个主元在第 i1 行，交换第1行与第 i1 行,P107,（相当于用 左乘矩阵 ）,然后再进行消元，用矩阵表示即为,同样道理，我们有,由此得到等价的三角形方程组为,这时,例 3 解方程组,(它的准确解为： ， ， ),解 第1步：首先选取

11、第1列绝对值最大的主元，对于本方程组，绝对值最大的主元为 -18，交换第1行和第2行 得,P107,将第1个方程的2/3倍加到第2个方程，第1个方程的1/18倍加到第3个方程得,P108,第2步：选取第2列（第2个元素以下）绝对值最大的主元，对于本方程组，绝对值最大的主元为 7/6，交换第2行和第3行消元后得,把上面方程进行回代，就可逐步解出,2.3 全主元高斯消去法,但找主元和交换行列次序要花费大量的机器时间，因此有时我们不用全主元高斯消去法而用列主元高斯消去法。,P108,全主元高斯消去法与列主元高斯消去法类似，所不同的是选取主元的范围不同。全主元高斯消去法在整个系数矩阵中找绝对值最大的元

12、素作为主元，这样可控制舍入误差的增长，将舍入误差控制在一个最小的范围。,2.3.1 全主元高斯消去法计算过程,若j11 , 则交换第1列和第j1列,（ 左乘 ），,第一步，求 i1 , j1 使得,并交换x1 与,交换完后再进行消元得到：,（1.13）,如果i11 , 则交换 第1行和第i1行及相应的右端项,( 左乘 ),P108,记,上式变为,一般地，第k步（ k = 1, 2, ，n-1），求 ik ，jk 使得,如果ikk , 则交换 第k行和第ik行及相应的右端项,若jkk , 则交换 第xk列和第jk列,交换完后再进行消元得到：,记,上式变为,归纳上述过程有以下算法:,2.3.2 全

13、主元高斯消去法算法,P109,2.3.3 全主元高斯消去法的计算量,全主元高斯消去法乘除及加减次数与顺序高斯消去法相同，增加了选主元时,比较以及交换方程次序所需的时间。,P110,2.3.4 全主元高斯消去法的矩阵解释,记,由上述过程我们知道，主元在第 i1 行，第 j1 列，交换第1行与第 i1 行，第1列与第 j1 列，然后再消元即可得到 A(2) ，用矩阵表示即为,P110,同理我们得,两边同乘,注意到,由于,由（1.14）知,（1.15）,通过 （1.17）可求得Z ，由（1.16）最终求得,2.2 LU分解法,而从消元过程我们可以看到从 A(1)转化到 A(n) 的过程实际是进行若干

14、次初等变换的过程，由顺序高斯消去法的矩阵解释我们得到,从而,注意到,P110,前面介绍过的顺序高斯消去法可以把方程组,转化为等价的一个上三角形方程组,均为下三角矩阵且对角线元素为1，因此，,为对角线元素为1的下三角矩阵，,由顺序高斯消去法我们看到：只要方程,的系数矩阵A的所有顺序主子式不为零，,其中L为对角线元素为1的单位下三角矩阵，,则A一定可以分解成,记,U为上三角矩阵,如果n阶矩阵A的所有顺序主子式不为零，,设,P111,定理 1 . 1,证明,则A有唯一的LU分解。 (第二部分),其中 L , L1 为单位下三角矩阵，,U , U1 为上三角矩阵，,因为 A 非奇异，,于是有,故 存在

15、，,上式左端为单位下三角矩阵，右端为上三角矩阵，从而上式两边均为单位矩阵，从而证得,如果能实现这种分解，则求解方程组,就相当简单。注意到 A = LU ，因此,则,用到上(下)三角阵的乘积、逆阵仍为 上(下)三角矩阵,U = U1 , L = L1 , A的三角形矩阵分解的唯一性获证。,记,i = 2，n,解,我们得,我们得,解,2.2.1 直接LU分解法,下面将详细讨论矩阵A的LU分解，设 A =，即,要两个矩阵相等，其对应元素应相等利用矩阵乘法，由矩阵的第行对应相等得,P111,由矩阵的第列对应相等得,这就求出了U矩阵的第1行和L矩阵的第1列元素。,一般地，设U矩阵的前k-1行和L矩阵的前

16、k-1列已经求出，则由,我们得,又由,我们得,综上所述，A的LU分解公式如下,k = 1, 2, , n,上述计算公式有如下特点：U的元素按逐行求，L的元素按逐列求；先求U的第k行元素，然后求L的第k列元素，U和L的元素一行一列交叉进行。,例4 P114,练习 P126 EX2,2.1.1 直接LU分解法算法,P113,每次将计算结果 和 仍存放在矩阵 A 的相应元素 和 所占的单元内，不必占用新的单元。 只要将上述算法中出现的 和 相应变换为 和 即可。,2.1.2 直接LU分解法的计算量（乘除次数）,LU分解计算U的计算量,计算L的计算量,总的计算量为,求解 的计算量,求解 的计算量,P1

17、14,例 4 用LU分解求解方程组,解 先把系数矩阵进行LU分解得,求解,得 y1 = 0, y2 = 3, y3 = 0.5,求解,解得 x1 =1 , x2 = -1 , x3 = 1,故 AX= b 的解为,P114,用直接LU分解法求解方程组所需要的计算量仍为 ， 和高斯消去法乘除法次数相同。但LU分解把对系数矩阵的计算和对右端项的计算分开了，这就使得求解系数矩阵相同而右端项不同的一系列方程组变得特别方便。,2.2.2 列主元LU分解法,假若第k-1步的分解已完成，在进行第 k 步分解时，为避免出现小的 作除数，计算,P115,选取行号 ,使,若 ,则对调 A 的第 k 行与第 行，再

18、求,于是完成了U的第 k 行和L的第 k列的计算,算法如下,P115,列主元LU分解算法,2.3 对称正定矩阵的平方根法和LDLT分解法,当A是对称正定矩阵时，存在一个实的非奇异下三角矩阵L1使,且当限定L1的对角线为正时，这种分解是唯一的。这种分解称为矩阵的Cholesky分解。,定理1.2 设A是对称正定矩阵，则A有如下分解：,其中，L是单位下三角阵，D为对角阵，且这种分解是唯一的。,P116,证明 因为A对称正定，从而A有唯一的LU分解,因为,故上式可化为,将右端三角矩阵分别记为 L , D 和R，因A对称有,A = L D R = R T D LT,因为A的LU分解唯一，故L = R

19、T , LT = R,A = L D LT,此分解显然是唯一的。,定理1. 3 n 阶对称正定矩阵A一定有Cholesky分解，当限定L1的对角线为正时，矩阵的Cholesky分解唯一。,P116,证明,使,又记 ，它是一个非奇异下三角阵，且对角线元素为正，故A有Cholesky分解,分解的唯一性显而易见。,下边推导Cholesky分解的计算公式,设 A = L LT , 即,其中 （ i,j = 1, 2, , n ）,第1步，由矩阵乘法有,且,故求得,(i = 2, 3, ， n ),一般地，设L矩阵的前 k1 列元素已求出，则,(i = 1, 2, , n ),第k步，由矩阵乘法得,P1

20、17,由于分解公式 ( 1.18 ) 中的每一步都有开方运算，故也称Cholesky方法为平方根法。,注意到,所以,这说明在分解过程中元素 的平方不会超过 A 的最大对角元，因而舍入误差的放大受到限制，所以平方根法求解对称正定方程组时可以不考虑选主元的问题。,可以证明，若用顺序高斯消去法求解对称正定方程组 AX=b，则有,( k =1, 2, , n ),其中 是第 k 步顺序高斯消去法过程所得到的元素，这说明高斯消去法求解对称正定方程组也可以不用选主元。从运算的角度看，平方根法是有利的，用平方根法 求解AX=b所需的乘除次数为,另外还有n次开平方运算。乘除次数只是高斯消去法的一半左右。,例

21、5 用平方根法解方程组,解：由平方根分解法可得,P117,求解,得,求解,得,平方根求解算法,P118,为避免平方根法的开方运算，利用定理1.2，也可以对A作 L D L T 分解，其分解方法可借助于前面所说的LU分解法。,从定理1.1我们可以看到，当A的所有顺序主子式不为零且对称时， A的LU分解法存在唯一，并且,U = D L T,其中， D = diag ( u11, u22 , , u nn) ,利用LU分解法先计算出U的第k行，由U的第k行我们得到D的第k个对角元素及L的第k列，,因此 L D L T分解的乘除法计算次数为,大概只是高斯消去法的一半，与平方根法大致相同，但避免了平开运

22、算。,例 6 用分解求方程组,解 由 分解法得,P119,解,得,解,得,最后解,得,这样就得到,LDLT分解算法,P120,P126 EX3,P126 EX3,25 向量与矩阵范数,用直接法求解线性方程组 AX=b 时，由于有舍入误差，只能得到近似解。为了进行解的误差分析和以后要介绍的迭代法的收敛性分析，我们首先介绍向量范数和矩阵范数的概念。,P120,定义1.1：设X R n， 表示定义在Rn上的一个实值函数， 称之为X的范数，它具有下列性质：,（3）三角不等式：即对任意两个向量X、Y R n，恒有,(1) 非负性：即对一切X R n，X 0, 0,(2) 齐次性：即对任何实数a R，X

23、R n，,向量范数是n维欧几里德空间中长度概念的推广。,P121,2.4.1 向量的范数,设X = (x1, x2, xn)T，则有,（4）更一般的p范数为,（1）向量的1范数,（2）向量的2范数,（3）向量的范数,容易验证，上面几种范数的定义满足范数的三个条件,常用的范数：,P121,定义1.2,在Rn上定义的任一向量范数 都与范数 等价，,即存在正数 M 与 m ( Mm ) 对一切XRn，不等式,成立。,定理1.4 Rn上定义的任何两个范数都是等价的。,证,对常用范数，容易验证下列不等式：,设给定Rn中的向量序列 ，即,其中,若对任何i (i = 1, 2, n )都有,则向量,称为向量

24、序列 的极限，或者说向量序列 依坐标收敛于向量，记为,定义2.3,定理1.5 向量序列Xk依坐标收敛于X*的充要条件是,向量序列依范数收敛与依坐标收敛是等价的。,P121,2.4.2 矩阵的范数,（1）当A = 0时， 0，当A 0时， 0,（2）对任意实数和任意A，有,（3）对任意两个n阶矩阵A、B有,（）对任意两个n阶矩阵A、B，有,定义,则称A为矩阵A的范数。,P122,设 ,若的非负实数 满足下列 四个条件,设 是 中矩阵序列，A是 中矩阵，,则,其中 和 分别表示 和 A 的分量。,定理 2.6,定义 1.4 对于给定的向量范数与矩阵范数， 如果 AR nn 和 X R n , 满足

25、,则称此矩阵范数与向量范数是相容的。,P122,并称其为由向量诱导出的矩阵范数。,现在考察由向量范数1 ，2 ， 产生的矩阵范数,P122,定理2.7,设A为n 阶方阵，Rn中已定义了向量范数 ，,则称 为矩阵A的范数或模，,记为,定理2.8 设n 阶方阵A = (aij)nn，则,（）与 相容的矩阵范数是,（）与 相容的矩阵范数是,其中1为矩阵ATA的最大特征值。,（）与 相容的矩阵范数是,（证明略）,这3种范数中，1范数与 范数的计算比较简单，2范数要计算矩阵的特征值，比较复杂。但2范数有一些好的性质，常用于理论分析.,P122,P123,此外，在 上的一个常用且易于计算的矩阵范数为,通常

26、称作Frobenius范数，,它是 上的2范数的自然推广。,A 的范数与A 的特征值之间的关系,矩阵A 的任一特征值绝对值不超过A的范数。,同向量范数一样，矩阵范数也有下面定理：,定理 1.6,设 是 中矩阵序列，A是 中矩阵，,则,其中 和 分别表示 和 A 的分量。,2. 5. 3谱半径,P123,定义1.5,矩阵A 的绝对值最大的特征值称为A的谱半径，,记为：,定理1.9,定理 1.10 设A是任意 阶矩阵，由A的各次幂所组成的矩阵序列,收敛于零,即 的充分必要条件是,P124,证明从略，读者可参阅冯康等编的数值计算方法,求解 时，A 和 b 的误差对解 x 有何影响？,设 A 精确b

27、有误差 ，得到的解为 ，即,绝对误差放大因子,相对误差放大因子,线性方程组的性态和解的误差分析,1.4.4 方程右端误差对解的影响,P124,P124,在这种情形，数 就是误差的放大率，,即解的相对误差不超过右端项相对误差 的倍。,我们称 为矩阵 A 的条件数，记作,条件数的性质：,）cond ( A )1,）cond ( kA )= cond ( A ) ,k 为非零常数,）若 ,则,2.4.5 系数矩阵误差对解的影响,P124,引理,而且,证,(1)设矩阵 I+F奇异,设矩阵F的范数小于1,即 则矩阵 I+F非奇异,矛盾,所以矩阵 I+F非奇异,(2),P125,设 b 精确，A有误差 ，

28、得到的解为 ，即,系数矩阵误差对解的影响,据前述引即可知，若,非奇异,且,*,所以,由于,(上页*),从而便得到不等式,P125,例：Hilbert (病态)阵,线性方程组,令,则准确解为,当n15时，各种方法都不能求到近似解,P124,cond (H2) =,27,cond (H3) ,748,cond (H6) =,2.9 106 ,cond (Hn) as n ,注：一般判断矩阵是否病态，并不计算A1，而由经验得出。 行列式很大或很小（如某些行、列近似相关）； 元素间相差大数量级，且无规则； 主元消去过程中出现小主元； 特征值相差大数量级。,*,近似解的误差估计及改善：,设 的近似解为

29、，则一般有,误差上限,cond (A),Step 1: 近似解,Step 2:,Step 3:,Step 4:,若 可被精确解出，则有就是精确解了。,经验表明：若 A 不是非常病态（例如： ），则如此迭代可达到机器精度；但若 A 病态，则此算法也不能改进。,改善方法：,EX2 用LU分解求解方程组,解 先把系数矩阵进行LU分解得,求解,得 y1 = 5, y2 = 24, y3 =33,求解,解得 x1 =109 , x2 = 27 , x3 = -66,P126,EX3 用平方根法解方程组,解：由平方根分解法可得,P126,求解,求解,P126,EX7 对矩阵,解：,P127 EX9,证,EX10,1. 上三角方程组的解法,直接法解方程组小结 （下标变量推导方法）,2. 下三角方程组的解法,3. 对称正定方程组的分解,当 时,列乘,U对角线元素,1. 矩阵范数,2. 矩阵条件数,3. 矩阵谱半经,矩阵几个重要参数,P127 EX8,几个例题,