线性方程组AX=B的数值解法(j)[优制课件]

1、第3章 线性方程组 AX=B的数值解法 1相关知识 引言 n在自然科学和工程技术中很多问题的解决常常归 结为解线性代数方程组。例如电学中的网络问题， 船体数学放样中建立三次样条函数问题，用最小 二乘法求实验数据的曲线拟合问题，解非线性方 程组问题，用差分法或者有限元法解常微分方程， 偏微分方程边值问题等都导致求解线性方程组， 而且后面几种情况常常归结为求解大型线性方程 组。 n线性代数方面的计算方法就是研究求解线性方程 组的一些数值解法与研究计算矩阵的特征值及特 征向量的数值方法。 2相关知识 线性方程组求解问题 n考虑线性方程组 Ax = b n其中A是一个是一个(n n)的非奇异矩阵, x

2、是要求解 的n维未知向量, b是n维常向量 nnnnnn n n b b b x x x aaa aaa aaa 2 1 2 1 21 22221 11211 nnnnnn nn nn b x a x a xa b x a x a xa b x a x a xa 2211 22222121 11212111 3相关知识 线性方程组的解的存在性和唯一性 n定理3.4 设A是NN方阵，下列命题等价： 给定任意N1矩阵B，线性方程组AX=B有唯 一解 矩阵A是非奇异的（即A-1存在） 方程组AX=0有唯一解X=0 det(A) 0 4相关知识 线性方程组的解 n最常见的求线性方程组Ax=b的解的方法

3、是在方 程组两侧同乘以矩阵A的逆 nGram法则： bAAxA 11 bAx 1 Ax = b 1,2,., Ddet(A)0det() i i iii D xin D DAAA ib ，其中 ，是的第 列用 代替所得。 5相关知识 线性方程组的解（续1） n求逆运算和行列式计算由于运算量大，实 际求解过程中基本不使用，仅作为理论上 的定性讨论 n克莱姆法则在理论上有着重大意义，但在 实际应用中存在很大的困难，在线性代数 中，为解决这一困难给出了高斯消元法 n还有三角分解法和迭代求解法 6相关知识 解法分类 n关于线性方程组的数值解法一般有两类 直接法：若在计算过程中没有舍入误差，经 过有限步

4、算术运算，可求得方程组的精确解 的方法 迭代法：用某种极限过程去逐步逼近线性方 程组精确解的方法 迭代法具有占存储单元少，程序设计简单， 原始系数矩阵在迭代过程中不变等优点，但 存在收敛性及收敛速度等问题 7相关知识 3.3 上三角线性方程组 n定义3.2 NN矩阵A=aij中的元素满足对所 有ij，有aij=0，则称矩阵A为上三角矩阵； 如果A中的元素满足对所有ij，有aij=0，则 称矩阵A为下三角矩阵。 n定理3.5（回代）设AX=B是上三角线性方程 组，如果akk0，其中k=1,2,N，则该方程 组存在唯一解。 8相关知识 3.3 上三角线性方程组（续1） n定理3.6 如果NN矩阵A

