线性方程组求解数值方法

1、线性方程组求解数值方法 Ch3 线性方程组求解的数 值方法 Numerical Solutions of Systems of Linear Equations 线性方程组求解数值方法 3.0 Introduction 线性方程组线性方程组 工程问题工程问题（电子网络，船体放样等）（电子网络，船体放样等）自然科学自然科学（实验数据的最小二乘法曲线拟合）（实验数据的最小二乘法曲线拟合） 解非线性方程组解非线性方程组解常解常/偏微分方程组偏微分方程组（差分法，有限元法）（差分法，有限元法） 线性方程组求解数值方法 3.0 Introduction 线性方程组的系数矩阵：线性方程组的系数矩阵： （1

2、）低价稠密矩阵（阶数）低价稠密矩阵（阶数 L,U,P=lu(A) 其中其中 nji a a l lll ll l L j jj j ij ij nnn , 1, 1 1 1 1 321 3231 21 【作业】【作业】P.67 题题1, 2 线性方程组求解数值方法 3.1.3 3.1.3 选主元消去法选主元消去法 /* Pivoting Strategies */ 例：单精度解方程组例：单精度解方程组 2 110 21 21 9 xx xx /* 精确解为精确解为 和和 */ .1000.00. 1 101 1 9 1 x 8个个 .8999.99. 02 12 xx 8个个 用用Gaussi

3、an Elimination计算：计算： 9 112121 10/ aam 999 2122 10101010.0 . 011 ma 8个个 9 212 1012 mb 99 9 10100 1110 0, 1 12 xx 小主元小主元 /* Small pivot element */ 可能导致计可能导致计 算失败。算失败。 3.1 Gaussian Elimination and LU Decomposition 线性方程组求解数值方法 全主元消去法全主元消去法 /* Complete Pivoting */ 每一步选绝对值最大的元素为主元素，保证每一步选绝对值最大的元素为主元素，保证 。

4、 1| ik m Step k: 选取选取;0|max| , ij njik ji aa kk If ik k then 交换第交换第 k 行与第行与第 ik 行行; If jk k then 交换第交换第 k 列与第列与第 jk 列列; 消元消元 列交换改变了列交换改变了 xi 的顺序，须记录交换次序，解完后再的顺序，须记录交换次序，解完后再 换回来。换回来。 列主元消去法列主元消去法 /* Partial Pivoting, or maximal column pivoting */ 省去换列的步骤，每次仅选一列中最大的元。省去换列的步骤，每次仅选一列中最大的元。 0|max| , ik

5、nik ki aa k 3.1 Gaussian Elimination and LU Decomposition 线性方程组求解数值方法 例：例： 211 1110 9 1,1 12 xx 110 211 1110 211 9 列主元法没有全主元法稳定。列主元法没有全主元法稳定。 0,1 12 xx 例：例： 211 10101 99 99 99 10100 10101 注意：这两个方程组注意：这两个方程组 在数学上严格等价。在数学上严格等价。 标度化列主元消去法标度化列主元消去法 /* Scaled Partial Pivoting */ 对每一行计算对每一行计算 。为省时间，。为省时间，

6、si 只在初始时计只在初始时计 算一次。以后每一步考虑子列算一次。以后每一步考虑子列 中中 最大的最大的 aik 为主元。为主元。 |max 1 ij nj i as nk kk a a . . . i ik s a 稳定性介于列主元法和全主元法之间。稳定性介于列主元法和全主元法之间。 3.1 Gaussian Elimination and LU Decomposition 线性方程组求解数值方法 实际应用中直接调用实际应用中直接调用Gauss Elimination 解解3 3阶线性方程阶线性方程 组的结果：组的结果： 结合全主元消去后的结合全主元消去后的 结果：结果： 3.1 Gauss

7、ian Elimination and LU Decomposition 线性方程组求解数值方法 3.2 Cholesky分解分解 【Matlab】 设设A是对称正定矩阵，则是对称正定矩阵，则A有以下形式的有以下形式的LU分解：分解： A=PTP （3.2.1） 其中其中P是一个上三角矩阵是一个上三角矩阵, 上式称为上式称为A的的Cholesky分解分解 (Choleskian decomposition). P=chol(A) 定理定理 【作业】【作业】P.68 题题4 线性方程组求解数值方法 3.3 向量范数与矩阵范数向量范数与矩阵范数 误差的度量误差的度量 3.3.1 向量范数向量范数

