WordWord資料追趕法求解三對角線性方程組實驗?zāi)康睦镁幊谭椒▽崿F(xiàn)追趕法求解三對角線性方程組實驗內(nèi)容1、學(xué)習(xí)和理解追趕法求解三對角線性方程組的原理及方法2、利用MATLAB編程實現(xiàn)追趕法；3、舉例進(jìn)行求解，并對結(jié)果進(jìn)行分。實驗原理設(shè)n設(shè)n元線性方程組Ax=d的系數(shù)矩陣A為非奇異的三對角矩陣alalclb2a2A=?…c2??????…a(n-1)c(n-1)b(n)a(n)這種方程組稱為三對角線性方程組。顯然，A是上下半寬帶都是1的帶狀矩陣。設(shè)A的前n-1個順序主子式都不為零，根據(jù)定理2.5的推論，A有唯一的Crout分解，并且是保留帶寬的。xx1(n-1)m(n)1(n)11m212A=LU=?…u11u2??????1u(n-1)1其中L是下三角矩陣，U是單位上三角矩陣。利用矩陣相乘法，可以得到：得到：11二al11*ul二clul二cl/11mi二bi(2<i<n)a(i+1)二m(i+l)ui+1(i+1)(l<i<n一l)ci二li?uiui二ci/li1(i+1)二a(i+1)一b(i+l)ui由上列各式可以得到L和U。引入中間量y,令y=Ux，則有：Ax=LUx=Ly=dLy=d已知L和d,可求得y?！鰕l■y2一dl一d2X???——???y(n-1)d(n-1)_y(n)__d(n)_1(n)1l1(n)m212???????…1(n-1)m(n)則可得到y(tǒng)的求解表達(dá)式：yl=d1/11i=2,3,…,nm(i)y(i一1)+1i*y(i)=diy(i)=[di一m(i)y(i一1)]/1i由y=Ux得:1u1x1■y1一1u2x2y2??????X???——???1u(n一1)x(n一1)y(n-1)1x(n)_y(n)_可得到X的求解表達(dá)式：x(n)=y(n)i=n一1,n一2,???,1x(i)=y(i)一u(i)x(i+1)從而得到Ax=d的解X。Matlab編程根據(jù)以上的實驗原理，在Matlab中編程如下函數(shù)X=Trid(A,d)TOC\o"1-5"\h\z；functionx=Trid(A,d)；；%追趕法求解三對角的線性方程組Ax=dI'%b為主對角線元素，a,c分別為次對角線元素，d為右端項:I:%A=[alclI-■%b2a2c2；I-;%:I-%b(n-l)a(n-l)c(n-l)I；%b(n)a(n)]；'%a=[a1???a(n)]%把系數(shù)矩陣的三對角轉(zhuǎn)變成3個列向量7I-:%b=[0b2???b(n)]%不足的元素用0代替I-i%c=[c1.?.c(n-1)0]:n二size(A,l);%n為系數(shù)矩陣的行數(shù)Iia⑴=A(1,1)1b⑴=0:c⑴=A(1,2)!fori=2:n-lIa(i)=A(i,i)I；b(i)=A(i,i-1);ic(i)=A(i,i+1);I；end!a(n)=A(n,n)ib(n)=A(n,n-1)Iic(n)=0IIIl(1)=a⑴；％開始求解L,Uim⑴=0Iifori=2:nim(i)=b(i)%求得m(i)|u(i-l)=c(i-l)/l(i-l);%求得u(i)；l(i)=a(i)-b(i)*u(i-l);%求得l(i)；endI；u(n)=0III；y(l)=d⑴/l⑴；fori=2:nI'y(i)=[d(i)-m(i)*y(i-1)]/l(i);%求得y(i)I;endI-x(n)=y(n);I；fori=n-1：-1：1x(i)=y(i)-u(i)*x(i+l);%求得x(i)I；end]x=x'%將x轉(zhuǎn)置，變?yōu)榱邢蛄縒ord資料WordWord資料在Matlab中新建Trid.m,其中程序為如上虛線框內(nèi)代碼，放在工作目錄。在CommandWindow輸入以下語句：l'■~~~~~~■■iclearall;clc;【fprintf('輸入非奇異三對角系數(shù)矩陣A\n');；A=input('Amatrix二')；%輸入系數(shù)矩陣[fprintf('系數(shù)矩陣');A|遷det(A)==O%判斷系數(shù)矩陣是否奇異ifprintf('系數(shù)矩陣A奇異?。?！請重新輸入！\n');；fprintf('重新輸入非奇異三對角系數(shù)矩陣A\n');;A=input('Amatrix=');'fprintf('系數(shù)矩陣');AI:end;fprintf('輸入矩陣d\n');I；d=input('dmatrix=');:fprintf('矩陣d');d'x=Trid(A,d)%調(diào)用Trid.m中的Trid函數(shù)進(jìn)行求解III■fid=fopen('Ax=d.txt','wt');%生成Ax=d.txt文件；fprintf(fid,'%s\r\n','利用三對角線追趕法求解Ax=d');:fprintf(fid,'%s\r\n','=================');:fori=l:size(A,l)Iifprintf(fid,'%.1f\t',A(i,:));%輸出Ax=d，以上二為分隔符；fprintf(fid,'%s','x');:fprintf(fid,'%u\t',i);ifprintf(fid,'%.1f\n',d(i));Iiend；fprintf(fid,'%s\r\n','=================');I

