实验报告哈尔滨工程大学教务处制实验四 解线性方程组一．解线性方程组的基本思想 1．直接三角分解法：将系数矩阵A 转变成等价两个矩阵L 和U 的乘积 ，其中L 和U 分别是下三角和上三角矩阵。

当A 的所有顺序主子式都不为0时，矩阵A 可以分解为A=LU ，且分解唯一。

其中L 是单位下三角矩阵，U 是上三角矩阵。

2．平方根法：如果矩阵A 为n 阶对称正定矩阵，则存在一个对角元素为正数的下三角实矩阵L ，使得：A=LL^T 。

当限定L 的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky ）分解。

3．追赶法：设系数矩阵为三对角矩阵1122233111000000000000000n n n nn b c a b c a b A a b c a b ---⎛⎫ ⎪ ⎪ ⎪=⎪ ⎪⎪⎪ ⎪⎝⎭则方程组Ax=f 称为三对角方程组。

设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记1122233110000100000001000000100,00000000000001n n nn b L U γαβγββγβ--⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪ ⎪ ⎪∂==⎪⎪ ⎪⎪ ⎪ ⎪⎪ ⎪⎪⎪∂⎝⎭⎝⎭ 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三角方程组Ux=y 。

4．雅克比迭代法：首先将方程组中的系数矩阵A 分解成三部分，即：A = L+D+U ，如图1所示，其中D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。

之后确定迭代格式，X )1(+k = BX )(k +f ，如图2所示，其中B 称为迭代矩阵，雅克比迭代法中一般记为J 。

（k = 0,1，......）再选取初始迭代向量X )0(，开始逐次迭代。

5．超松弛迭代法（SOR ）它是在GS 法基础上为提高收敛速度，采用加权平均而得到的新算法。

选取分裂矩阵M 为带参数的下三角矩阵M ＝ω1（D －L ω）， 其中ω＞0 为可选择的松弛因子，一般当1<ω<2时称为超松弛。

二.实验题目及实验目的1．(第五章习题8)用直接三角分解（杜利特尔（Doolittle ）分解）求线性方程组141x +251x +361x = 9， 131x +241x +351x = 8，121x + 2x +32x = 8 的解。

2．(第五章习题9)用追赶法解三对角方程组Ax=b ，其中A=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------2100012100012100012100012，b=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛00001. 3．（第五章习题10）用改进的平方根法解线性方程组⎪⎪⎪⎭⎫ ⎝⎛---131321112⎪⎪⎪⎭⎫ ⎝⎛321x x x = ⎪⎪⎪⎭⎫ ⎝⎛654 4．（第六章习题7）用SOR 方法解线性方程组（分别取松弛因子ω=1.03，ω=1，ω=1.1）41x - 2x = 1， -1x +42x - 3x = 4，-2x +43x = -3.精确解x *=（21，1，-21）T .要求当)(*k x x -∞<5×106-时迭代终止，并且对每一个ω值确定迭代次数.5.（第六章习题8）用SOR 方法解线性方程组（取ω=0.9） 51x -22x + 3x = -12， -1x +42x - 23x = 20， 21x -32x +103x = 3.要求当)()1(k k x x -+∞<104-时迭代终止.6．（第六章习题9）设有线性方程组Ax=b ，其中A 为对称正定阵，迭代公式)()1(k k x x =++ω(b- A )(k x )，k=0，1，2…，试证明当0<ω<β2时上述迭代法收敛（其中0<α≤λ(A)≤β）. 7．（第六章计算实习题1）给出线性方程组H n x=b ，其中系数矩阵H n 为希尔伯特矩阵：H n x=(h ij )∈R n n ⨯, h ij =11-+j i ，i ，j=1，2，…，n.假设x *=（1，1，…，1）T ∈R n ，b= H n x *.若取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:1h=s(i,m+1);for j=m:-1:1;if j~=ih=h-x(j)*s(i,j);endendx(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ÊǶԽÇÏßÉÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö%cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö%fÊdz£ÊýÏîbN=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:Na1(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:1x(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ΪÁÐά£»% LDL'·Ö½â£»d(1)=A(1,1);for i=2:nfor j=1:i-1sum1=0;for k=1:j-1sum1=sum1+t(i,k)*l(j,k);endt(i,j)=A(i,j)-sum1;l(i,j)=t(i,j)/d(j);endsum2=0;for k=1:i-1sum2=sum2+t(i,k)*l(i,k);endd(i)=A(i,i)-sum2;endfor i=1:nl(i,i)=1;enddisp('µ¥Î»ÏÂÈý½Ç¾ØÕóLΪ£º'); %½â³öµ¥Î»ÏÂÈý½Ç¾ØÕóL£»ldisp('¶Ô½Ç¾ØÕóDΪ£º'); %½â³ö¶Ô½Ç¾ØÕóD£»d%ÓÉLDL'x=bÇó½âx£»%ÓÉLy=b£¬Çóy£»%ÓÉL'x=inv£¨D£©y£¬Çó½âx£»y(1)=b(1);for i=2:nsum3=0;for k=1:i-1sum3=sum3+l(i,k)*y(k);endy(i)=b(i)-sum3;endx(n)=y(n)/d(n);for i=n-1:-1:1sum4=0;for k=i+1:nsum4=sum4+l(k,i)*x(k);endx(i)=(y(i)/d(i))-sum4;enddisp('ÓÉLy=bÇó½âyµÃ£º');ydisp('Ax=bµÄ½â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('²ÎÊý·¶Î§´íÎó');return;endeps=5.0e-6;D=diag(diag(A)); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏÈý½ÇÕó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)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=200)disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡');return;endenddisp('Ax=bµÄ½âΪ£º');xdisp('µü´ú´ÎÊýΪ£º');n②控制台输入代码：>> 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)5．①SOR方法（文件SOR_2.m）function SOR_2(A,b,x0,w,eps)if(w<=0 || w>=2)error('²ÎÊý·¶Î§´íÎó');return;endD=diag(diag(A)); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏÈý½ÇÕó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)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=200)disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡');return;endenddisp('Ax=bµÄ½âΪ£º');xdisp('µü´ú´ÎÊýΪ£º');n②控制台输入代码：>> A=[5 2 1;-1 4 2;2 -3 10];>> b=[-12;20;3];>> x0=[0;0;0];>> w=0.9;>> eps=10e-4;>> SOR_2(A,b,x0,w,eps)6．此题为证明题，无程序代码。