數(shù)值分析上機(jī)實(shí)驗(yàn)——解線性方程組

1、實(shí) 驗(yàn) 報(bào) 告課程名稱數(shù)值分析實(shí)驗(yàn)項(xiàng)目名稱解線性方程組實(shí)驗(yàn)類型上機(jī)實(shí)驗(yàn)學(xué)時(shí)4班級(jí)20111131學(xué)號(hào)2011113130姓名張振指導(dǎo)教師沈艷實(shí)驗(yàn)室名稱理學(xué)樓407實(shí)驗(yàn)時(shí)間2013.12.9實(shí)驗(yàn)成績(jī)預(yù)習(xí)部分實(shí)驗(yàn)過程表現(xiàn)實(shí)驗(yàn)報(bào)告部分總成績(jī)教師簽字日期哈爾濱工程大學(xué)教務(wù)處 制實(shí)驗(yàn)四 解線性方程組一解線性方程組的基本思想1直接三角分解法：將系數(shù)矩陣A轉(zhuǎn)變成等價(jià)兩個(gè)矩陣L和U的乘積 ，其中L和U分別是下三角和上三角矩陣。當(dāng)A的所有順序主子式都不為0時(shí)，矩陣A可以分解為A=LU，且分解唯一。其中L是單位下三角矩陣，U是上三角矩陣。2平方根法：如果矩陣A為n階對(duì)稱正定矩陣，則存在一個(gè)對(duì)角元素為正數(shù)的下三角實(shí)

2、矩陣L，使得：A=LLT。當(dāng)限定L的對(duì)角元素為正時(shí)，這種分解是唯一的，稱為平方根法（Cholesky）分解。3追趕法：設(shè)系數(shù)矩陣為三對(duì)角矩陣則方程組Ax=f稱為三對(duì)角方程組。設(shè)矩陣A非奇異，A有Crout分解A=LU，其中L為下三角矩陣，U為單位上三角矩陣，記可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程組Ly=f得出y，再求解上三角方程組Ux=y。4雅克比迭代法：首先將方程組中的系數(shù)矩陣A分解成三部分，即：A = L+D+U，如圖1所示，其中D為對(duì)角陣，L為下三角矩陣，U為上三角矩陣。之后確定迭代格式，X = BX +f ，如圖2所示，其中B稱為迭代矩陣，雅克比迭代法中一般記為J

3、。（k = 0,1，.）再選取初始迭代向量X，開始逐次迭代。5超松弛迭代法（SOR）它是在GS法基礎(chǔ)上為提高收斂速度，采用加權(quán)平均而得到的新算法。選取分裂矩陣M為帶參數(shù)的下三角矩陣M（D），其中0 為可選擇的松弛因子，一般當(dāng)12時(shí)稱為超松弛。二.實(shí)驗(yàn)題目及實(shí)驗(yàn)?zāi)康?(第五章習(xí)題8)用直接三角分解（杜利特爾（Doolittle）分解）求線性方程組 + += 9， + += 8， + += 8的解。2(第五章習(xí)題9)用追趕法解三對(duì)角方程組Ax=b，其中A=，b=.3（第五章習(xí)題10）用改進(jìn)的平方根法解線性方程組 = 4（第六章習(xí)題7）用SOR方法解線性方程組（分別取松弛因子=1.03，=1，=1.

4、1）4 - = 1，- +4- = 4，- +4= -3.精確解x=（，1，-）.要求當(dāng)510時(shí)迭代終止，并且對(duì)每一個(gè)值確定迭代次數(shù).5.（第六章習(xí)題8）用SOR方法解線性方程組（取=0.9）5 -2+ = -12，- +4- 2= 20，2 -3+10= 3.要求當(dāng)10時(shí)迭代終止.6（第六章習(xí)題9）設(shè)有線性方程組Ax=b，其中A為對(duì)稱正定陣，迭代公式+(b- A)，k=0，1，2，試證明當(dāng)0時(shí)上述迭代法收斂（其中0 A=1/4,1/5,1/6;1/3,1/4,1/5;1/2,1,2; b=9;8;8; x=ZJsanjiao(A,b)2. 追趕法（文件ZG_SDJ.m）function x=

5、ZG_SDJ(a,b,c,f)%a%ba%ca%fb N=length(a);b=b,0;c=0,c; a1=zeros(N,1);b1=zeros(N,1);y=zeros(N,1);x=zeros(N,1); a1(1)=a(1);b1(1)=b(1)/a1(1);y(1)=f(1)/a1(1);for j1=2:N a1(j1)=a(j1)-c(j1)*b1(j1-1); b1(j1)=b(j1)/a1(j1); temp1=f(j1)-c(j1)*y(j1-1); y(j1)=temp1/a1(j1);endj1=N;x(j1)=y(j1);for j1=N-1:-1:1 x(j1)=

6、y(j1)-b1(j1)*x(j1+1);end控制臺(tái)輸入代碼： a=2 2 2 2 2; b=-1 -1 -1 -1; c=-1 -1 -1 -1; f=1;0;0;0;0; x=ZG_SDJ(a,b,c,f)3.改進(jìn)的平方根法（文件GJPFG.m）function GJPFG(A,b)n=length(b);% n% LDLd(1)=A(1,1);for i=2:n for j=1:i-1 sum1=0; for k=1:j-1 sum1=sum1+t(i,k)*l(j,k); end t(i,j)=A(i,j)-sum1; l(i,j)=t(i,j)/d(j); end sum2=0;

