数值分析上机实验——解线性方程组

1、实 验 报 告课程名称数值分析实验项目名称解线性方程组实验类型上机实验学时4班级20111131学号2011113130姓名张振指导教师沈艳实验室名称理学楼407实验时间2013.12.9实验成绩预习部分实验过程表现实验报告部分总成绩教师签字日期哈尔滨工程大学教务处 制实验四 解线性方程组一解线性方程组的基本思想1直接三角分解法：将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且分解唯一。其中L是单位下三角矩阵，U是上三角矩阵。2平方根法：如果矩阵A为n阶对称正定矩阵，则存在一个对角元素为正数的下三角实

2、矩阵L，使得：A=LLT。当限定L的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky）分解。3追赶法：设系数矩阵为三对角矩阵则方程组Ax=f称为三对角方程组。设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。4雅克比迭代法：首先将方程组中的系数矩阵A分解成三部分，即：A = L+D+U，如图1所示，其中D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，X = BX +f ，如图2所示，其中B称为迭代矩阵，雅克比迭代法中一般记为J

3、。（k = 0,1，.）再选取初始迭代向量X，开始逐次迭代。5超松弛迭代法（SOR）它是在GS法基础上为提高收敛速度，采用加权平均而得到的新算法。选取分裂矩阵M为带参数的下三角矩阵M（D），其中0 为可选择的松弛因子，一般当1<<2时称为超松弛。二.实验题目及实验目的1(第五章习题8)用直接三角分解（杜利特尔（Doolittle）分解）求线性方程组 + += 9， + += 8， + += 8的解。2(第五章习题9)用追赶法解三对角方程组Ax=b，其中A=，b=.3（第五章习题10）用改进的平方根法解线性方程组 = 4（第六章习题7）用SOR方法解线性方程组（分别取松弛因子=1.0

4、3，=1，=1.1）4 - = 1，- +4- = 4，- +4= -3.精确解x=（，1，-）.要求当<5×10时迭代终止，并且对每一个值确定迭代次数.5.（第六章习题8）用SOR方法解线性方程组（取=0.9）5 -2+ = -12，- +4- 2= 20，2 -3+10= 3.要求当<10时迭代终止.6（第六章习题9）设有线性方程组Ax=b，其中A为对称正定阵，迭代公式+(b- A)，k=0，1，2，试证明当0<<时上述迭代法收敛（其中0<(A)）.7（第六章计算实习题1）给出线性方程组Hx=b，其中系数矩阵H为希尔伯特矩阵：Hx=(h)R, h=，

5、i，j=1，2，n.假设x=（1，1，1）R，b= Hx.若取n=6，8，10，分别雅克比迭代法及SOR迭代（=1，1.25，1.5）求解.比较计算结果.三实验手段：指操作环境和平台:win7系统下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控制台输入代码：>> 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控制台输入代码：>> 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.改进的平方根法（文件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控制台输入代码：>> 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控制台输入代码：>> 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控制台输入代码：>> 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控制台输入代码：>> 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)五 实验结果比较与分析1.2.3.4.5.9.证：+b，(k=0，1，2)故迭代矩阵B=I-A，其特征值=1-