计算方法线性方程组解法.ppt

《计算方法线性方程组解法.ppt》由会员分享，可在线阅读，更多相关《计算方法线性方程组解法.ppt（122页珍藏版）》请在装配图网上搜索。

1、计算方法 (力学系本科生 ) 第三章 线性方程组解法 (Solution for Linear Algebraic Equations ) 3.1 问题的提出 第三章 线性方程组解法 n阶线性方程组 3.1 问题的提出 1 1 1 1 2 2 1 3 3 1 1 2 1 1 2 2 2 2 3 3 2 2 1 1 2 2 3 3 . . . . nn nn n n n n n n n a x a x a x a x b a x a x a x a x b a x a x a x a x b 3.1 问题的提出 线性方程组 Ax=b，其中 A是 n维 方阵， x是 n维未知数向量， b是 n维常

2、数向量。 1 1 1 2 1 1 1 2 1 2 2 2 2 2 12 . . , . . . . . . . . . . . . . n n n n n n n n a a a b x a a a b x a a a b x A b x 3.1 问题的提出 如果 A是非奇异阵时 ，方程组有唯一解， 且可以用 克莱姆 (Grammer)法则表示 : , ( 1 , 2 , ., )ii Dx i n D 其中 xi是解向量 x*的第 i个分量， D=detA, Di是用 b代替 A的第 i列后得到矩阵的行列 式。 3.1 问题的提出 克莱姆方法求解计算量太大，需要计 算 (n+1)个 n阶行列

3、式，共需要 (n+1)!次乘 法运算。 3.1 问题的提出 求解线性方程组的数值方法有两大类： 1)直接法 (direct methods)。 经过有限次 算术运算可求方程组 精确解 的方法 (实 际上，由于舍入误差不可避免，一般 得不到精确解 )。适合于求解低阶稠密 阵方程组。 3.1 问题的提出 2) 迭代法 (iterative methods)。采用极限过 程去逐步逼近线性方程组精确解的方 法。迭代法需要计算机存储单元较少， 对计算机要求不高，程序设计简单， 但有收敛性和收敛速度方面的问题。 迭代法是求解大型稀疏矩阵方程的重 要方法。 3.1 问题的提出 我们在本章将要学习迭代法有：

4、雅可比 (Jacobi)迭代法 高斯塞德尔 (Gauss-Seidel)迭代法 超松弛迭代法 (Successive overrelaxation method, SOR)。 3.1 问题的提出 追赶法 (Forward elimination and backward substitution)。 我们在本章将要学习直接解法有： 高斯消去法 (Gauss Elimination)， 高斯主元素消去法 (Gauss Elemination with pivoting), 三角分解法 (LU decomposition)， 3.1 问题的提出 Jac ob i Ga us s Sei del 迭

5、 代 法 迭 代 法 迭 代 法 超 松 弛 迭 代 法 不 选 主 元 高 斯 消 去 法 选 主 元 （ 列 选 ， 全 选 ） 直 接 法 三 角 分 解 法 追 赶 法 【 历史注记 】 线性代数方程组数值解法有着 悠久的历史。我国古代数学著作 九章算术 (公元 1世纪 )的“方程”章中就有了较好的线 性方程组数值解法相当于现代对方程组 的增广矩阵进行初等变换、消去未知数的方 法。中世纪的印度数学家也可以求解线性方 程组。例如 12世纪的 婆什迦罗 的著作中，也 有求解线性方程组的内容。 3.1 问题的提出 在欧洲， 16世纪的 比特奥 在其 算术 （ 1559）中采用了与 九章算术

6、类似的消 元法。日本数学家 关孝和 在其 解伏题之法 一书 (1683)中首先采用了类似于现代行列式 法求解了三元线性方程组。稍后， 莱布尼茨 提出关于行列式解线性方程组的思想 (1693)。 1721年 马可劳林 用行列式展开式的方法给出 了二元、三元、四元线性方程组的解法， 3.1 问题的提出 但他的符号记法不完善。 1750年， 克莱姆 给 出了现在比较通用的线性方程组行列式解法， 即克莱姆法则。 1764年， 贝祖 用行列式建立 了线性方程组的一般理论。但由于当时计算 的效率很低，这一理论几乎只有理论的意义， 实际上只能求出未知数很少的线性代数方程 组的解。只是在 20世纪中叶电子计算

