數(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)成績預(yù)習(xí)部分實(shí)驗(yàn)過程表現(xiàn)實(shí)驗(yàn)報(bào)告部分總成績教師簽字日期哈爾濱工程大學(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)1<<2時(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.0

4、3，=1，=1.1）4 - = 1，- +4- = 4，- +4= -3.精確解x=（，1，-）.要求當(dāng)<5×10時(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)）.7（第六章計(jì)算實(shí)習(xí)題1）給出線性方程組Hx=b，其中系數(shù)矩陣H為希爾伯特矩陣：Hx=(h)R, h=，

5、i，j=1，2，n.假設(shè)x=（1，1，1）R，b= Hx.若取n=6，8，10，分別雅克比迭代法及SOR迭代（=1，1.25，1.5）求解.比較計(jì)算結(jié)果.三實(shí)驗(yàn)手段：指操作環(huán)境和平臺(tái):win7系統(tǒng)下MATLAB R2009a程序語言：一種類似C語言的程序語言，但比C語言要寬松得多，非常方便。 四.程序1. 直接三角分解（文件ZJsanjiao.m）function x=ZJsanjiao(A,b)m,n=size(A);l u=lu(A);s=inv(l)*A,b;x=ones(m,1);for i=m:-1:1 h=s(i,m+1); for j=m:-1:1; if j=i h=h-x(j

6、)*s(i,j); end end x(i)=h/s(i,i);end控制臺(tái)輸入代碼：>> 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=ZG_SDJ(a,b,c,f)%aÊǶԽÇÏßÔªËØ%bÊǶԽÇÏß

7、1;Ï·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö%cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±Èa

8、ÉÙÒ»¸ö%fÊdz£ÊýÏîb 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);

9、y(j1)=temp1/a1(j1);endj1=N;x(j1)=y(j1);for j1=N-1:-1:1 x(j1)=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ΪÁÐά

10、63;»% LDL'·Ö½â£»d(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; 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('µ¥Î

11、;»ÏÂÈý½Ç¾ØÕóLΪ£º'); %½â³öµ¥Î»ÏÂÈý½Ç¾ØÕóL£»ldisp('¶Ô½Ç¾ØÕóDΪ£

12、º'); %½â³ö¶Ô½Ç¾ØÕóD£»d%ÓÉLDL'x=bÇó½âx£»%ÓÉLy=b£¬Çóy£»%ÓÉL'x=inv£¨D£©y£¬Çó½

13、6;x£»y(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(i)/d(i)-sum4;enddisp('ÓÉLy=bÇó½âyµÃ£º');ydisp('Ax=b&#

14、181;ĽâxΪ£º');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_aΪ¾«È·½âif(w<=0 | w>=2) error('²ÎÊý·¶&

15、#206;§´íÎó'); return;endeps=5.0e-6; D=diag(diag(A); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏ

16、0;ý½ÇÕóB=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_a-x)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200) disp('Warning: µü´ú´ÎÊýÌ«¶à£

17、¬¿ÉÄܲ»ÊÕÁ²£¡'); return; endend disp('Ax=bµÄ½âΪ£º');x disp('µü´ú´ÎÊýΪ£º');n控制臺(tái)輸入代碼：>> A=4 -1 0;-1 4 -1;0 -1

18、 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(w<=0 | w>=2) error('²ÎÊý·&

19、#182;Χ´íÎó'); return;end D=diag(diag(A); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏÈ

20、53;½ÇÕóB=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: µü´ú´ÎÊýÌ«¶à£¬

21、¿ÉÄܲ»ÊÕÁ²£¡'); return; endend disp('Ax=bµÄ½âΪ£º');x disp('µü´ú´ÎÊýΪ£º');n控制臺(tái)輸入代碼：>> A=5 2 1;-1 4 2;2 -3 10;>

22、> b=-12;20;3;>> x0=0;0;0;>> w=0.9;>> eps=10e-4;>> SOR_2(A,b,x0,w,eps)6此題為證明題，無程序代碼。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(

23、i)/A(i,i); end for i=1:n sum=sum+(y(i)-x(i)2; end if sqrt(sum)<eps break; else for i=1:n x(i)=y(i); end endend SOR方法（文件SOR_3.m）function SOR_3(n,w)%x_aΪ¾«È·½âif(w<=0 | w>=2) error('²ÎÊý·¶Î§´íÎ

24、ó'); 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); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄ

25、ÉÏÈý½ÇÕóB=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>=2000) disp('Warning: µü´ú´ÎÊýÌ«&#

26、182;࣬¿ÉÄܲ»ÊÕÁ²£¡'); 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-A，其特征值=1-