7、for k=1:i-1 sum2=sum2+t(i,k)*l(i,k); end d(i)=A(i,i)-sum2;endfor i=1:n l(i,i)=1;enddisp(L); %Lldisp(D); %Dd%LDLx=bx%Ly=by%Lx=invDyxy(1)=b(1);for i=2:n sum3=0; for k=1:i-1 sum3=sum3+l(i,k)*y(k); end y(i)=b(i)-sum3;endx(n)=y(n)/d(n);for i=n-1:-1:1 sum4=0; for k=i+1:n sum4=sum4+l(k,i)*x(k); end x(i)=(y

8、(i)/d(i)-sum4;enddisp(Ly=by);ydisp(Ax=bx);x控制臺(tái)輸入代碼： A=2 -1 1;-1 -2 3;1 3 1; b=4;5;6; GJPFG(A,b)4. SOR方法（文件SOR_1.m）function SOR_1(A,b,x0,x_a,w)%x_aif(w=2) error(); return;endeps=5.0e-6; D=diag(diag(A); %AL=-tril(A,-1); %AU=-triu(A,1); %AB=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1; % while

9、 norm(x_a-x)=eps x0=x; x =B*x0+f; n=n+1; if(n=200) disp(Warning: ); return; endend disp(Ax=b);x disp();n控制臺(tái)輸入代碼： A=4 -1 0;-1 4 -1;0 -1 4; b=1;4;-3; x0=0;0;0; x_a=0.5;1;-0.5; w=1.03; SOR_1(A,b,x0,x_a,w) w=1; SOR_1(A,b,x0,x_a,w) w=1.1; SOR_1(A,b,x0,x_a,w)5SOR方法（文件SOR_2.m）function SOR_2(A,b,x0,w,eps)if

10、(w=2) error(); return;end D=diag(diag(A); %AL=-tril(A,-1); %AU=-triu(A,1); %AB=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1; % while norm(x-x0)=eps x0=x; x =B*x0+f; n=n+1; if(n=200) disp(Warning: ); return; endend disp(Ax=b);x disp();n控制臺(tái)輸入代碼： A=5 2 1;-1 4 2;2 -3 10; b=-12;20;3; x0=0;0;0; w

11、=0.9; eps=10e-4; SOR_2(A,b,x0,w,eps)6此題為證明題，無(wú)程序代碼。7. 雅克比迭代法（文件Jocabi.m）function x=Jocabi(n)A=hilb(n);x_a=ones(n,1);b=A*x_a;eps=1e-4; n=length(b);N=50;x=zeros(n,1);y=zeros(n,1); for k=1:N sum=0; for i=1:n y(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i); end for i=1:n sum=sum+(y(i)-x(i)2; end if sqrt(su

12、m)eps break; else for i=1:n x(i)=y(i); end endend SOR方法（文件SOR_3.m）function SOR_3(n,w)%x_aif(w=2) error(); return;end x0=zeros(n,1);A=hilb(n);x_a=ones(n,1);b=A*x_a;eps=1e-4; D=diag(diag(A); %AL=-tril(A,-1); %AU=-triu(A,1); %AB=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1; % while norm(x-x0)=

13、eps x0=x; x =B*x0+f; n=n+1; if(n=2000) disp(Warning: ); return; endend disp(Hx=b);x disp();n控制臺(tái)輸入代碼： x=Jocabi(6) x=Jocabi(8) x=Jocabi(10) SOR_3(6,1) SOR_3(6,1.25) SOR_3(6,1.5) SOR_3(8,1) SOR_3(8,1.25) SOR_3(8,1.5) SOR_3(10,1) SOR_3(10,1.25) SOR_3(10,1.5)五 實(shí)驗(yàn)結(jié)果比較與分析1.2.3.4.5.9.證：+b，(k=0，1，2)故迭代矩陣B=I-

14、A，其特征值=1-.由|1，|1-|1得0故當(dāng)0時(shí)，更有0，從而有|1，1，迭代格式收斂。7.雅克比迭代法：可以看到用雅克比迭代法求希爾伯特陣方程組的解是病態(tài)的，這是因?yàn)橄柌仃嚨淖V半徑大于1，并不收斂。SOR迭代法：松弛因子迭代次數(shù)近似解1620（1.0005 ，1.0045，0.9626，1.0441，1.0285，0.9583）1.25588（ 0.9997，1.0134，0.9362，1.0706，1.0230，0.9555）1.5539（ 0.9991，1.0211，0.9082，1.1146，0.9899，0.9656）n=6 ,10時(shí)，n=8 ,10時(shí)松弛因子迭代次數(shù)近似解14

15、26（0.9970，1.0417，0.8967，1.0155，1.0654，1.0505，0.9991， 0.9309）1.251157（0.9971，1.0353，0.9174，1.0167，1.0408，1.0378，1.0022，0.9508）1.51701（0.9980，1.0232，0.9484，1.0045，1.0260，1.0324，1.0037，0.9623）n=10 ,10時(shí)松弛因子迭代次數(shù)近似解11216（0.9977，1.0203，0.9797，0.9654，1.0010，1.0300，1.0367，1.0223，0.9924，0.9525）1.251379（0.9985，1.0103，1.0049，0.9