7、机问世 并投入应用之后，大型线性代数方程组的数 值求解才成为可能。 3.1 问题的提出 如何利用计算机更精确、更有效地求解大型 线性方程组，是计算数学中最重要的课题之 一。 现代计算实践中，常用的线性代数方程组 数值解法有 直接法 和 迭代法 两大类。直接法 是在没有舍入误差的假设下，经过有限次运 算就可得出方程组的精确解的方法，如各种 消元法。迭代法则采用逐次逼近的方法 ,即从 一个初始值出发，按照一定的计算格式 (迭代 公式 )，构造一个向量的无穷序列，其极限才 3.1 问题的提出 是方程组的精确解，用有限次运算得不到精 确解。迭代法是牛顿最先提出来的， 1940年 经 司威尔 提出的松弛

8、法也是一种迭代法，共 轭梯度法则是另一种迭代法，是 弗莱彻 等人 于 20世纪 60年代提出来的。 3.1 问题的提出 例 3.1 3.1 问题的提出 5 2 8 3 2 0 2 6 xy xy 精确解为 *2 , 1xy 0 . 4 1 . 6 0 . 1 5 1 . 3 xy yx 将方程写为 取 ( 0 ) ( 0 ) 0 xy ( 1 ) ( 0 ) ( 1 ) ( 0 ) 0 . 4 1 . 6 1 . 6 0 . 1 5 1 . 3 1 . 3 xy yx 3.1 问题的提出 重复以上过程得 ( 2 ) ( 1 ) ( 2 ) ( 1 ) 0 . 4 1 . 6 2 . 1 2 0

9、 . 1 5 1 . 3 1 . 0 6 xy yx k 0 1 2 3 4 5 6 x(k) 0 1.6 2.12 2.024 1.9928 1.99856 2.000432 y(k) 0 -1.3 -1.06 -0.982 -0.9964 -1.00108 -1.000216 3.1 问题的提出 如果把原方程写为 6 . 6 6 7 8 . 6 6 7 2 . 5 4 . 0 xy yx 构造 ( 1 ) ( ) ( 1 ) ( ) 6 . 6 6 7 8 . 6 6 7 2 . 5 4 . 0 kk kk xy yx k 0 1 2 3 x(k) 0 8.667 35.335 -109.

10、126 y(k) 0 4.0 -17.668 -84.358 5 2 8 3 2 0 2 6 xy xy 3.1 问题的提出 例 3.2 1 2 3 1 2 3 1 2 3 8 2 1 2 4 1 0 2 1 3 2 5 1 6 x x x x x x x x x 构造 ( 1 ) ( ) ( ) 1 2 3 ( 1 ) ( ) ( ) 2 1 3 ( 1 ) ( ) ( ) 3 1 2 0. 12 5 0. 25 1. 5 0. 4 0. 1 2. 1 0. 6 0. 4 3. 2 k k k k k k k k k x x x x x x x x x ( 0 ) ( 0 ) ( 0 )1

11、2 3 0 xxx 3.1 问题的提出 得 3.1 问题的提出 构造 ( 1 ) ( ) ( ) 1 2 3 ( 1 ) ( ) ( ) 2 1 3 ( 1 ) ( ) ( ) 3 1 2 2.5 0.25 5.25 ( 2) 1.5 2.5 8.0 ( 3 ) 4 0.5 6.0 ( 1 ) k k k k k k k k k x x x x x x x x x ( 0 ) ( 0 ) ( 0 )1 2 3 0 xxx 1 2 3 1 2 3 1 2 3 8 2 1 2 4 1 0 2 1 3 2 5 1 6 x x x x x x x x x 由原方程 3.1 问题的提出 得 3.1 问题

12、的提出 如果 A非奇异，则线性方程组 Ax=b有 唯一解 x*，将上式化为 x=Bx+f，给出初 始向量 x0，则有： ！ xk+1=Bxk+f, k=0,1,2 可以构成一向量序列 xk，若向量序列 xk收敛于 x*，则 x*=Bx*+f, 即 x*是方程 组的解 。 这种方法称为 迭代法 , B称为 迭代矩阵。 3.1 问题的提出 构造迭代法的中心问题是建立一个由 本次近似值计算下一次近似值的规则。 用迭代法求解线性方程组时要解决的问 题有： 1) 构造一种迭代格式，由 xk计算 xk+1 3) 证明向量序列 xk的 收敛性 2) 给出初始向量 x0 4) 如果序列收敛，证明是原方程组的解

13、 5) 给出估计误差和迭代停止判据。 3.1 问题的提出 定义：在 n维空间中给定一个向量序 列 ， ，如果对每一个分 量 ,当 时都有极限 xi, 即 , 则称向量序列 有极限 , 或称 收敛于 x。 kx 12( , , . )k k k k Tnx x xx k ix k li m k iik xx kx 1 2 2( , , . . . ) Tx x xx kx 3.2 雅可比迭代 (Jacobi iteration) 第五章 线性方程组解法 3.2 雅可比迭代 最简单的迭代方法是从第 i个方程解出 未知数 xi， i=1,2, n 1, 1 ()n i i ij j j j iii

