追趕法求解三對角線性方程組

1、. word 資料追趕法求解三對角線性方程組一實驗?zāi)康睦镁幊谭椒▽崿F(xiàn)追趕法求解三對角線性方程組。二實驗內(nèi)容1、 學(xué)習(xí)和理解追趕法求解三對角線性方程組的原理及方法；2、 利用 matlab 編程實現(xiàn)追趕法；3、 舉例進行求解，并對結(jié)果進行分。三實驗原理設(shè) n 元線性方程組ax=d的系數(shù)矩陣 a 為非奇異的三對角矩陣11222=(1)(n 1)( )( )acbacaa ncb na n這種方程組稱為三對角線性方程組。顯然，a 是上下半寬帶都是1 的帶狀矩陣。設(shè) a 的前 n-1 個順序主子式都不為零， 根據(jù)定理 2.5 的推論，a 有唯一的 crout 分解，并且是保留帶寬的。其中 l 是下三

2、角矩陣， u 是單位上三角矩陣。利用矩陣相乘法，可以1112212(1)1u(n1)( )( )1lumlualul nm nl n. word 資料得到：由上列各式可以得到l 和 u。引入中間量 y，令yux，則有：已知 l 和 d ，可求得 y。則可得到 y 的求解表達式：11/ 12,3,( ) (1)*y( )=( )( ) (1)/ydlinm i y iliidiy idim i y ili1111111/ 1(2)(1)(1)u(1)(11)/(1)(1)(1)lalucuclmibiina im iil iinciliuiuici lil ia ib iui?axluxlydl

3、yd1112222(1)(n 1)(n 1)( )( )(n)(n)lydmlydl nydm nl nyd. word 資料由yux得：111112221u(n 1)(n 1)(n 1)1(n)(n)uxyuxyxyxy可得到 x 的求解表達式：( )( )1,2,1( )( )u( ) (1)x ny ninnx iy ii x i從而得到ax=d的解 x。四matlab編程根據(jù)以上的實驗原理，在matlab中編程如下函數(shù)x=trid（a,d ） function x=trid(a,d) % 追趕法求解三對角的線性方程組 ax=d % b 為主對角線元素， a，c 分別為次對角線元素， d

4、 為右端項% a= a1 c1 % b2 a2 c2 % . % b(n-1) a(n-1) c(n-1) % b(n) a(n) % a=a1.a(n) % 把系數(shù)矩陣的三對角轉(zhuǎn)變成3 個列向量% b=0 b2.b(n) % 不足的元素用 0 代替% c=c1.c(n-1) 0 . word 資料n=size(a,1); % n為系數(shù)矩陣的行數(shù)a(1)=a(1,1) b(1)=0 c(1)=a(1,2) for i=2:n-1 a(i)=a(i,i) b(i)=a(i,i-1); c(i)=a(i,i+1); end a(n)=a(n,n) b(n)=a(n,n-1) c(n)=0 l(1)

5、=a(1); %開始求解 l，u m(1)=0 for i=2:n m(i)=b(i) %求得 m(i) u(i-1)=c(i-1)/l(i-1); %求得 u(i) l(i)=a(i)-b(i)*u(i-1); %求得 l(i) end u(n)=0 y(1)=d(1)/l(1); for i=2:n y(i)=d(i)-m(i)*y(i-1)/l(i); %求得 y(i) end x(n)=y(n); for i=n-1:-1:1 x(i)=y(i)-u(i)*x(i+1); %求得 x(i) end x=x %將 x 轉(zhuǎn)置，變?yōu)榱邢蛄? word 資料在 matlab中新建 trid.m

6、, 其中程序為如上虛線框內(nèi)代碼，放在工作目錄。在 command window輸入以下語句：clear all;clc; fprintf(輸入非奇異三對角系數(shù)矩陣an); a=input(a matrix=); % 輸入系數(shù)矩陣fprintf(系數(shù)矩陣 );a if det(a)=0 % 判斷系數(shù)矩陣是否奇異fprintf(系數(shù)矩陣 a奇異！ ！ ！請重新輸入！ n); fprintf(重新輸入非奇異三對角系數(shù)矩陣an); a=input(a matrix=); fprintf(系數(shù)矩陣 );a end fprintf(輸入矩陣 dn); d=input(d matrix=); fprintf

7、(矩陣 d);d x=trid(a,d) % 調(diào)用 trid.m 中的 trid 函數(shù)進行求解fid = fopen(ax=d.txt, wt); % 生成 ax=d.txt文件fprintf(fid,%srn,利用三對角線追趕法求解ax=d); fprintf(fid,%srn,=); for i=1:size(a,1) fprintf(fid, %.1ft, a(i,:); % 輸出 ax=d，以上 =為分隔符fprintf(fid, %s, x); fprintf(fid, %ut, i); fprintf(fid, %.1fn, d(i); end fprintf(fid,%srn,=

8、); . word 資料五舉例計算及分析以課本（ 數(shù)值分析第4 版，顏慶津，北京航空航天大學(xué)出版社）27 頁例 3 為例進行計算，輸入系數(shù)矩陣a 和 d：41141=4114a1411，10.5=132d調(diào) 用x=trid(a,d) 后， 并 生 成 ax=d.txt文 件 ， command window同時也會輸出解為0.20.2=0.50.80.3x，與課本答案一致 ; ax=d.txtfprintf(fid,%srn,求解得到結(jié)果如下： ); for i=1:size(a,1) fprintf(fid,%s,x); % 輸出解向量 x（i ）fprintf(fid, %u, i); f

9、printf(fid, %st, =); fprintf(fid,%.5frn,x(i); end fclose(fid) . word 資料文件內(nèi)容如下圖：再次對程序進行驗證，輸入矩陣如下：12213=120.531a，21=13d. word 資料利用 matlab計算 ax 得到 ax=d ， 驗證所得到的 x 即為方程組的解。由以上兩組計算可表明，該程序能滿足追趕法解三對角線性方程組，其中要求系數(shù)矩陣滿足要求，即非奇異，三對角且前（n-1 ）個順序主子式都不為零。在又換了一組系數(shù)矩陣后，系數(shù)矩陣如下：a=1 1 0 0 0;2 2 1 0 0;0 2 3 1 0;0 0 2 4 1;0 0 0 2 5 即11221=23124125a，11= 111d運行程序后出錯，查找問題發(fā)現(xiàn)a 進行 crout 分解時，1112212(1)1u(n1)( )( )1lumlualul nm nl n. word 資料l1=1 m2=2 u1=1 l2=a2-m2*u1=2-21=0 再計算 u2=c2/l2時無法計算，分析發(fā)現(xiàn)a 的 2 階順序主子式為零，無法進行crout分解，所以該程序還無法判斷輸入的系數(shù)矩陣的（