ifprintf(fid,'%s\r\n','求解得到結(jié)果如下:');I；fori=l:size(A,l)|fprintf(fid,'%s','x');%輸出解向量x(i)ifprintf(fid,'%u',i);I;fprintf(fid,'%s\t','二')；I-fprintf(fid,'%?5f\r\n',x(i));IiendI:fclose(fid)II五舉例計算及分析以課本(《數(shù)值分析》第4版，顏慶津，北京航空航天大學(xué)出版社)27頁例3為例進(jìn)行計算，輸入系數(shù)矩陣A和d：「41-_1_1410.5A=141d=-1141，3142調(diào)用x=Trid(A,d)后，并生成Ax=d.txt文件，CommandWindow同時也會輸出解為0.20.2-同時也會輸出解為0.20.2-0.50.80.3，與課本答案致;Ax=d?txt文件內(nèi)容如下圖：rai?xrai?x文件編輯⑥格式查看?幫助?利用三對角線追趕法求解服"4.91.04.91.00.90.90.91..00..00..00..0X11.0It..01..00..00..0X20.51..0h.01..00..0X3-1■0..01..0h.01..0xU3.00..00..01..0h.0x52.00-2盹盹-0.50080氛0-2盹盹-0.50080氛3盹盹12A=2-2-13-1d=120.5，131__3_求解得到結(jié)果如下:x1=x2=x3=x4=x5=Ln1,Col1再次對程序進(jìn)行驗證，輸入矩陣如下：）利用三對角線追趕法求解抽"r==========================1.02.00.00.0X12.0c2.01.03.00.0X2-1.00.01.02.00.5x31.0-0.00.03.01.0x43.0求解得到結(jié)果如下:x1=1.55556x2=0.22222x3=-1.44444x4=7.33333利用matlab計算Ax得到Ax=d,驗證所得到的x即為方程組的解。由以上兩組計算可表明，該程序能滿足追趕法解三對角線性方程

組，其中要求系數(shù)矩陣滿足要求，即非奇異，三對角且前(n-1)個順序主子式都不為零。在又換了一組系數(shù)矩陣后，系數(shù)矩陣如下：A=[11000;22100;02310;00241;00025]_11-丁2211A=231d=1即241，1_25_1運行程序后出錯，查找問題發(fā)現(xiàn)A進(jìn)行Crout分解時,l1_l1__m2l21A=LU=??????X…l(n-1)m(n)l(n)u11u2??????1u(n-1)1l1=1m2=2u1=112二a2-m2*ul=2-2X1=0再計算u2=c2/12時無法計算，分析發(fā)現(xiàn)A的2階順序主子式為零，無法進(jìn)行C