14、x b a x a 于是雅可比迭代式为 1 1 1 ( ) , 1 , 2 , . . . , n kk i i i j j j j iii x b a x i n a 把系数矩阵分解为 A=U+D+L，其中 U 为由 A上三角部分构成的上三角阵， L为 由 A下三角元素构成的下三角阵， D为由 A对角线元素构成的对角阵。 3.2 雅可比迭代 显然，所有 aii, i=1,2, n不为 零 上式才 有意义，从线性代数知，对于任何系数 方阵非奇异的方程组，通过适当交换方 程的顺序总可以使所有方程的 0iia 3.2 雅可比迭代 12 13 1 23 2 0 . 0 . 0 . . . . 0 n

15、 n a a a aa U 21 3 1 3 2 12 0 0 0 . . . . . . 0 nn a aa aa L 11 22 33 . nn a a a a D 于是原方程组为 (U+D+L)x=b 3.2 雅可比迭代 上式两边左乘 D-1得 x= D-1 b- D-1(U+L)x=Bx+f 其中 B=-D-1(U+L), f=D-1b 于是有迭代格式 xk+1=Bxk+f 例 3.3 用 Jacobi迭代格式解下面方程组。 3.2 雅可比迭代 解： Jacobi迭代格式为 1 2 3 8 1 1 1 2 1 0 1 4 1 1 5 3 x x x 3.2 雅可比迭代 取初始向量 x(

16、0)=(0,0,0)T，各次迭代结果 ( 1 ) ( ) ( ) 1 2 3 ( 1 ) ( ) ( ) 2 1 3 ( 1 ) ( ) ( ) 3 1 2 ( 1 ) / 8 ( 4 2 ) / 10 ( 3 ) / ( 5 ) k k k k k k k k k x x x x x x x x x k 0 1 2 3 4 5 6 x1(k) 0.0000 0.1250 0.2500 0.2263 0.2235 0.2251 0.2250 x2(k) 0.0000 0.4000 0.3150 0.3005 0.3060 0.3058 0.3056 x3(k) 0.0000 -0.6000 -

17、0.4950 -0.4870 -0.4946 -0.4941 -0.4948 3.3 高斯塞德尔迭代 (Gauss-Seidel iteration ) 第五章 线性方程组解法 在雅可比迭代中 , 计算第 k+1次迭代近 似值时用的是上一次即第 k次的近似值， 从式 3.3 高斯塞德尔迭代 1 1, 1 ( ) , 1 , 2 , . . . , n kk i i i j j j j iii x b a x i n a 可以看出， 在计算第 i个 xik+1分量时，前 面 i-1个分量 x1k+1, x2k+1 xi-1k+1已经从上 式中计算出来了，于是很自然会想到如 果把它们代入用来计算

18、xik+1可能会改进 迭代 ,于是就得到 Gauss-Seidel迭代格式： 写成矩阵形式为： 3.3 高斯塞德尔迭代 1 11 11 1 ( ) , 1 , 2 . . . ,ink k k i i i j j i j j j j iii x b a x a x i n a 1 1 1()k k k x D b L x U x 或 1() kk D L x b U x 如果 (D+L)-1存在，则 1 1 1( ) ( )k k k x D L U x D L b B x f 其中 3.3 高斯塞德尔迭代 11( ) , ( ) B D L U f D L b 【 注记 】 通常高斯塞得尔方

19、法比雅可比方 法有更快的收敛速度，但不是总这样，对于 某些方程组，雅可比迭代收敛，而高斯塞 得尔方法发散。即，并不是任何时候高斯 塞得尔方法都比雅可比方法好。 例 3.4 用 Gauss-Seidel迭代格式解下面方 程组，精确到 3位有效数。 3.3 高斯塞德尔迭代 1 2 3 8 1 1 1 2 1 0 1 4 1 1 5 3 x x x 解： Gauss-Seidel迭代格式如下 3.3 高斯塞德尔迭代 取初始近似值 x0=(0,0,0)T，各次迭代结果 1 1 2 3 11 2 1 3 1 1 1 3 1 2 ( 1 ) / 8 ( 4 2 ) / 10 ( 3 ) / ( 5 ) k

