第二章 解线性方程组的直接法 2.3 Gauss列主元 消去法§ 2.3 Gauss列主元 消去法§ 例 1. 用 Gauss消去法解线性方程组 (用 3位十进制浮 点数计算 ） ?? ? =+ =+ 2 10001.0 21 21 xx xx 解 : 本方程组的精度较高的解为 Tx )99989999.0,00010001.1(* = 用 Gauss消去法求解 (用 3位十进制浮点数计算 ) 一 、 Gauss列主元消去法的引入 ),( bAA = ?? ? ? ??? ?= 2 1 11 1000100.0 ??? →? =1000021m ?? ? ? ??? ? ×?×? 44 1000.1 1 1000.10 1000100.0 9999 00.1,00.0 21 == xx回代后得到 与精确解相比 ,该结果相当糟糕 究其原因 ,在求行乘数时用了很小的数 0.0001作除数 主元 ),( bAA = ?? ? ? ??? ?→ 1 2 1000100.0 11 ??? →? = 0001.021m ?? ? ? ??? ? 00.1 2 00.10 11 如果在求解时将 1,2行交换 ,即 0.9999 00.1,00.1 21 == xx 回代后得到 这是一个相当不错的结果 ),( bAA = 例 2. 解线性方程组 (用 8位十进制尾数的浮点数计算 ) ?? ? ? ? ?? ? ? ? = ?? ? ? ? ?? ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? 3 2 1 643.5072.12 623.4712.31 3210 3 2 1 8 x x x 解 : 这个方程组和例 1一样 ,若用 Gauss消去法计算会有 小数作除数的现象 ,若采用换行的技巧 ,则可避免 ? ? ? ? ? ? ? ? ? ? ? ?= ? 3 2 1 643.5072.12 623.4712.31 3210 8 行交换因此 的列元素为 绝对值最大很小 3,1 ,2 ,10 13 8 ?= ? a ?? →? ? 31 rr ? ? ? ? ? ? ? ? ? ? ? ? ? 1 2 3 3210 623.4712.31 643.5072.12 8 ???? →? ?×?= = 8 31 21 105.0 5.0 m m ? ? ? ? ? ? ? ? ? ? ××× ×× ? 101.0 5.0 3 103.0102.00 1018015.0103176.00 643.5072.12 绝对值最大 不需换行 ???? →? = 92722629.032m ? ? ? ? ? ? ? ? ? ? × ×× ? 54138685.0 5.0 3 1041555186.000 1018015.0103176.00 643.5072.12 ),( )1()1( bA= ),( )2()2( bA= ),( )3()3( bA= )3( )3( 3 3 33a bx = = 经过回代后可得 )1( 11 3 )1( 132 )1( 12 )1( 1 1 a xaxabx ??= 54138685.0 1041555186.0 × 39257367.0= )2( 22 3 )2( 23 )2( 2 2 a xabx ?= 103176.0 1018015.05.0 3 × ××?= x 05088607.0?= 49105820.0?= 事实上 ,方程组的准确解为 Tx )367257384.0,050886075.0,491058227.0(* ??= ),( )1()1( bA 例 2所用的方法是在 Gauss消去法的基础上 ,利用换行 避免小主元作除数 ,该方法称为 Gauss列主元消去法 二 、 Gauss消元过程与系数矩阵的分解 1.Gauss消去法消元过程的矩阵描述 ? ? ? ? ? ? ? ? ? ? ? ? = )1()1()1( 2 )1( 1 )1( 2 )1( 2 )1( 22 )1( 21 )1( 1 )1( 1 )1( 12 )1( 11 nnnnn n n baaa baaa baaa L MMMM L L niaam ii ,,3,2)1( 11 )1( 1 1 L== 行变换相 当于左乘 初等矩阵 由于 ? ? ? ? ? ? ? ? ? ? ? ? ? ?= 1 1 1 1 21 1 nm mL OM令 则 ),( )1()1(1 bAL ? ),( )2()2( bA= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= + 1 1 1 1 , ,1 kn kk k m mL OM O 显然若令 则有 ),( )()( kkk bAL ? ),( )1()1( ++= kk bA 1,,3,2,1 ?= nk L ),( )1()1(1 bAL ??? 2LL??1nL ),( )()( nn bA=因此 从而 AL ?1?? 2LL??1nL )(1 1 1 2 1 1 n n ALLL ? ? ??= LA故 为上三角矩阵)(nAU = 为单位下三角矩阵1 11211 ????= nLLLL L LU= )(nA= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? ??? 1 1 1 1 1 1,3,2,1, 2,12,11,1 3231 21 nnnnn nnn mmmm mmm mm m L L L OMMM即 ?? ? ? ? ? ? ?? ? ? ? ? ? = )( )2( 2 )2( 22 )1( 1 )1( 12 )1( 11 n nn n n a aa aaa MO L L )(nAU = 且 Udet=Adet 顺序主元∏ = = n i i iia 1 )( 定义 1. 不带行交换的 Gauss 消去法的消元过程 ,产生 一个单位下三角矩阵 L和一个上三角矩阵 U,即 该过程称之为 LUA = .分解的矩阵 LUA 由上述分析不难得到 ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? = nnnkn knkkk nk aaa aaa aaa A LL MOMM LL MMOM LL 1 1 1111 ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? = 1 1 1 1 1 LL OMM L OM nkn k mm m ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? )( )()( )1( 1 )1( 1 )1( 11 n nn k kn k kk nk a aa aaa MO L MMO LL 阶顺序主子式k kA kkk ULA = kAdet kUdet= ∏ = = k i i iia 1 )( kL kU nk Ak ,,2,1 0det L= ≠ ni a iii ,,2,1 0)( L= ≠Gauss消去法 可以执行 定理 1. ,0det ≠= kk ADAn 的顺序主子式阶方阵若 ,1,,2,1 ?= nk L 存在且唯一分解结果的则 LUALUA = 在定理中 ,可能注意到 0det == nn AD 0)( =nnna即可能存在 分解并不影响的这对 LUA 起决定作用消去法的回代能否进行用但对方程组 GaussbAx = 2.Gauss列主元消去法消元过程的矩阵描述 由于 Gauss列主元消去法每一步都要选取列主元 ,因 此不可避免要进行行交换 行交换后再消元行与的第要将步消元时设第 kkk ikbAk ),(, )()( 即 ),( )()( kk bA ),( )1()1( ++= kk bA?kL ?kikI , kikI , ?? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? = 1 01 10 1 O L MOM L O 行第 k 行第 ki kik ≤一般 时kik = II ki k =, 表示不换行初等矩阵 因此 ,Gauss列主元消去法的消元过程为 : ),( )1()1( bA ),( )()( nn bA=?? 1,1 1iIL?? 2,2 2iILL?? ?? ? 1,1 1 nin nIL A )( nA=?? 1,1 1i IL?? 2,2 2i ILL?? ?? ? 1,1 1 nin n IL显然 U= 上三 角阵消元过程为时当例如 ,4, =n A )4(A=?? 1,1 1i IL?? 2,2 2i IL?? 3,3 3i IL U= 3L 3,23, 33 ii ILI? 3,2,12,3, 3223 iiii IILII? 1,2,3, 123 iii III? A? U= 3,23,2 33 ~ ii ILIL =设 3,2,12,3,1 3223 ~ iiii IILIIL = 仍然为单位 下三角矩阵 1,2,3, 123 iii IIIP = 初等矩阵的乘积 ,称为 排列阵 则 3L 2~L? 1~L? P? A? U= 推广到一般情形 1?nL 2 ~L? 1 ~L? P? A? U=L 2 ~ ?? nL =PA 1( ?nL 2~L? 11 )~ ?? LL2~ ?? nL U? 1( ?nL 2 ~L? 1 1 ) ~ ?? LL 2 ~ ?? nL令 =L 仍然为单位 下三角矩阵 =PA LU则 单位下三角阵与上三角阵的乘积 分解的上述过程称为矩阵 LUPA 综合以上讨论 ,有 定理 2. ,0det, ≠AAn 即为非奇异矩阵阶方阵若 使得阵和一个非奇异上三角矩 、 一个单位下三角矩阵则必存在一个排列矩阵 ,U LP =PA LU 请作出 Gauss列主元消去法的程序 (用 Matlab语言 ) 已编程序 gaussliezhuyuan.m 开始 EPSnbA ,,,输入 输出无解信息 k?1 x输出解 EPSnnA ≤|),(| nk = kk ?+ 1 P选取主元素 EPSP ≤|| 消元换行 停机 回代求解 T T T F F F 三 、 Gauss 列主元消去法的算法设计 (一 ) 流程图 (二 ) 自然语言 n输入方程组的维数.1 1,,2,1,,,2,1,),( +== njniabA ij LL的元素增广矩阵 EPS控制条件转移精度 1,2,1.2 ?= nk L对于 PkkA ?),(1.2 nkki ,,1,2.2 L+=对于 |||),(| PkiA >如果 ).(7,||3.2 输出无解信息则转到如果 EPSP ≤ 5.2,,),( 0 否则转则 IiPkiA ?? 选主元 1,,1,4.2 ++= nkkj L对于 wjkA ?),( ),(),( 0 jkAjIA ? ),( 0 jIAw ? ,,,16.2 nki L+=对于 ,1,,1 ++= nkj L对于 ),(),(/),( kimkkAkiA ? ),(),(*),(),( jiAjkAkimjiA ?? ).(7,|),(|.3 无解则转到如果 EPSnnA ≤ ).(7,||5.2 输出无解信息则转到如果 EPSP ≤ 换行 消元 )(),(/)1,(.4 nxnnAnnA ?+ 1,2,,1.5 L?= nk对于 )(n nn n n a bx = )( 1 )()( i ii n ij j i ij i i i a xab x ∑ += ? = ),( )(),()1.( 1 kkA jxjkAnkA x n kj k ∑ += ?+ = 0=S nkj ,,1 L+=对于 SjxjkAS ?+ )(*),( SSnkA ??+ )1,( )(),(/ kxkkAS ? .8)),(,),2(),1((.6 并转输出解 nxxxxT L= 回代 输出无解信息.7 停机.8 上述过程中的储存空间需要 : 注意 : 1:1,:1),,( +== njnijiA增广矩阵 1:1,:1),,( ?=+= nknkikim行乘数矩阵 nkkx :1),( =解向量 SIPkji ,,,,, 0 储存空间需要较大 (三 ) 自然语言的改进 选主元 n输入方程组的维数.1 EPS控制条件转移精度 1,,2,1,,,2,1,),( +== njniabA ij LL的元素增广矩阵 1,2,1.2 ?= nk L对于 PkkA ?),(1.2 nkki ,,1,2.2 L+=对于 |||),(| PkiA >如果 5.2,,),( 0 否则转则 IiPkiA ?? .7,||3.2 则转到如果 EPSP ≤ 1,,1,4.2 ++= nkkj L对于 wjkA ?),( ),(),( 0 jkAjIA ? ),( 0 jIAw ? ,,,16.2 nki L+=对于 ,1,,1 ++= nkj L对于 ),(),(/),( kimkkAkiA ? ),(),(*),(),( jiAjkAkimjiA ?? .7,|),(|.3 则转到如果 EPSnnA ≤ .7,||5.2 则转到如果 EPSP ≤ ),(),(/),( kiAkkAkiA ? ),(),(*),(),( jiAjkAkiAjiA ?? 换行 消元 )(),(/)1,(.4 nxnnAnnA ?+ )1,(),(/)1,( +?+ nnAnnAnnA 1,2,,1.5 L?= nk对于 0=S 1:1,:1),,( +== njnijiA增广矩阵 1:1,:1),,( ?=+= nknkikim行乘数矩阵 nkkx :1),( =解向量 1:1,:1),,( +== njnijiA ?? ??? 回代 nkj ,,1 L+=对于 SjxjkAS ?+ )(*),( SnjAjkAS ?++ )1,(*),( SSnkA ??+ )1,( )(),(/ kxkkAS ? )1,(),(/ +? nkAkkAS .8)),(,),2(),1((.6 并转输出解 nxxxxT L= 输出无解信息.7 停机.8