5、=aij是上三角矩阵或下 三角矩阵，则 1122 1 det( ) N NNii i Aa aaa n 条件akk0很重要，因为回代算法中包含对akk的除 法。如果条件不满足，则可能无解或有无穷解 n 联系定理3.4，可知要条件akk0成立才能保证方 程组存在唯一解 9相关知识 3.3 上三角线性方程组（续2） n求解上三角线性方程组的回代算法 . . n n n a b x 11 2 11 11 13132121 1 )( )( a xab a x a x a xab x n k kk nn . 1 . , 1 1 . 1, 1 n n nn n nn b x + axa a x - ab

6、x nn n nnn n . 1, 1 . , 1 . 1 1 最后最后 . 33333 22323222 11313212111 nnnn nn nn nn b xa b xa xa b xa xa xa b x a x a x a xa 10相关知识 上三角线性方程组的求解 n基本算法： 12,.1 1 ,ni )/uxu(bx /ubx ii n ij jijii nnnn 11相关知识 上三角线性方程组的求解（续1） nn n n u uu uuu U Ubx 222 11211 (2)其中式可简写成， 12相关知识 3.4 高斯消去法和选主元 n求解有N个方程和N个未知数的一般方程组

7、 AX=B的一般做法：构造一个等价的上三角 方程组UX=Y，并利用回代法求解 n如果两个NN线性方程组的解相同，则称 二者等价 n对一个给定方程组进行初等变换，不会改 变它的解 13相关知识 3.4 高斯消去法和选主元（续1） n考虑一个简单的例子： n求解第二个方程，得 14 723 21 21 x x x x n第二个方程减去第 一个方程除以3再乘 以4得到的新方程， 得到新的方程组： 3 25 3 5 723 2 21 x x x 1 72 5 1 3 x 2 5x n 回代到第一个方程，得 14相关知识 3.4 高斯消去法和选主元（续2） n考虑包含n个未知数的方程组 or nnnnn

8、nn nn nn nn b x a x a x a xa b x a x a x a xa b x a x a x a xa b x a x a x a xa 332211 33333232131 22323222121 11313212111 n作如下行变换之后方程组的解向量 x 不变 对调方程组的两行 用非零常数乘以方程组的某一行 将方程组的某一行乘以一个非零常数，再加到另一行上 n 通过对增广矩阵A|B进行如上的行变换求解 15相关知识 3.4 高斯消去法和选主元（续3） 1121/a a nnnnnnn nn nn nn b x a x a x a xa b x a x a x a x

9、a b x a x a x a xa b x a x a x a xa 332211 33333232131 22323222121 11313212111 11112213311 22223322 31132233333 112233 0 nn nn nn nnnnnnn ax ax ax ax b ax a x ax b ax ax ax ax b ax ax ax ax b 16相关知识 3.4 高斯消去法和选主元（续4） 1131/a a 111/a an nnnnnn nn nn nn b xa xa xa b xa xa xa b xa xa xa b x a x a x a xa

10、 3322 33333232 22323222 11313212111 nnnnnnn nn nn nn b x a x a x a xa b x a x a x a xa b x a x a x a xa b x a x a x a xa 332211 33333232131 22323222121 11313212111 17相关知识 3.4 高斯消去法和选主元（续5） . 33333 22323222 11313212111 nnnn nn nn nn b xa b xa xa b xa xa xa b x a x a x a xa nnnnnn nn nn nn b xa xa xa

11、b xa xa xa b xa xa xa b x a x a x a xa 3322 33333232 22323222 11313212111 222 / j aa 18相关知识 3.4 高斯消去法和选主元（续6） . 33333 22323222 11313212111 nnnn nn nn nn b xa b xa xa b xa xa xa b x a x a x a xa . . nn n n a b x 利用3.3节的回代法求解上述上三角方程组 19相关知识 3.4 高斯消去法和选主元（续7） 1633 1023 1642 321 321 321 xxx x xx x xx 16

12、33 145 2 1 1642 321 32 321 xxx -x x x xx 8 2 5 145 2 1 1642 32 32 321 xx -x x x xx 7826 2810 1642 3 32 321 x -x x x xx 消去过程 20相关知识 3.4 高斯消去法和选主元（续8） 3 26 78 3 x 23 28 10 28 10 3 2 x -x - 23 1 16(4) 2 1624 3 2 1 x x x 回代过程 21相关知识 3.4 高斯消去法和选主元（续9） n上述消去过程中，如果akk=0，则不能使用第k行 消除第k列的元素，而需要将第k行与对角线下的 某行进行

13、交换，以得到一个非零主元。如果不能 找到非零主元，则线性方程组的系数矩阵是奇异 的，因此线性方程组不存在唯一解 n选主元以避免 ，如果此主元非零，则不 换行；如果此主元为零，则寻找第p行下满足 的第1行，将此行与第p行互换，使新主元非零。 ( ) 0 p pp a ( ) 0 p kp a 平凡选主元策略平凡选主元策略 22相关知识 3.4 高斯消去法和选主元（续10） n选主元以减少误差：把元素中的最大绝对值移到 主对角线上 例3.17和3.18 n 偏序选主元策略 |akp|=max|app|,|app+1|,|aN-1p|,|aNp| n 按比例偏序选主元（平衡）策略 sr=max|ar

14、p|,|arp+1|,|arN| 其中r=p,p+1,N 1 1 | | max, kpppppNp kppN aaaa ssss 23相关知识 3.4 高斯消去法和选主元（续11） n病态问题：矩阵A中元素的微小变化引起解 的很大变化 47. 1 3 99. 048. 0 21 2 1 x x 1 1 2 1 x x 47. 1 3 99. 049. 0 21 2 1 x x 0 3 2 1 x x cond(A)=207.012 24相关知识 -2-101234 -0.5 0 0.5 1 1.5 2 2.5 图形解释 25相关知识 3.4 高斯消去法和选主元（续12） n一个线性方程组称

15、为是病态的，如果 其系数矩阵接近奇 异且它的行列式接 近0 n 矩阵条件数 cond(A)=|A|A-1| 26相关知识 3.5 三角分解法 nA=LU:下三角矩阵L的主对角线为1，上三角矩阵 U的对角线元素非零 定义3.4 如果非奇异矩阵A可表示为下三角矩阵L和上 三角矩阵U的乘积：A=LU，则A存在一个三角分解 1112131411121314 2122232421222324 3132333431323334 4142434441424344 1000 1000 1000 1000 aaaauuuu aaaamuuu aaaammuu aaaammmu A非奇异蕴含着对所有的k有ukk0

16、，k=1,2,3,4. 27相关知识 矩阵的LU分解 n是否所有的非奇异矩阵A都能作LU分解呢？ n一个例子： nN阶方阵A有唯一LU分解的充要条件是A的各阶 顺序主子式均不为零 01 10 A 10 10 bd ac 28相关知识 3.5 三角分解法（续1） n 利用前代/回代算法求解形如Lx=b或Ux=b的线性 方程组是容易的 n 如果对一个给定的矩阵A，能够找到一个下三角 矩阵L和一个上三角矩阵U，使A=LU n 则求解线性方程组Ax=b的问题可以分解成两个简 单的问题： p Ly=b p Ux=y n 易见：Ax=(LU)x=L(Ux)=Ly=b 29相关知识 3.5 三角分解法（续2

17、） n 假设已有矩阵A: n 对A作LU分解: n 检验分解结果: 23 12 A 10 0.51 L 23 00.5 U 10231*20*01*30*0.523 0.5100.50.5*21*00.5*3 1*0.512 LU 30相关知识 3.5 三角分解法（续3） n 构造一系列乘数矩阵M1, M2, M3, M4, MN-1使得: (MN-1M4M3M2M1)A是上三角矩阵，把它重新记成U. n对44矩阵A，M1可取: 21 11 1 31 11 41 11 1000 -100 -010 -001 a a Ma a a a 31相关知识 3.5 三角分解法（续4） nM2可取:nM3

18、可取： 1 32 2 1 22 1 42 1 22 1000 0100 0-10 0-01 M A M M A M A M A 3 21 43 21 33 1000 0100 0010 00-1 M M M A M M A 32相关知识 3.5 三角分解法（续5） n 则U=(M3M2M1)A是上三角形矩阵 n 每个M矩阵都是下三角形矩阵 n 如M2的逆为: n 注意到每个M矩阵的逆只是它自身下三角部分元素取相反数 n A = (M3M2M1)-1 U = (M1)-1 (M2)-1 (M3)-1 U n 定义L = (M1)-1 (M2)-1 (M3)-1，则L就是一个对角元素全为1的下三

19、角矩阵，因为所有的M矩阵的逆都是对角元素全为1的下三角矩阵 1 -132 2 1 22 1 42 1 22 1000 0100 010 001 M A M M A M A M A 33相关知识 3.5 三角分解法（续6） n计算复杂性：高斯消去法与三角分解法的三角化 过程是一样的，都需要 3 1 1 ()(1) 3 N p NN Np Np 32 1 1 23 ()() 6 N p NNN Np Np 次乘法和除法 次减法 求解LUX=B又需要N2次乘法和除法，以 及(N2-N)次减法 34相关知识 3.5 三角分解法（续7） n 每一个M矩阵中都需要计算1/A(i,i) n 当第i个对角元素

20、为0或者很接近0时就没法计算M， 这时A的直接LU分解就没法继续进行 n 可以将第i行与它下面的某一行互换，该行的第i列 元素非零 n 带选主元过程的LU分解 35相关知识 3.5 三角分解法（续8） 之前我们构造了一系列的M矩阵使得 是上三角矩阵 现在我们构造一系列的M矩阵和P矩阵使得 是上三角矩阵 (MN-1. M4 M3 M2 M1)A (MN-1 PN-1 . M4 P4 M3 P3 M2 P2 M1 P1)A 36相关知识 3.6 求解线性方程组的迭代法 n考虑线性方程组b Ax nnnnnn n n b b b x x x aaa aaa aaa 2 1 2 1 21 22221

21、11211 nnnnnn nn nn b x a x a xa b x a x a xa b x a x a xa 2211 22222121 11212111 37相关知识 3.6 求解线性方程组的迭代法（续1） n 高斯消去法高斯消去法 受限于舍入误差和病态性 n 迭代法迭代法 另一种求解线性方程组的方法 n 给出初始估计值，通过迭代得到更好的解的近 似值 n 迭代法对求解大型线性方程组非常有效 n Jacobi（雅可比）（雅可比）和Gauss-Seidel（高斯（高斯-赛德赛德 尔）尔）方法方法 38相关知识 3.6 求解线性方程组的迭代法（续2） n将方程组改写成每个方程的左边只有一个

22、未知 数的形式： nn nnnnnn n nn nn a x a x a xab x a x a xab x a x a xab x )( )( )( 112211 22 21212 2 11 12121 1 n给出初始估计值 和迭代规则 00 2 0 1 , n xxx 39相关知识 Jacobi迭代法 nn n nnnnn n nn nn a x a x a xab x a x a xab x a x a xab x )( )( )( 0 1 1 0 2 2 0 1 1 1 22 0 2 0 1 2121 2 11 0 1 0 2 1211 1 n初始估计值 00 2 0 1 , n xx

23、x n迭代一步后的结果： 40相关知识 Jacobi迭代法（续1） nn k n nn k n k nn k n k nn k k k nn k k a x a x a xab x a x a xab x a x a xab x )( )( )( 1 1 2 2 1 1 1 22 2 1 2121 2 11 1 2 1211 1 nk步迭代后的结果： 41相关知识 Jacobi迭代法（续2） n例： 52 83 12 32 321 21 x x x x x x x 2 5 2 )(5 3 8 3 )(8 2 1 2 )(1 221 3 31311 2 221 1 kk k kkkk k kk

24、k xx x xxxx x xx x nJacobi迭代公式： 42相关知识 Jacobi迭代法（续3） n初始迭代值 0, 0, 0 0 3 0 2 0 1 xxx 5 . 2 2 05 2 5 667. 2 3 008 3 8 5 . 0 2 01 2 1 0 21 3 0 3 0 11 2 0 21 1 x x xx x x x 1667. 1 2 833335. 15 2 5 2 3 )5 . 2(5 . 08 3 8 833335. 1 2 6667. 21 2 1 1 22 3 1 3 1 12 2 1 22 1 x x xx x x x 1, 3, 2 321 xxx n20步迭

25、代后 43相关知识 Jacobi迭代法（续4） n迭代会不会收敛到方程组的解？ n迭代到何时会终止？终止的判断条 件是什么？ 两个必须考虑的问题： 44相关知识 3.6 求解线性方程组的迭代法（续3） n定义3.6 设有NN维矩阵A，如果 1, N kkkj jj k aa 其中k1,2,N 则称A具有严格对角优势具有严格对角优势（严格对角占优严格对角占优）。 n 定理3.15（雅可比迭代）设矩阵A具有严格对角优 势，则AX=B有唯一解X=P。利用雅可比迭代可产生 一个向量序列Pk，而且对于任意初始向量P0，向 量序列都将收敛到P。 45相关知识 3.6 求解线性方程组的迭代法（续4） n向量

26、之间的距离可以用来判断Pk是否收敛到P。 因为两个向量P=(x1,x2,xN)和Q=(y1,y2,yN)之间 的欧几里德距离 1 2 2 2 1 () N jj j PQxy 计算复杂；而1范数具有度量的数学结构，也适 合作为一个一般化的“距离公式”。而且根据线性 代数的理论可知，如果两个向量的|1范数接近， 则它们的欧几里德范数|2也接近。所以定义两个 N维向量的距离为|1范数，用来确定N维空间中 的收敛性 46相关知识 3.6 求解线性方程组的迭代法（续5） n1-范数： 满足一般向量范数的性质 1 1 N j j Xx 定理3.16 设X和Y是N维向量，c是一个标量。则函数 |X|有如下

27、性质： p正定性：|X|0，|X|=0当且仅当X=0 p齐次性：|cX|=|c|X| p三角不等式：|X+Y|X|+|Y| 47相关知识 Gauss-Seidel迭代法 00 1 11221 1 11 10 1 221 12 2 22 111 1 1 12211 () () () nn nn nnnnnn n nn ba x a x x a ba x a x x a ba x a x ax x a n初始估计值 00 2 0 1 , n xxx n迭代一步后的结果： 48相关知识 1 112 21 1 11 1 1 221 12 2 22 111 1 1 12 211 () () () kk k n n kk k n n kkk k nnnnnn n nn ba x a x x a ba x a x x a ba x a x ax x a n 每一次迭代新产生的 被 认为是比 更好的xj的近似 值，所以在计算xj+1时用 来 替换 是合理的 Gauss-Seidel迭代法（续1） 1k j x k j x k j x nk步迭代后的结