20、 k k k k k k k k x x x x x x x x x k 0 1 2 3 4 xk1 0.0000 0.1250 0.2344 0.2245 0.2250 xk2 0.0000 0.3750 0.3031 0.3059 0.3056 xk3 0.0000 -0.5000 -0.4925 -0.4939 -0.4936 3.4 逐次超松弛迭代法 (SOR) (Successive Overrelaxation Method) 第五章 线性方程组解法 逐次超松弛迭代简称 SOR方法，是高 斯塞得尔法的一种加速方法。 3.4 逐次超松弛迭代法 高斯塞得尔法迭代格式得到 1 11 11

21、 1 ( ) , 1 , 2 , . . . ,ink k k i i i j j i j j j j iii x b a x a x i n a 把 xik+1改进为 xik与 的加权平均，即 1k ix 11 1 1 1 ( 1 ) ( ) , 1 , 2 , . . . , k k k i i i in k k k i i i j j i j j j j iii x x x x b a x a x i n a 3.4 逐次超松弛迭代法 上式中 时 , 就是高斯塞得尔方法， 为保证迭代过程收敛 ,要求 1 02 当 时叫 低阶松弛法 ； 当 时叫 超松弛法 。 1 21 SOR方法收敛时，

22、希望选择一个最佳 的使收敛速度最快。目前还没有最佳松 弛因子 的一般算法理论，实际大都 由计算经验通过试算确定 的近似值 opt opt 超松弛迭代式的矩阵形式 3.4 逐次超松弛迭代法 1)直接由分量形式公式写 由 1 11 1 () in k k k k i i i i j j i j j j j iii x x b a x a x a 有 1 1 1()k k k k k x x D b L x U x D x 所以 1 1 1( ) ( ) ( )kk x D L D D U x D L b (证明 ) 3.4 逐次超松弛迭代法 1 1 1()k k k k k x x D b L x

23、 U x D x 1 1 1( ) ( ) ( )kk x D L D D U x D L b 11 ()k k k k k D x D x b L x U x D x 11k k k k k D x Lx D x U x D x b 1( ) ( )kk D L x D U D x b 2)由高斯塞德尔公式推导。 3.4 逐次超松弛迭代法 高斯塞德尔迭代公式的矩阵形式是 11k k k D x b Lx U x 加权平均 11( 1 )k k k x x x 11( 1 ) ( )k k k k D x D x b L x U x 1 1 1( ) ( 1 ) ( )kk x D L D U

24、 x D L b (证明 ) 3.4 逐次超松弛迭代法 11( 1 ) ( )k k k k D x D x b L x U x 1 1 1( ) ( 1 ) ( )kk x D L D U x D L b 11 ( 1 )k k k k D x L x D x U x b 1( ) ( 1 ) kk D L x D U x b 3.5 迭代法的收敛性 (convergence) 第五章 线性方程组解法 1 ( ) m a x i in A 矩阵 的特征值 的 绝对值最大值称为矩阵 A的谱半径，即 3.5 迭代法的收敛性 nnR A , 1 , 2 , . . ,i in 定理 (迭代法基本定

25、理 )：设有方程组 x=Bx+f，对于任意初始向量 x0及任意 f， 迭代公式 xk+1=Bxk+f收敛的充要条件是 ( ) 1 B 3.5 迭代法的收敛性 定理 (迭代收敛的充分条件 )：设有迭 代式 xk+1=Bxk+f，如果 ，则对于 任意初始向量 x0，这个迭代过程收敛于 方程组 x=Bx+f的唯一解 x*，并且有事后 估计 1qB *1 1 1 k k k q x x x x 以及事前估计 * 1 0 1 k k q q x x x x 3.5 迭代法的收敛性 定义：如果对于方阵 ，有 nnA 1, , 1 , 2 , . . . , n i i i j j j i a a i n

26、则称方阵 对角占优 。 12 13 1 2 3 87 98 97 xx xx x x x 定理：如果方程组 Ax=b的系数阵 对角 占优 ，则方程组有唯一解且对任意初始 向量 x0雅可比迭代 和 高斯塞德尔迭代 都收敛于真解。？、 3.5 迭代法的收敛性 【 思考题 3.1】 如何 对方程组进行调整， 使用 Gauss-Seidel迭 代格式求解时收敛？ 定理：如果方程组 Ax=b的系数阵 对称 正定 ，则方程组有唯一解且对任意初始 向量 x0高斯塞德尔迭代 收敛于真解。 3.5 迭代法的收敛性 Jacobi迭代格式的收敛性 3.5 迭代法的收敛性 Jacobi迭代矩阵 1 () J D U

