第三章-解线性方程组的直接法ppt课件

第三章第三章 解线性方程组的直接法解线性方程组的直接法第三章第三章 解线性方程组的直接法解线性方程组的直接法3.1 引言引言n小行星轨道问题：小行星轨道问题：天文学家要确定一小行星的轨道，在轨道平面建天文学家要确定一小行星的轨道，在轨道平面建立以太阳为原点的直角坐标系。在坐标轴上取天文测立以太阳为原点的直角坐标系。在坐标轴上取天文测量单位量单位(一天文单位为地球到太阳的平均距离：一天文单位为地球到太阳的平均距离：9300万万英里，约英里，约1.5亿千米亿千米)，对小行星作，对小行星作5次观察，测得轨道次观察，测得轨道上上5个点的坐标数据如下：个点的坐标数据如下：x 5.7640 6.2860 6.7590 7.1680 7.4800 y 0.6480 1.2020 1.8230 2.5260 3.3600椭圆的一般方程椭圆的一般方程:a1x2+a2xy+a3y2+a4x+a5y+1=0将数据逐个代入，可得五个方程的方程组，求解该线性将数据逐个代入，可得五个方程的方程组，求解该线性方程组即可得行星轨道方程。方程组即可得行星轨道方程。对一般线性方程组对一般线性方程组:A x=b,其中其中nnnnnnaaaaaaaaaA212222111211nbbbb2112nxxxx由以前所学内容知，当且仅当矩阵由以前所学内容知，当且仅当矩阵A行列式不为行列式不为0时，即时，即A非奇异时，方程组存在唯一解，可根据非奇异时，方程组存在唯一解，可根据Cramer法则求解。法则求解。其算法设计如下：其算法设计如下：（1）输入系数矩阵输入系数矩阵A和右端向量和右端向量b；（2）计算系数矩阵）计算系数矩阵A的行列式值的行列式值D，如果，如果D=0，则输出，则输出错误信息，结束，否则进行第错误信息，结束，否则进行第(3)步；步；（3）对对k=1,2,n，用用b替换替换A的第的第k列数据，并计算替列数据，并计算替换后矩阵的行列式值换后矩阵的行列式值Dk；（4）计算并输出计算并输出x1=D1/D，x2=D2/D，xn=Dn/D，结束。结束。n但但Cramer法则只适用于低阶方程组，高阶方程组工作法则只适用于低阶方程组，高阶方程组工作量太大，故一般用数值方法求解。数值方法分两类：量太大，故一般用数值方法求解。数值方法分两类：u1.直接法直接法u2.迭代法迭代法3.2 Gauss消元法消元法u第二步第二步:回代过程回代过程，解上三角形方程组，得原解上三角形方程组，得原方程组的解。方程组的解。u第一步第一步:消元过程消元过程，将方程组消元化为等价的，将方程组消元化为等价的上三角形方程组上三角形方程组;)1()1()1(2)1(22)1(2211212111nnnnnnnnnnbxabxaxabxaxaxaGauss消元的目的：消元的目的：原始方程组原始方程组约化方程组约化方程组11 11221121 1222221 122nnnnnnnnnna xa xa xba xa xa xba xa xa xb消元过程消元过程(化一般方程组为上三角方程组化一般方程组为上三角方程组)4444343242141343433323213124243232221211414313212111bxaxaxaxabxaxaxaxabxaxaxaxabxaxaxaxa444434241334333231224232221114131211baaaabaaaabaaaabaaaaA)()()()()()()()()()()()()(1414414314213134133132121241231221141312111baaabaaabaaabaaaaA)0(11 a第一轮消元：第一轮消元：u计算计算3个数个数:m21 m31 m41T=a21 a31 a41T/a11u用用-m21乘矩阵第一行后加到矩阵第二行乘矩阵第一行后加到矩阵第二行;u用用-m31乘矩阵第一行后加到矩阵第三行乘矩阵第一行后加到矩阵第三行;u用用-m41乘矩阵第一行后加到矩阵第四行乘矩阵第一行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()()()(2424424323234233121241231221141312112baabaabaaabaaaaA)0()1(22 a计算计算2个数：个数：m32 m42T=a32(1)a42(1)T/a22(1)u 用用-m32乘矩阵第二行后加到矩阵第三行乘矩阵第二行后加到矩阵第三行;u用用-m42乘矩阵第二行后加到矩阵第四行乘矩阵第二行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()()(3434423234233121241231221141312113babaabaaabaaaaA)0()2(33 a第三轮消元：第三轮消元：u计算计算:m43=a43(2)/a33(2)u用用-m43乘矩阵第三行后加到矩阵第四行乘矩阵第三行后加到矩阵第四行;其系数增广矩阵变为：其系数增广矩阵变为：)()()()()()()()()(3443442342343233124124312321221414313212111bxabxaxabxaxaxabxaxaxaxa),1(),1,(),1()1()1()()1()1()()1(nkibmbbnkjiamaankiaamkkikkikikkjikkijkijkkkikik）（记ijijaa)0(回代过程回代过程(解上三角方程组解上三角方程组)a11ann0nnnnnnnnbxabxaxabxaxaxa2222211212111)1,2,1()(1niaxabxabxiinijjijiinnnnn回代过程的计算公式回代过程的计算公式：(n3-n)/3t总工作量：总工作量：S1=n(n-1)/2+(n3-n)/3t总工作量：总工作量：S2=n+n(n-1)/2=n(n+1)/2),1(),1,(),1()1()1()()1()1()()1(nkibmbbnkjiamaankiaamkkikkikikkjikkijkijkkkikik1(),(1,2,1)niijjj inninniiba xbxxinaa=n(n-1)/2+(n3-n)/3+n(n+1)/2n2+(n3-n)/3若用克莱姆法则求解，则工作量为：若用克莱姆法则求解，则工作量为：“”：（n+1个个n阶行列式的值）阶行列式的值）(n+1)(n-1)n!“”：n故总工作量为：故总工作量为：(n+1)(n-1)n!+n如当如当n=6时时,Gauss消元法工作量为消元法工作量为106；而克莱姆法；而克莱姆法则求解工作量为则求解工作量为25206。前前一次课内容回顾一次课内容回顾 约化的主元素约化的主元素ak+1,k+1(k)0(k=0,1,n-1)的充分必要条件是矩阵的充分必要条件是矩阵A的各阶顺序主子式不的各阶顺序主子式不为零。即为零。即01111 kkkkkaaaaD111()1,11/(2,3,.,1)kkkkkaDaDDkn 32563034253432321xxx 32303452536432 202230445.006432 4200445.006432 446245.0432321xxxx1=-13,x2=8,x3=2m21=3/2m31=4/2m32=-3/0.5 1111211nmmL 11111,1,knkkkmmLn将将A 分解为单位下三角矩阵分解为单位下三角矩阵 L 和上三角矩阵和上三角矩阵 U的乘积的算法称为矩阵的乘积的算法称为矩阵A的的三角分解算法三角分解算法。)1(1221 nnnAALLLL)1(1221 nnnbbLLLLLUULLLAn 111211 1112121nnmmmLDi 0（i=1,2,n-1），则，则A可分解为一个单位下三角矩阵可分解为一个单位下三角矩阵L和一个上三角矩阵和一个上三角矩阵U的乘积，且这种分解是唯一的。的乘积，且这种分解是唯一的。2131321231111nnnmLmmmmm由由Gauss消元过程可推得消元过程可推得U即为即为Gauss消元后所得的上三角方程的系数矩阵。消元后所得的上三角方程的系数矩阵。111041221Am21=0，m31=2，m32=-1故故112010001200140111A=LU如果已经有如果已经有 A=L U 则则 AX=b =L U X=b，（1）求解方程组）求解方程组 LY=b 得向量得向量 Y 的值；的值；（L 是下三角矩阵，用顺代算法）是下三角矩阵，用顺代算法）（2）求解方程组）求解方程组 UX=Y 得向量得向量 X 的值。的值。（U 是上三角矩阵，用回代算法）是上三角矩阵，用回代算法）akk(k)太小太小会使误差增大，故应避免采用绝对值小的元素会使误差增大，故应避免采用绝对值小的元素作主元。最好每一步选取系数矩阵中（或消元作主元。最好每一步选取系数矩阵中（或消元后的低阶矩阵中）绝对值最大的元素作主元，后的低阶矩阵中）绝对值最大的元素作主元，以具较好的数值稳定性。以具较好的数值稳定性。3.3 Gauss列主元素消元法列主元素消元法 000.3000.2000.1643.5072.1000.2623.4712.3000.1000.3000.2001.0321xxx0.0012.0003.0001.0001.0003.7124.6232.0002.0001.0725.6433.000 20036006400101002300520040000.1000.3000.2001.0（A/b）=000.2000.5001002300520040000.1000.3000.2001.0 000.1000.3000.2001.0000.2623.4712.3000.1000.3643.5072.1000.2 002.1003.3001.205000.0801.1176.30000.3643.5072.1000.2（A/b）=6870.0868.1005000.0801.1176.30000.3643.5072.1000.2m21=0.5000 m31=-0.0005 m32=0.6300n解为解为nnnnnnnbaaabaaabaaa21222221111211在在A的第一列中选绝对值最大的元素作主元，设该元素的第一列中选绝对值最大的元素作主元，设该元素所在行为第所在行为第i1行，交换第一行与第行，交换第一行与第i1行，进行第一次消行，进行第一次消元；再在第元；再在第2n行的第二列中选绝对值最大的元素作行的第二列中选绝对值最大的元素作主元，设该元素所在行为第主元，设该元素所在行为第i2行，交换第二行与第行，交换第二行与第i2行，行，进行第二次消元，进行第二次消元，直到消元过程完成为止。直到消元过程完成为止。用列主元法解用列主元法解 6745150710623321xxx 6515707104623 6515462370710 5.255.21.661.070710 2.62.65.255.270710列主元矩阵的三角分解：列主元矩阵的三角分解：2)11032(EEE 3)11053(EEE 21EE 5150710623A 51562307101APA 55.261.00710111APFAP 1000010101P 100100131211mmF 0101000012P 10010001322mF3)25.21.03(EEE 61.055.2071011211APFPAPF32EE 2.655.207101122APFPF2121PFPF 1211 FFL若若记记APFFU 12LUPA 可得可得12PPPPn)(1122111nnFPFPFPPL)(nAU LUPA )1()1(max kijnjknikkjiaakk)1(kjikka)1(kjikka)1(kkka3.4 矩阵三角分解法矩阵三角分解法（一）（一）Doolittle分解法分解法（二（二）Crout分解法分解法（三）（三）对称正定矩阵的对称正定矩阵的Cholesky分解法分解法（四）（四）三对角方程组的数值解法三对角方程组的数值解法 nnnnnnaaaaaaaaa212222111211 1112121nnmmm nnnnuuuuuu22211211设设A非奇异，且非奇异，且ALU，L为单位下三角阵，为单位下三角阵，U为上为上三角阵，即三角阵，即一、一、Doolittle分解法：分解法：u11=a11,u1n=a1nm21u11=a21,mn1u11=an1对对A的元素的元素aij，当，当 jk 和和 ik时时 kjkrrjkrkjuuma 11kkikkrrkirikumuma 11m21u12+u22=a22,m21u1n+u2n=a2nm21=a21/u11,mn1=an1/u11 u22=a22-m21u12,u2n=a2n-m21u1n m31u12+m32u22=a32,mn1u12+mn2u22=an2 m32=(a32-m31u12)/u22,mn2=(an2-mn1u12)/u22 11krrjkrkjkjumaukkkrrkirikikuumam/)(11 矩阵矩阵L和矩阵和矩阵U的元素计算公式：的元素计算公式：计算秩序如图所示：计算秩序如图所示：)1,2,1(/)(/)5(),3,2()4(),1,;,3,2(/)()3(),1,;,3,2()2(),3,2(/),2,1()1(111111111111111nniuxuyxuyxniymbybynkkinkuumamnkkjnkumauniuamnjauiinirririinnnnirririikkkrrkirikikkrrjkrkjkjiijjju1ja1 1111uamii 725101391444321131243301024321xxxx14131211uuuuTmmm413121130102T25.05.112423220uuu5.812110 11krrjkrkjkjumau343300uu11211300Tm43100T910044000u4000Tyyyy4321T1611/172010Txxxx4321nnnnuyx iinirririiuxuyx1T4321Tmm423210T11611310kkkrrkirikikuumam1111irririiymby11by 432144434241343332312423222114131211bbbbaaaaaaaaaaaaaaaa 4535251544434241343332312423222114131211aaaaaaaaaaaaaaaaaaaa 4321444342413433323124232221141312111bbbyaaamaaamaaamuuuuk 4321444342413433323124232221141312112bbyyaammaammuuumuuuuk 4321444342413433323124232221141312113byyyammmuummuuumuuuuk 4321444342413433323124232221141312114yyyyummmuummuuumuuuukULy 4535251544434241343332312423222114131211aaaaaaaaaaaaaaaaaaaaA 72510139144432113124330102 725101391444321131243301024321xxxx22100310317111220221334221162913711k 164911621117112113113212021712112310301024kULy 72510139144432113124330102 71391422432215131242310301021k 713911621117112113113212021712112310301023k 44911623112113113212217121123130102yUx解解x 43214321xxxxx 164911621117112113113212021712112310301024kULy 111211221222111nnnnnnuuullllllULAnrriulalrkkrikirir,1,11 nrrjlulaurkrrkjrkrjrj,2,1,/)(11 例例：求矩阵求矩阵A的的Crout分解：分解：303021112A)21()1(222l)21(3032l)23/()21()1(023u)31()23()21(3333l 1003/1102/12/1152/3302/31002A所以所以 52/333/12/312/12/12前前一次课内容回顾一次课内容回顾n若若n阶实矩阵阶实矩阵A为对称正定矩阵，则：为对称正定矩阵，则：u（1）ATA；u（2）对任意的）对任意的X0，有，有XTAX0；u（3）A的各阶顺序主子式大于零。的各阶顺序主子式大于零。nnnnnnnnuuuuuummmA,122112111,121111故故A可进行可进行LU分解：分解：nnnnuuuuuuU222112111/1/11,1,111111122211nnnnnnnuuuuuuuuu1DUTTTLLLDLDLDLDUDLDLDULUA)(21212121121211 ),(2211nnuuudiagD 22211),(nnuuudiag 则则 A=LDU1121211UDDDUU =AT=U1TDLT =U1=LT对对A的第的第i行元素行元素,有有 (j=1,2,i)ijjkjkikall 1jjjkjkikijijlllal/)(11 21112)(ikikiiiilal(j=1,2,i-1)nnnnnnnnnnnnnnaaaaaaaaallllllllllll21222211121122212111212221111111al 对于对于 i=2,3,n 1111112121212121nnnnnllldddlllA nnnnnnnnnnnaaaaaaaaalllddldlddld2122221112112121221121211111jjkjkkikijijdldlal/)(11 11ikikkikiiildladd1=a11,l21=a21/d1,ln1=an1/d1(i=2,3,n；j=1,2,i-1)yxLbLyTyDxLbLyT1 151531140150523121A 81621bd1=a11,l21=a21/d1,ln1=an1/d1jjkjkkikijijdldlal/)(11 11ikikkikiiildlad 13213012100120001L 1000090000100001D 1150134324214144232131331212211ylylylbyylylbyylbyby 1/1/1/1/441331221111442332222443333444xlxlxldyxxlxldyxxldyxdyx四、四、三对角方程组的数值解法三对角方程组的数值解法 nnnnnnnnnffffxxxxbacbacbacb121121111222110011|nniiiabcabcb 且且 ai ci 0，(i=2，3，n-1)其中其中三对角方程组形式如下：三对角方程组形式如下：111112212211nnnnnbabacb )1,3,2(/),3,2(,/,111111nkcnkabacbkkkkkkkkk ),2(,/)(,/1111niyafyfyiiiii nnnyyyxxx212111111 )1,1(,1nixyxyxiiiinn nnnnfffyyy2121221 例：例：求解三对角方程组求解三对角方程组 101653242231124321xxxx 5324223112A 2/536/55/1225/42/512/12 2/535/1222/512L 16/515/412/11U 1016Ly y=3,8/5,4/3,2T 23/45/83Uxx=5,4,3,2Tnnnyyyy1111112211nnnnnnnfbafcbafcbafcb11112222111 nnnnnnnfbafcbafcbay11112222111 1=c1/b1 y1=f1/b1 ni=ci/(bi ai i 1),(i=2,3,n-1)yi=(fi ai yi 1)/(bi aii 1),(i=2,3,n)nnnyyyy1111112211 例：例：用追赶法求上例三对角方程组用追赶法求上例三对角方程组 101653242231124321xxxx 15302421231612其系数增广矩阵为其系数增广矩阵为 15302421231612追赶法求解追赶法求解 153024212313211 1530242422503211主元单位化主元单位化消元消元主元单位化主元单位化 15302425854103211 15302425854103211消元消元 153516251205854103211主元单位化主元单位化 1533465105854103211消元消元 561503465105854103211主元单位化主元单位化 2103465105854103211 561503465105854103211回代可得回代可得 2345x3.5 向量和矩阵的范数向量和矩阵的范数Rn,若若X 的某个的某个非负实值函数非负实值函数N(X)=|X|满足条件：满足条件：u（1）非负性非负性：|X|0且且|X|=0的充要条件的充要条件是是X=0；u（2）齐次性齐次性：|kX|=|k|X|（k R）；）；u（3）三角不等式三角不等式：对：对 XRn,有有|X+Y|X|+|Y|；则称则称N(X)=|X|为为Rn 上的上的向量向量X的范数的范数。3.5.1 向量的范数向量的范数(x1,x2,xn)T Rn,则定义：则定义：222212nxxxX （1）向量的）向量的“2范数范数”：（2）向量的）向量的“范数范数”：inixX 1max（3）向量的）向量的“1范数范数”：（4）向量的）向量的“p范数范数”：niixX11pnipipxX11 421011 X62)1(0122222 X 22,1,0,1max XRn,有有212XnXX XnXX1 XnXX2|211nxxxx 证明（证明（2）:|max|max111knkknkxnxx|1xnxx所以所以0lim)(kkXX*)(limikikxx *lim)(XXkk Rnn,若若A的某个的某个非负实值函非负实值函数数N(A)=|A|满足条件：满足条件：u（1）非负性非负性：|A|0且且|A|=0的充要条件的充要条件是是A=0；u（2）齐次性齐次性：|kA|=|k|A|（k R）；）；u（3）三角不等式三角不等式：对：对 Rnn,有有|A+B|A|+|B|；u（4）相容性相容性：|AB|A|B|；则称则称N(A)=|A|为为Rnn 上的上的矩阵矩阵A的范数的范数。3.5.2 矩阵的范数矩阵的范数定义定义 设向量设向量X Rn，矩阵，矩阵A Rnn，且给定一，且给定一种向量范数种向量范数|X|p，则定义，则定义(p=1,2,)为矩阵为矩阵A的范数，并称为的范数，并称为A的的算子范数算子范数。ppxpXAXA0maxniijnjaA111max)max(2AAATnjijniaA11maxRnn,A=(aij)nn ,则则（1）列范数列范数（2）（3）行范数行范数3.5.3 谱半径、谱范数与方阵的谱半径、谱范数与方阵的F范数范数iniA1max)(AA)(Rnn,则则即即A的谱半径不超过的谱半径不超过A的任何一种算子范数的任何一种算子范数。)(2AA21121njijniFaA为为A的的Frobenius范数，简称范数，简称F范数范数。3.6 误差分析误差分析 220001.111121xx 0001.220001.111121xxbbAAxx 1AAAAxxx 1A-1A-1 A A-1 A(27.0000 -192.0000 210.0000)T(30.0000 -210.0000 228.0000)T 1211111131211211nnnnnHn 5/14/13/14/13/12/13/12/11A 321b 1.321bb AxbrbAx ,*brAcondxxx )(*rAxx1*rAxx1*xAAxb brAcondbrAAxxx)(*1 本章作业