8、( (vector norms) ) 常用向量范数：常用向量范数： , | 1 1 n i xx,| 2/1 1 2 2 n i xx . |max 1 i ni xx 这些范数满足：这些范数满足： ; 0 ifonly and if 1 and , 0 ) 1 (xxx ;constantany is , )2(xx . )3(yxyx |limxx p p 一般范数：一般范数： 线性方程组求解数值方法 如，设向量如，设向量 x=(2, -4, -1)T, 则则 ,21| 2/1 1 2 2 n i xx, 7| 1 1 n i xx . 4|max 1 i ni xx 3.3 3 向量范数

9、与矩阵范数向量范数与矩阵范数 【Matlab】 x=(1:4)/5 N1=norm(x,1) N2=norm(x) N7=norm(x,7) Ninf=norm(x,inf) x = 0.2000 0.4000 0.6000 0.8000 N1 = 2 N2 = 1.0954 N7 = 0.8153 Ninf = 0.8000 线性方程组求解数值方法 定义定义 矩阵矩阵A A的范数定义为的范数定义为 3.3 3 向量范数与矩阵范数向量范数与矩阵范数 3.3.2 矩阵范数矩阵范数 ( (matrix norms) ) .max 0 x Ax A x 命题：矩阵的常用范数：命题：矩阵的常用范数：

10、; |max ) 1 ( 1 1 1 n i ij nj aA ; |max )2( 1 1 n j ij ni aA , )(max )3( 1 2 AAA T i ni 其中，其中，1, , n是是 AAT 的特征值，的特征值， i ni T AA 1 max: )(称为称为 的谱半径。的谱半径。 AAT 线性方程组求解数值方法 例例 设矩阵设矩阵 3.3 3 向量范数与矩阵范数向量范数与矩阵范数 , 43 21 A ; 6|max ) 1 ( 1 1 1 n i ij nj aA ; 7|max )2( 1 1 n j ij ni aA .46. 522115)(max )3( 1 2

11、AAA T i ni 则则 可以证明可以证明，矩阵范数满足：，矩阵范数满足： ; 0 ifonly and if 0 ; 0 (1)AAA ;|c| (2)AcA ; (3)BABA (4) .ABAB 线性方程组求解数值方法 3.4 经典迭代法经典迭代法 Jacobi迭代法与迭代法与Gauss-Seidel迭代法迭代法 3.4.1Jacobi迭代法迭代法 设有设有n阶方程组阶方程组 nnnnnn nn nn bxaxaxa bxaxaxa bxaxaxa 2211 22222121 11212111 (3.4.1) 线性方程组求解数值方法 若系数矩阵非奇异，且若系数矩阵非奇异，且 (i =

12、1, 2, n)，将方程组，将方程组 0 ii a 11,22111 23231212 22 2 13132121 11 1 1 1 1 nnnnnn nn n nn nn xaxaxab a x xaxaxab a x xaxaxab a x （4.1）改写成）改写成 3.4 经典迭代法经典迭代法 线性方程组求解数值方法 然后写成迭代格式然后写成迭代格式 )( 11, )( 22 )( 11 )1( )( 2 )( 323 )( 1211 22 )1( 2 )( 1 )( 313 )( 2121 11 )1( 1 1 1 1 k nnn k n k nn nn k n k nn kkk k

13、nn kkk xaxaxab a x xaxaxab a x xaxaxab a x （3.4.2） （3.4.2）式也可以简单地写为）式也可以简单地写为 ), 2, 1( 1 )( 1 )1( nixab a x k j n ij j iji ii k i （3.4.3） 3.4 经典迭代法经典迭代法 线性方程组求解数值方法 写成矩阵形式：写成矩阵形式： A = L U D Bf Jacobi 迭代阵迭代阵 bDxULDx kk 1)(1)1( )( (3.4.6) bxULxD bxULDbxA )( )( bDxULDx 11 )( 3.4 经典迭代法经典迭代法 线性方程组求解数值方法

14、)( 1 1 )( 1 )( 414 )( 313 )( 212 11 )1( 1 bxaxaxaxa a x k nn kkkk )( 1 2 )( 2 )( 424 )( 323 )1( 121 22 )1( 2 bxaxaxaxa a x k nn kkkk )( 1 3 )( 3 )( 434 )1( 232 )1( 131 33 )1( 3 bxaxaxaxa a x k nn kkkk )( 1)1( 11 )1( 33 )1( 22 )1( 11 )1( n k nnn k n k n k n nn k n bxaxaxaxa a x 只存一组向量即可。只存一组向量即可。 写成矩

15、阵形式：写成矩阵形式： bDxUxLDx kkk 1)()1(1)1( )( bxUxLD kk )()1( )( bLDxULDx kk 1)(1)1( )()( Bf Gauss-Seidel 迭代阵迭代阵 3.4.2 Gauss-Seidel迭代法迭代法 (3.4.7) 3.4 经典迭代法经典迭代法 线性方程组求解数值方法 定理定理3.1 对任意初始向量对任意初始向量X(0)及常向量及常向量F，迭代格式迭代格式(3.5.1) FBXX kk )1( (3.5.1) 推论推论：若迭代矩阵若迭代矩阵的的某种范数某种范数 ，则则（3.5.1） 1B 收敛的充分必要条件是迭代矩阵收敛的充分必要条

16、件是迭代矩阵B B的谱半径的谱半径 (B) 1) （此时矩阵也称为“病态矩阵”）的线性方程组Ax=b是病态的； 反之，系数矩阵的条件数很小的线性方程组Ax=b是良态的。 A=hilb(n); c=cond(A) Example: Example: 著名的Hilbert 矩阵是病态的。 【Matlab】 3.7 7条件数、病态方程组与算法的稳定性条件数、病态方程组与算法的稳定性 线性方程组求解数值方法 定理定理3.6 ( (条件数的性质条件数的性质) ) 设A 是一个非奇异矩阵. (1) cond(A)1, cond(A)=cond(A-1), con(aA)=cond(A), a 0. (2) 如果 Q是正