27、L 特征方程 0 IJ 1 ( ) 0 I D U L 1 ( ) 0 D D U L 3.5 迭代法的收敛性 由于 1 0 D 0 D U L 所以 注意到 A D U L 所以， 将 A的对角线元素乘以 后取行列 式，令其等于零，就是 Jacobi迭代矩阵的 特征方程 。 例 3.5 讨论 用 Jacobi迭代格 式解方程组的 收敛性。 3.5 迭代法的收敛性 1 2 3 8 1 1 1 2 1 0 1 4 1 1 5 3 x x x 解： Jacobi迭代矩阵的特征方程为 8 1 1 2 1 0 1 0 1 1 5 展开得 3.5 迭代法的收敛性 Jacobi迭代格式收敛。 34 0 0

28、 1 2 3 0 解得 1 2 , 30. 14 60 84 , 0. 07 30 41 0. 21 44 87 i 1 2 3( ) m a x , , 0 . 2 2 6 5 8 1 B Gauss-Seidel迭代格式的收敛性 3.5 迭代法的收敛性 Gauss-Seidel迭代矩阵 1() G D L U 特征方程 0 IG 1( ) 0 I D L U 1( ) ( ) 0 D L D L U 3.5 迭代法的收敛性 由于 1( ) 0DL ( ) 0 D L U 所以 注意到 A D U L 所以， 将 A的对角线以下元素乘以 后取 行列式，令其等于零，就是 Gauss-Seide

29、l 迭代矩阵的特征方程 。 3.5 迭代法的收敛性 【 思考题 3.2】 Gauss-Seidel迭代矩阵 1() G D L U 至少有 1个特征值为零，为什么？ 例 3.6 讨论用 Gauss-Seidel迭 代格式解方程组 的收敛性。 3.5 迭代法的收敛性 1 2 3 8 1 1 1 2 1 0 1 4 1 1 5 3 x x x 解： Gauss-Seidel迭代矩阵特征方程为 8 1 1 2 1 0 1 0 5 展开得 3.5 迭代法的收敛性 Gauss-Seidel迭代格式收敛。 2( 4 0 0 1 0 1 ) 0 解得 1 2 , 30 , ( 1 1 7 ) / 8 0 1

30、 2 3( ) m a x , , 0 . 0 6 4 0 3 8 1 G 3.5 迭代法的收敛性 作业： (1) 写出 下面方程组的 Jacobi迭代格 式和 Gauss-Seidel迭代格式的算法 (2) 讨论它们的收敛性。 1 2 3 2 1 1 1 1 1 1 1 1 1 2 1 x x x 3.6高斯消去法 (Gauss Elimination) 第三章 线性方程组解法 高斯消去法是最古老的数值方法之一， 现在仍然是一个很有用的方法，它在计 算机上容易实现。 3.6 高斯消去法 其基本思想是在各个方程之间进行乘 法和加减运算，逐步消去方程中的未知 数。它分为 消去 和 回代 两个过程

31、。 A x b 给定线性方程组 3.6 高斯消去法 1 1 1 2 1 1 1 2 1 2 2 2 2 2 12 . . . . . . . . . . . . . . . n n n n n n n n a a a x b a a a x b a a a x b 第一步消元，令 ， 用 乘第一个方程再加到第 i个方程上 作为第 i个新方程，消去 x1的项，变为 1 1 1 1/ , 2 , 3 , . . . ,iim a a i n 1im 逐次进行同样过程，最后，经过 n-1次消 元，得到： 3.6 高斯消去法 1 1 1 2 1 1 1 ( 1 ) ( 1 ) ( 1 ) 2 2 2

32、2 2 ( 1 ) ( 1 ) ( 1 ) 2 . 0 . . . . . . . . . . . . . . 0 . . n n n n n n n a a a x b a a x b a a x b 以上这些步骤叫 消元过程 。 3.6 高斯消去法 1 1 1 2 1 1 1 ( 1 ) ( 1 ) ( 1 ) 2 2 2 2 2 ( 1 ) ( 1 ) . 0 . . 0 . . . . . . . . . . 0 0 . . n n nn n n n n a a a x b a a x b a x b 然后，从第 n个方程开始，依次解出 xn, xn-1, x1 3.6 高斯消去法 (

33、 1 ) ( 1 )/nn n n nnx b a ( 1 ) ( 1 ) ( 1 ) 1 ( ) / n i i i i i i j j i i ji x b a x a 1 , . ., 1in 高斯消去法的计算量 3.6 高斯消去法 消去过程的计算量。第一步计算乘数 mi1,(i=2,3, n)需要 i-1次除法运算，计算 aij(2)(i,j=2,3, n)需要 (i-1)2次乘法运算以 及 (i-1)2次加减法运算。 利用求和公式 2 11 ( 1 ) / 2 , ( 1 ) ( 2 1 ) / 6 , 1 nn ii i n n i n n n n 得到 3.6 高斯消去法 第 k

