数值分析(05)高斯消元法.pdf

数值分析数值分析 数值分析数值分析 第二节第二节 高斯消元法高斯消元法 第三节第三节 矩阵的三角分解法矩阵的三角分解法 第四节第四节 误差分析和解的精度改进误差分析和解的精度改进 第五节第五节 大型稀疏方程组的迭代法大型稀疏方程组的迭代法 第三章第三章 线性代数方程组的数值解法线性代数方程组的数值解法 n 第一节第一节 引言引言 第六节第六节 极小化方法极小化方法 数值分析数值分析 数值分析数值分析 线性代数方程组的一般形式线性代数方程组的一般形式 (1)mn Axb AR 用用矩矩阵阵形形式式表表示示为为 其其增增广广矩矩阵阵记记为为 11112211 21122222 1122 nn nn mmmnnm a xa xa xb a xa xaxb axaxaxb 111211 212222 12 , n n mmmnm aaab aaab AA b aaab 第一节第一节 引引 言言 数值分析数值分析 数值分析数值分析 AxbAA 线线性性方方程程组组 有有解解的的充充分分必必 结结论论1 1 （线线性性代代数数方方程程组组 要要条条件件是是：秩秩（ ）= =秩秩 有有解解判判定定理理 （ 别别） ） (1)( )( ), Axb AArnAxb 线线性性方方程程组组有有解解（即即相相容容）时时， 秩秩 结结论论2 2 秩秩则则方方程程组组存存在在唯唯一一解解。 (2) ( )( ),r Ar ArnAxb 方方程程组组有有无无穷穷多多解解。 通通解解原原方方程程组组一一个个特特解解对对应应齐齐次次方方程程组组的的基基 础础解解系系的的线线性性组组合合。 22 2 ,|min| mn Axb xxxxAxb 常常见见是是，称称为为欠欠定定方方程程组组（方方程程数数少少于于未未知知数数） 此此时时，从从的的无无穷穷多多个个解解中中需需求求出出范范数数最最小小的的解解。 即即求求使使， 满满足足。 数值分析数值分析 数值分析数值分析 2 2 ()() () |min r Ar AAxb mn bAR A x xbAx 方方程程组组无无解解（即即不不相相容容）。 常常见见是是，称称为为超超定定方方程程组组（又又称称矛矛盾盾方方程程组组） 此此时时，向向量量 不不在在 的的列列空空间间之之中中，原原方方程程组组 无无解解，但但可可求求出出最最小小二二乘乘意意义义下下的的解解 。 即即求求 使使 MATLAB实现实现: x=Ab 11112211 21122222 1122 nn nn nnnnnn n a xa xa xb a xa xaxb a xaxaxb 本本章章介介绍绍求求解解 阶阶线线性性方方程程组组的的数数值值方方法法 数值分析数值分析 数值分析数值分析 数值求解方法有以下三条途径数值求解方法有以下三条途径 直接法：利用直接法：利用Gauss消元或消元或矩阵分解，通过有限次运算矩阵分解，通过有限次运算 可求出精确解。可求出精确解。 迭代法：构造迭代格式，产生迭代序列，通过无限迭代法：构造迭代格式，产生迭代序列，通过无限 次迭代过程求解。有限次截断得近似解。次迭代过程求解。有限次截断得近似解。 极小化方法：构造二次模函数，用迭代过程求二次极小化方法：构造二次模函数，用迭代过程求二次 模函数的极小化问题，即变分法（经模函数的极小化问题，即变分法（经n 次运算，理论上得精确解）要求次运算，理论上得精确解）要求A 对称正定对称正定(S.P.D) 数值分析数值分析 数值分析数值分析 用增广矩阵表示为用增广矩阵表示为 同解同解 初等变换初等变换 组组化为同解的上三角方程化为同解的上三角方程将原方程组将原方程组 求解求解 gUxbAx gUxbAx RAbAx nn 第二节第二节高斯消元法高斯消元法 数值分析数值分析 数值分析数值分析 三角形方程组包括上三角形方程组和下三角形三角形方程组包括上三角形方程组和下三角形 方程组方程组，是最简单的线性方程组之一是最简单的线性方程组之一。上三角方上三角方 程组的一般形式是程组的一般形式是： ),.,2 , 1(0 . . . 11111 22222 11212111 nia bxa bxaxa bxaxa bxaxaxa ii nnnn nnnnnnn nn nn 其中其中 一、三角形方程组的解法一、三角形方程组的解法 数值分析数值分析 数值分析数值分析 124 234 34 4 4 57 31313 1313 xxx xxx xx x 用用回回代代法法求求解解线线性性方方例例程程组组 4 34 243 142 1234 1 (1313)/30 ( 75)( 750)2 44121 ,)(1 ,2,0 : ,1) TT x xx xxx xxx xxxx 所所以以，解解为为（ 解解 数值分析数值分析 数值分析数值分析 1 / ()/1,1 nnnn n iiikkii k i xba xba xain 为求解上三角方程组，从最后一个方程入手，先为求解上三角方程组，从最后一个方程入手，先 解出解出 xn=bn/ann, 然后按方程由后向前的顺序，从方程然后按方程由后向前的顺序，从方程 中依次解出中依次解出xn-1,xn-2,x1。 。这样就完成了上三角方程组 这样就完成了上三角方程组 的求解过程。这个过程被称为回代过程其计算步骤如的求解过程。这个过程被称为回代过程其计算步骤如 下：下： 11112211 22222 11111 . . . 0(1,2,., ) nn nn nnnnnnn nnnn ii a xa xa xb a xaxb axaxb a xb ain 其其中中 数值分析数值分析 数值分析数值分析 function X=backsub(A,b) %InputA is an nn upper- triangular nonsingullar matrix % -b is an n1 matrix %OutputX is the solution to the system AX=b 函数名返回变量参数表 n=length(b); X=zeros(n,1); X(n)=b(n)/A(n,n); for i=n-1:-1:1 X(i)=(b(i)-A(i,i+1:n)* X(i+1:n)/A(i,i); end A的第i行、第i+1到n列元素 构成的行向量 1 / ()/1,1 nnnn n iiikkii k i xba xba xain 数值分析数值分析 数值分析数值分析 124 234 34 4 4 57 31313 1313 xxx xxx xx x 求求解解上上三三角角方方程程组组例例 4 : 1x 解解 2 2x 12 23 3 3 2 30 xx xx x 3 0 x 12 2 3 2 xx x 1 1234 1 ,)(1 ,2,0 ,1) TT x xxxx 所所以以，解解为为 （ 数值分析数值分析 数值分析数值分析 for i= n : 1 : 2 b ( i ) = b ( i ) / A( i , i ); b (1: i - 1 ) = b (1: i - 1 ) - b ( i ) *A(1: i - 1 , i ) ; end b ( 1 ) = b ( 1 ) / A( 1 ,1 ); 求解上三角方程组求解上三角方程组 Ax=b 1 2 1 1 1)(1) 2 1 1)(1) 2 n i n i nin n in nn 求求解解一一个个三三角角形形方方程程组组需需 次次除除法法, , （次次乘乘法法, , （次次加加法法。总总计计算算量量约约为为 。 数值分析数值分析 数值分析数值分析 1 12 123 22 2 3249 x xx xxx 用用回回代代法法求求解解线线性性方方程程组组例例 1 2 3 123 2/21 (21)/11 (93 12 1)/41 ,)(1 ,1 , ) : 1 x x x xxx 所所以以，解解为为（ 解解 数值分析数值分析 数值分析数值分析 12 1111 1 1 , / ()/(2,3, ) n i iiikkii k xxx xba xba xain 下下三三角角形形方方程程组组可可以以参参照照上上三三角角形形方方程程组组的的解解法法来来求求解解， 下下三三角角形形方方程程组组的的求求解解顺顺序序是是从从第第一一个个方方程程开开始始，按按从从上上到到下下 的的顺顺序序，依依次次解解出出：其其计计算算公公式式为为： 如如上上解解三三角角形形方方程程组组的的方方法法称称为为回回代代法法. . 1111 2112222 1122 0,1,2, nnnnnn ii a xb a xa xb a xaxaxb ain 下三角方程组的一般形式为：下三角方程组的一般形式为： 其中其中 数值分析数值分析 数值分析数值分析 高斯消元法是一个古老的直接法高斯消元法是一个古老的直接法, ,由它改进得到由它改进得到 的选主元法的选主元法, ,是目前计算机上常用于求低阶稠密矩阵是目前计算机上常用于求低阶稠密矩阵 方程组的有效方法方程组的有效方法, ,其特点就是通过消元将其特点就是通过消元将一般线性一般线性 方程组方程组的求解问题转化为的求解问题转化为三角方程组三角方程组的求解问题。的求解问题。 高斯消元法的求解过程高斯消元法的求解过程, ,可大致分为两个阶段可大致分为两个阶段: : 首先首先, ,把原方程组化为上三角形方程组把原方程组化为上三角形方程组, ,称之为称之为“消消 元元”过程”过程; ;然后然后, ,用逆次序逐一求出上三角方程组用逆次序逐一求出上三角方程组( (原原 方程组的等价方程组方程组的等价方程组) )的解的解, ,称之为称之为“回代”过程“回代”过程. . 高斯高斯“消消元元”过程”过程可通过矩阵运算来实现。具可通过矩阵运算来实现。具 体过程如下：体过程如下： 二、顺序高斯消元法二、顺序高斯消元法 数值分析数值分析 数值分析数值分析 123 123 123 236 2349 31 326 Gauss xxx xxx xxx 用用顺顺序序消消元元法法求求解解方方例例程程组组 11/1/ 21/2/ 01, 3 6231 9432 6321 113131 112121 11 )1( aam aam an bAA增广矩阵：增广矩阵： 解：解： 11 1 21,: 11 L AxL b 1 1 L =,L =,完完成成第第一一步步消消元元 得得 数值分析数值分析 数值分析数值分析 (2)(2)(2) 22323222 22121 10,/1/( 1)1 1 1, 11 amaa LL L AxL L b =,=,完完成成第第二二步步消消元元 得得 33 32 632 3 32 321 x xx xxx 3 23 123 123 3/31 ( 32)( 32 1)1 62362 13 11 1,1,1 x xx xxx xxx 回回代代求求得得 故故所所求求解解为为 0110 3210 6321 )2( A 3300 3210 6321 )3( A 数值分析数值分析 数值分析数值分析 将方程组将方程组Ax=b的系数矩阵与右端项合并为的系数矩阵与右端项合并为 111211 212222 12 , n n nnnnn aaab aaab A bA aaab (1)(1)(1) 1111 (1) (1)(1)(1)(1) 12 (1)(1)(1) 1 . ,., . n n nnnn aab AAb aab 记记 (1)(1) 1 (1) 11111 ,0,.,0. T A LLa 对的第一列构对的第一列构 造使造使 1 (1)1 1111 11 0,2,., i i a amin a （ ）（ ） （ ）（ ） ：设取：设取第一步第一步 21 1 1 1 1 1 n m L m 数值分析数值分析 数值分析数值分析 (1)(1)(1)(1) 111211 (2)(2)(2) (1)(2) (2)(2)(2)(2)2222 112 (2)(2)(2) 2 . 0. ,., 0. n n n nnnn aaab aab L AAb aab (2)(1)(1) 11 (2)(1)(1) 1 1 2, ,2, ,2, ijijij iii aam ain jn bbm bin (1)(1) 1 (1)(1) 11 AxbL L AxLb 对对方方程程组组从从左左边边乘乘以以 (1)(1)(1) 1111 (1) 21 1 (1)(1)(1) 1 1 1 . 1 . 1 n nnnn n aab m L A aab m 数值分析数值分析 数值分析数值分析 (2) (2)2 222(2) 22 03,., i i a amin a ：设，：设，第二步第二步取取 (2) (2) 3222 2 (1)(1)(1)(1)(1) 1112131,1 (2)(2)(2)(2) 22232,2 (2)(1) (3)(3)(3)(3) 221333,3 (3)(3)(3) 3, 1 1 1, 1 0 00 00 n n n n nn nn mAL m aaaab aaab L AL L AaaAb aab - -对的第二列构造对的第二列构造 - - 使使 数值分析数值分析 数值分析数值分析 (2) 2 2(2) 22 , i i a m a (1)(1) 2121 L L A xL Lb (3)(2)(2) 22 (3)(2)(2) 22 ,3, 3,4,., ijijij iii aam ai jn bbm bin (1)(1)(1)(1) 111211 (2)(2)(2) (2) 2222 322 (2)(2)(2) 2 2 (1)(1)(1)(1) 1112131, (2)(2)(2) 22232, (3)(3) 333, (3 3 1 . 1 0. 1 0. 1 0 00 00 n n nnnn n n n n n aaab aab mL A aab m aaaa aaa aa a - - - - (1) 1 (2) 2 (3)(3) 3 )(3)(3) ,n nn b b Ab ab 数值分析数值分析 数值分析数值分析 进行到第进行到第k步消元时步消元时 ( )(1)( )kkk AAAk 下下一一步步消消元元，从从，将将的的第第 列列的的对对角角元元 以以下下的的元元素素化化为为零零。 (1)(1)(1)(1)(1) 11121311 (2)(2)(2)(2) 222322 (3)(3)(3) 3333 ( ) ( )( )( ) 1 ( )( )( ) 1,1,11 ( )( )( )( ) ,1 . . . . . . n n n k kkk kkkkk kkk kkkkk kkkk nkn knnn aaaab aaab aab A aab aab aaab 数值分析数值分析 数值分析数值分析 ( ) ( ) ( ) 0, (1,., ) k k ik k kkik kk k a am a ikn GaussL 设取设取 ， 构造变换阵构造变换阵 , , 1 1 1 1 1 1 T k k kk nk Il e m m (1)( )kk k AL A 消元计算递推公式：消元计算递推公式： ( ) (1,2,1) k kk a kn 称为主元素.称为主元素. ( )( ) (1)( )( ) (1)( )( ) 1, 1/ 21, 3 kk ikikkk kkk ijijikkj kkk iiikk ikn maa aam ajkn bbm b （ ）（ ） （ ），（ ）， （ ）（ ） 数值分析数值分析 数值分析数值分析 (1)(1)(1)(1) 11112211 (2)(2)(2) 22222 ( )( ) nn nn nn nnnn axaxaxb axaxb axb 即即 用回代过程求解上三角方程组，即可得解向量用回代过程求解上三角方程组，即可得解向量 ( x1*,x2*, ,xn*) )T T. . 是高斯消元的前提。是高斯消元的前提。)1,2, 1( ,0 )( nka k kk (1)(1) 121121nn LL L A xLL Lb (1)(1)(1)(1) 111211 (2)(2)(2) ( )2222 ( )( ) 0 00 n nn nn nnn aaab aab A ab 最后得最后得 数值分析数值分析 数值分析数值分析 求解的全过程包括两个步骤：消元和回代求解的全过程包括两个步骤：消元和回代 1 .1 .顺序消元顺序消元 2 . 回代求解回代求解 ( )( ) ( )( )( ) 1 / ()/,1,2,1 nn nnnn n kkk kkkjjkk j k xba xbaxaknn ( )( ) (1)( )( ) (1)( )( ) 1,1 1, 1/ 21, 3 kk ikikkk kkk ijijikkj kkk iiikk kn ikn maa aam ajkn bbm b （ ）（ ） （ ），（ ）， （ ）（ ） 数值分析数值分析 数值分析数值分析 步消元计算后步消元计算后，第，第的二维数组存放的二维数组存放一个一个 用用用动态存储方式。最初用动态存储方式。最初在计算机中计算时，采在计算机中计算时，采 存储方式存储方式 kAnn ), 1;, 1( ), 1( )1()( )( nkjnkiaa nkima k ij k ij ik k ik )()( 1, )( )( 1,1 )( ,1 )( 1 )( )3( 3 )3( 33 )2( 2 )2( 23 )2( 22 )1( 1 )1( 13 )1( 12 )1( 11 )( . . . . . . k nn k kn k nk k kk k kk k kk k kk n n n k aaa aa aa aa aaa aaaa A ik m ) 1( k ij a 数值分析数值分析 数值分析数值分析 UIL 消元过程全部完成后，原来的二维数组中存放的消元过程全部完成后，原来的二维数组中存放的 元素实际上是一个新的矩阵，记为元素实际上是一个新的矩阵，记为 F A F AA用动态形式表示为用动态形式表示为 )( 1,321 )3( 3 )3( 1, 3 )3( 333231 )2( 2 )2( 1, 2 )2( 23 )2( 2221 )1( 1 )1( 1, 1 )1( 13 )1( 12 )1( 11 n nnnnnnn nn nn nn F ammmm aaamm aaaam aaaaa A 数值分析数值分析 数值分析数值分析 function X=gauss(A,b) %InputA is an nn nonsingullar matrix % -b is an n1 matrix %OutputX is the solution to the system AX=b MATLAB For Gaussian Elimination n n=size(A); % 确定确定A的维数的维数 X=zeros(n,1); for k=1:n-1 for i=k+1:n % % 消元过程消元过程 A(i,k) =A(i,k)/ A(k,k); % A(k,k) 0 A(i,k+1:n)= A(i,k+1:n)-A(i,k) *A(k,k+1:n); b(i)= b(i)-A(i,k) *b(k); end end X=backsub(A, b); %回代求解回代求解 数值分析数值分析 数值分析数值分析 高斯消元法的计算量分析高斯消元法的计算量分析 33 2 1 333 nnn Nn 加减法计算量加减法计算量为为 乘除乘除法计算量法计算量为为 3 1 (1)(25)/6 3 n Nn nn 高斯消元法的计算量高斯消元法的计算量约约为为 3 2 3 n 数值分析数值分析 数值分析数值分析 定理定理 3-1 ),.2 , 1( )( nka k kk 全不为零的充分必要条件是 A 的顺序主子式 k D全不为零，即),.2 , 1( , 0nkDk 。 )()2( 22 )1( 11 )( )2( 2 )2( 22 )1( 1 )1( 12 )1( 11 1 111 . . . . . . k kk k kk k k kkk k k aaa a aa aaa aa aa D 证明证明：由于对系数矩阵 A 的消元计算不改变 A 的行列式的值， 因此 A 的顺序主子式为 当),.,2 , 1(0 )( nka k kk 时，显然),.2 , 1( , 0nkDk ， 必要性得证。 用归纳法证充分性。 数值分析数值分析 数值分析数值分析 定理定理3 3- -2 2 对线性代数方程组对线性代数方程组Ax=bAx=b，其中，其中A A非奇异，非奇异， 若系数矩阵若系数矩阵A A的顺序主子式全不为零，则可用顺序高的顺序主子式全不为零，则可用顺序高 斯消元法求解线性代数方程组斯消元法求解线性代数方程组Ax=bAx=b。 消元法是解线性方程组的基本方法，具有计算简消元法是解线性方程组的基本方法，具有计算简 单的优点，但有时由于主元过小，使得计算结果严重单的优点，但有时由于主元过小，使得计算结果严重 失真失真, ,实际中常采用选主元高斯消元法。实际中常采用选主元高斯消元法。 数值分析数值分析 数值分析数值分析 例例3-2 讨论下面方程组的解法讨论下面方程组的解法 0.0001x1+x2=1 x1+x2=2 假设求解是在四位浮点十进制数假设求解是在四位浮点十进制数 的计算机上进行的计算机上进行 0.1000 10-3x1+ 0.1000 101 x2= 0.1000 101 0.1000 101 x1 +0.1000 101 x2= 0.2000 101 解解: :本题用机器数系表示为本题用机器数系表示为 a11 =0.0001, m21=a21/a11=1/0.0001= 104, 消元消元得得 回代解得回代解得 x2 2=1=1, x1 1=0=0严重失真严重失真! ! (本题的准确解为本题的准确解为 x1 1= = 10000/9999, x2 2= =9998/9999 ） a22(2)= 0.1000 101- 104 0.1000 101 = 0.00001 105- 0.1000 105(对阶计算对阶计算) = - 0.1000 105 0.1000 10-3 x1+ 0.1000 101 x2= 0.1000 101 -0.1000 105 x2= -0.1000 105 主元主元a11过小过小 数值分析数值分析 数值分析数值分析 选主元选主元高斯消元法的高斯消元法的基本思想基本思想 用高斯消元法求解线性方程组时用高斯消元法求解线性方程组时, ,为避免小的主元为避免小的主元. . 在进行第在进行第k k步消元前步消元前, ,应该在第应该在第k k列元素列元素( (i=k,ni=k,n) ) 中找出第一个出现的绝对值最大者中找出第一个出现的绝对值最大者, ,例如例如再再 把第把第i ik k个方程与第个方程与第k k个方程组进行交换个方程组进行交换, ,使使成为主元成为主元. . 我们称这个过程为选主元我们称这个过程为选主元. .由于只在第由于只在第k k列元素中选主元列元素中选主元, , 通常也称为通常也称为按列选主元按列选主元. . )(k ik a ( )( ) max k kk i kik k i n aa ( ) k k i k a 三、选主元三、选主元高斯消元法高斯消元法 )(k ij a 如果在第如果在第k k步消元前步消元前, ,在第在第k k个方程到第个方程到第n n个方程所有个方程所有 的的x xk k到到x xn n的系数的系数( (i=k,n;j=k,ni=k,n;j=k,n) )中中, ,找出绝对值找出绝对值 最大者最大者, ,例如例如，再交换第，再交换第k k, ,i ik k两个两个 ( )( ) , max k k kk i jij k i j n aa 数值分析数值分析 数值分析数值分析 方程和第方程和第k k, ,j jk k列列, ,使使成为主元成为主元. . 称这个过程为称这个过程为完全完全 选主元选主元. 不论是哪种方式选出主元不论是哪种方式选出主元, ,而后再按上面介绍的计算步而后再按上面介绍的计算步 骤进行消元的计算骤进行消元的计算, ,一般都称为选主元的高斯消元法一般都称为选主元的高斯消元法. .在在 实际计算中实际计算中, ,常用按列选主元的高斯消元法常用按列选主元的高斯消元法. . ( ) k k k i j a 数值分析数值分析 数值分析数值分析 ( )( ) ( ) ( )( ) () ()| max| , ,1, , (1) ) , ( , , 2 k k kk kkk kk kki kik k i n k kk k k kji j kjkji ji j kikkii k ikinaa AAbik Tik aajk kn TaaaaT bbTb bbbT 对对每每一一步步 第第 步步 消消元元，分分两两步步 确确定定使使 对对增增广广矩矩阵阵使使 列列主主元元高高斯斯消消元元法法具具体体做做法法是是： 选选列列主主元元换换行行 在在计计算算机机上上，用用一一个个工工作作单单元元来来完完成成，对对， 包包括括 消消元元计计算算 数值分析数值分析 数值分析数值分析 列主元高斯消元法解线性方程组列主元高斯消元法解线性方程组Ax=b的算法步骤的算法步骤 停机。停机。信息信息 输出失败输出失败则认为则认为如果如果 使使确定确定 、选列主元、选列主元 步。步。循环执行到第循环执行到第对对 、置、置 , , 0det , 0 , max 2 51, 2 , 1 1det 1 kiki ik nik kik kk k aa aai nk detdet , ), 1,( , 4 , 3 k k ik jikj k bb nkkjaa ki 否则交换行否则交换行步步转出执行第转出执行第、如果、如果 具体执行行交换要通过工作单元具体执行行交换要通过工作单元 T。 TbbbbT TaaaaT kk kk iikk jijikjkj ; ; ; ; 数值分析数值分析 数值分析数值分析 。、输出解向量、输出解向量 、 否则否则 停机。停机。输出失败信息输出失败信息则认为则认为如果如果 、回代求解、回代求解 、 FT n nn ii n ij jijii nnnn nnnn kk AAbbbx a nniababb abb aa a det, ,),( 8 detdet7 )1 , 2 , 2, 1( / )( / , 0, 6 detdet 5 21 1 (3) ), 1( (2) /(1) , 2, 1 4 kikii kjikijij kkikikik babb nkjaaaa aama nkki 、消元计算、消元计算 数值分析数值分析 数值分析数值分析 算法算法3 3- -2 2列主元高斯消元法解线性方程组列主元高斯消元法解线性方程组Ax = bAx = b function X,det=gauss(A,b)function X,det=gauss(A,b) %输入输入A A 是一个是一个 n nn n 非奇异矩阵，非奇异矩阵，b b 是一个列向量。是一个列向量。 %输出输出X X是是 AX=bAX=b的解，的解，detdet是是A A的行列式。的行列式。 n n=size(A); % n n=size(A); % 确定确定A A的维数的维数 det=1; det=1; for k=1:nfor k=1:n-