34、步 加减法次数 乘法次数 除法次数 1 (n-1)2 (n-1)2 n-1 2 (n-2)2 (n-2)2 n-2 N-1 1 1 1 合计 n(n-1)(2n-1)/6 n(n-1)(2n-1)/6 n(n-1)/2 于是 消去过程 乘除法次数为 2( 1 ) / 3nn 消去过程 加减法次数为 ( 1 ) ( 2 1 ) / 6n n n 3.6 高斯消去法 计算 b(n-1)的计算量 ( 1 ) ( 2 ) . . 2 1 ( 1 ) / 2n n n n 乘除法次数为 加减法次数为 ( 1 ) / 2nn 3.6 高斯消去法 回代计算量 乘除法次数为 加减法次数为 ( 1 ) / 2n

35、n 1 2 . . . ( 2 ) ( 1 ) ( 1 ) / 2n n n n n 3.6 高斯消去法 总的计算量 33 2 3 3 3 n n nn 乘除法次数为 加减法次数为 335 3 2 6 3 n n n n 高斯消去法还有没有办法进行，为什么？ 3.6 高斯消去法 ( 1 ) 0 , 1 , 2 , . , ( 1 ) i iia i n 【 思考题 3.3】 如果遇到 作业：写出高斯消去法的 Fortran程序。 3.7高斯主元消去法 (Gauss Elimination with pivoting ) 第三章 线性方程组解法 3.7 高斯主元消去法 从高斯消去法我们已经看出，

36、为使高 斯消去法能顺利进行，必须在每一步消 去步都满足条件 ，但若 ( 1 ) 0i iia ( 1 ) ( 1 ) ,ii ii k ia a k i 相应的系数 ( 1 ) ( 1 )/,ii k i ik iim a a k i 计算可能引起很大的舍入误差。 为此，需要改进高斯消去法。有两种改 进方法： 列选主元， 完全选主元。 3.7 高斯主元消去法 列选主元 通过交换方程而使得 aki(i-1), k=i,i+1, n, 中绝对值最大的一个换到 (i, i)位置而成 为第 i步的消去主元。 完全选主元法就是在系数矩阵的子块 1 , 1 1 , , 1 , . . . . . k k

37、k n n k n n aa aa 中找出绝对值最大的元素，作为第 k+1次 消去过程的主元。 3.7 高斯主元消去法 假设 1 1 m a xp q i j k i n k j n aa 则 k+1行与 p行交换，第 k+1列与第 q列交 换，右端也同时交换， 在做列交换时， 要注意未知量也作交换 ，即把 xk+1与 xp交 换。 3.7 高斯主元消去法 列选主元算法： 1. for k=1,2,n -1 2. 找出满足 元素的行位置 p , m a xp k ikk i naa 3. if , error,无唯一解， stop error 4. if ( )换行 ,pka pk 5. fo

38、r j=k,k+1, n 6. ak,j与 ap,j交换 3.7 高斯主元消去法 7. bk与 bp交换 8. Endfor 9. Endif 消去计算： / , 1 , . . . ,i k i k k km a a i k n , , 1 , . ,ij ij ik k ja a m a i j k n , 1 , . . . ,i i i k kb b m b i k n 10. Endfor 3.7 高斯主元消去法 回代计算： /n n n nx b a 11. for i=n-1,n-2,2,1 1 ( ) / n i i i j j i i ji x b a x a 12. End

39、for 3.7 高斯主元消去法 3.8 三角分解法 (LU decomposition ) 第三章 线性方程组解法 高斯消去法 3.7 三角分解法 第一次消元过程为 1 1 1 1 1 1/ , 2 , 3 , . . . , , 0iim a a i n a 用 -mi1乘上面第一个方程再加到第 i个方 程，即可消去第 i个方程中的未知量 x1。 这个过程实际上是给系数矩阵 A左乘这 样一个下三角阵 L1： 作第二次消元相当于给系数矩阵 A(1)左乘 了一个下三角阵 L2： 21 1 31 1 1 0 0 . 0 1 0 . 0 0 1 0 0 . . . . . 0 0 . 1 n m m

40、 m L 3.7 三角分解法 作第 k次消元相当于系数矩阵 A(k-1)左乘一 个下三角阵 Lk 2 32 2 1 0 0 . 0 0 1 0 . 0 0 1 . 0 . . . . . 0 0 . 1 n m m L ( 1 ) ( 1 ) ( 1 ) 2 2 22 22/ , 2 , . , , 0iim a a i n a 3.7 三角分解法 | | 1 1, 2, , 1 0 0 . 0 1 0 . . . . . . . . . 1 . 0 . . 0 . 0 . . . . . . . . . . . k K k kk kk nk m m m L 3.7 三角分解法 于是第 n-1

41、次消元过程可表示为： () 1 2 2 1. n nn L L L L A A () 1 2 2 1. n nn L L L L b b 于是 1 1 1 1 ( ) 1 2 2 1. n nn A L L L L A 下三角阵乘积也是一个下三角阵。 记 1 1 1 1 2 1. n L L L L 3.7 三角分解法 L是单位下三角阵 (即对角元素全为 1的 下三角阵 )。 21 12 1 0 . . . 0 1 . . . 0 . . . . . . . . . . . . . . . 1 nn l ll L 3.7 三角分解法 ()nA L A 叫矩阵 A的 三角分解 ，或 LU分解 。

42、如果 L为单位下三角阵，则叫 Doolittle分解 ， 若 U为单位上三角阵，则叫 Crout分解 。 A L U 只要 A的各顺序主子式不为零，则 A可唯 一分解成一个 单位下三角阵 L与一个 上 三角阵 U的乘积。 3.7 三角分解法 设 Ax=b， A=LU，则 Ax=LUx=b 于是令 Ux=y，则 Ly=b 1 1 1 2 1 2 1 2 2 2 12 1 0 . . 0 . . 1 . . 0 0 . . , . . . . . . 0 . . . . . . . . . . 1 0 0 . . n n n n n n u u u l u u l l u LU 3.7 三角分解法

43、 这样原来方程能化为两个简单方程组 由 求 y 1 1 , 1 , 2 , . . . , i i i i j j j y b l y i n 然后由 求 x 1 ( ) / n i i i j j i i ji x y u x u 3.7 三角分解法 1 , 2 , . .,in 现在显式地给出 lij和 uij的表达式 1 1 1 2 1 1 1 1 2 1 2 1 2 2 2 2 1 2 2 2 1 2 1 2 . . . 1 . . . . . . 1 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44、 . . . 1 nn n n n n n n n n a a a u u u a a a l u u a a a l l u 从第一行可以看出， u1j=a1j, j=1,2, n 3.7 三角分解法 从第一列可以看出 ai1=li1u11, i=2,3, n, 即 li1=ai1/u11, i=2,3, n 。 从第二行可以看出， a2j=l21u1j+u2j , j=2,3, n，则 u2j=a2j-l21u1j, j=2,3, n 从第二列可以看出， ai2=li1u12+li2u22 , i=2,3, n，则 li2=(ai2-li1u12)/u22, i=2,3, n 3.7 三角

45、分解法 1 1 1 2 1 1 1 1 2 1 2 1 2 2 2 2 1 2 2 2 1 2 1 2 . . . 1 . . . . . . 1 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 nn n n n n n n n n a a a u u u a a a l u u a a a l l u 一般地，如果 U的前 k-1行， L的前 k-1 行已经求出，则第 k (k=2,3, n)步 1 1 k k j k m m j k j m a l u u (行 ) 3.7 三角分解法

46、1 1 , , 1 , . k k j k j k m m j m u a l u j k k n (列 ) 1 1 k ik im m k ik k k m a l u l u 1 1 ( ) / , 1 , . k ik ik im m k k k m l a l u u i k n 直到第 n步， A全部分解成 LU。 3.7 三角分解法 3.9 追赶法 ( Forward elimination and backward substitution ) 第三章 线性方程组解法 在很多实际问题中，方程组 Ax=b的系 数矩阵 A是一个稀疏矩阵。 3.6 追赶法 1 1 1 2 2 2 2

47、111 0 . 0 0 0 0 0 0 0 . 0 0 0 0 0 0 0 . . . . . . . . . 0 0 . 0 0 0 0 0 0 0 0 . . . . . . . . . . . . . 0 0 0 0 0 0 0 0 . 0 0 0 0 0 0 0 0 0 iii nnn n n n b c x a b c x a b c a b c a b x 1 2 . . . . n d d d 假设非零元素集中分布在对角线及其相 邻两个次对角线上，且系数阵为对角占 优阵，即有： 11 , 2 , 3 , . 1 i i i nn bc b a c i n ba 3.6 追赶法 把

48、系数矩阵三角分解有 11 22 2 11 0 . . . 0 1 0 0 0 0 . . . 0 1 0 0 , . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 1 0 0 0 0 nn n n qr qr p qr p q LU 3.6 追赶法 利用 LU分解公式，写出 1 1 1 1 1 1 , , 2 , 3 , . . . , kk k k k k k k k q b r c rc p q a q p c b k n 3.6 追赶法 得到 11 , kkq b r c 1 1 , , 2 , 3 , .,k

49、k k k k k k a p q b p c k n q 于是 ,L y d U x y 得到 11 1 , 2 , 3 , . . . ,k k k k yd y d p y k n 3.6 追赶法 1 / 1 ( ) , 1 , 2 , . , 1 n n n k k k k k x y q x y c x k n n q 3.6 追赶法 3.10 其它应用 第三章 线性方程组解法 计算行列式值 3.6 其他应用 A L U 1 1 2 2d e t ( ) d e t ( ) . . . nnu u uAU 求矩阵逆 1 A A I 1 1 1 2 1 2 1 2 2 2 12 . 1

50、 0 0 0 . 0 1 0 0 . . . . . . . . . . . . . . . . . 0 0 0 1 n n n n n n x x x x x x x x x A 分别由 1 2 0 . , 1 , 2 , . ., . 1 0 i i ni x x in i x A 解出 xij, ij=1,2, n 于是 1 ijx A 3.6 其他应用 3.11 误差分析 (Error analysis) 第三章 线性方程组解法 3.6 误差分析 定义： 设 为方程 组 Ax=b的精确解，近似解与精确解之间 的 误差向量 定义为 * * * * 12 . T nx x x x * 1

51、11 * 2 *22 * . . n nn e xx e xx e xx e x x 一种衡量近似解精确度的方法是把近似 解 代入方程组，看 残差向量 的大小。 *x r b A x 于是 1 1 * 1 x A b A r x A r *1 e x x A r 3.6 误差分析 由范数性质有 11e A r A r 而 * Ax b *b A x 1* A x b 上面两个式子相乘得 1 * () er AA bx 3.6 误差分析 这个式子左端计算的是近似解 x的相对误 差，而 是 相对残差 。 /rb 近似解的相对误差小于相对残差的一个 倍数： 。这个倍数对于近似解的误 差估计很重要，定

52、义为矩阵 A的 条件数 1AA 1()c o n d A A A 11 ( ) 1c o n d I A A A A A 3.6 误差分析 线性方程组 Ax=b的解 x*是由系数矩阵 A和向量 b决定的，那么如果 A和 b有微小 的变化 (扰动 )时，解会如何变化呢？ (1) 假设 b有微小变化 ，此时解有变 化 ，即 b *x *() A x x b b 注意到 * Ax b 3.6 误差分析 于是 *A x b 即 *1 x A b 两边取范数得 *1 x A b *b A x 利用 得 *1 1 * x A b A b AA bbx 3.6 误差分析 (2) 假设 A有微小变化 ，此时解

53、有变 化 ，即 A *x *( ) ( ) A A x x b 整理并注意到 得 * Ax b * * *() A x x A x 0 即 * 1 * *() x A A x x 两边取范数得 3.6 误差分析 即 * 1 * * x A A x x * 11 * x A A A A A Axx 于是从 (1),(2)两种情况可以看出，如果 b 或者 A有误差，则解的相对误差不超过 它们相对误差的 倍。 1AA 3.6 误差分析 Remark: 条件数刻画了线性方程组的解 对于初始数据 (b,A)误差的灵敏度，条件 数是方程组的固有属性，与求解该方程 组的方法无关。对于病态方程组，一般 不能得

54、到令人满意的结果。 3.6 误差分析 定义 为线性方程组的 条件数 1AA 当 时，方程“病态”。 ( ) 1c ond A 通常使用的条件数有 1( 1 ) ( )c o n d A A A 1 m a x 2 2 2 m in () ( 2 ) ( ) () T Tc o n d AAA A A AA 其中 是 ATA的绝对值最大、最小 的特征值。 m a x m in, 3.6 误差分析 当 A对称时 1() n c ond A 是 A的最大特征值， 是 A的最小特 征值。 1 n 3.6 误差分析 Remark: 为了计算方程组的条件数，需 要计算矩阵 A的逆阵，求逆阵的工作量 很大，而且， 如果 A是病态时，其逆阵 的计算也不准确 。 (2)系数阵某些行 (列 )近似线性相关 3.6 误差分析 判断方程组 病态 的参考： (1)用选主元消去法求解时出现 小主元 (3)系数阵元素量级相差很大 , 且无规则 本章小结 第三章 线性方程组解法 3.6 小结 迭代法： 雅克比迭代法 高斯塞德尔迭代法 超松弛迭代法 SOR。 追赶法。 直接解法： 高斯消去法 高斯主元素消去法 三角分解法 LU 其它应用 条件数 3.6 小结