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

实 验 报 告

课程名称 数值分析 解线性方程组 上机 张振 理学楼407 预习部分 实验过程 表现 实验学时 学号 &实验项目名称 ' 4 30 沈艳 总成绩 实验类型 班级 姓名 实验室名称 指导教师 实验时间 [ 实验成绩 教师签字 实验报告 部分 《 日期

哈尔滨工程大学教务处 制

实验四 解线性方程组

一．解线性方程组的基本思想 1．直接三角分解法：

将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且分解唯一。其中L是单位下三角矩阵，U是上三角矩阵。

.

2．平方根法：

如果矩阵A为n阶对称正定矩阵，则存在一个对角元素为正数的下三角实矩阵L，使得：A=LL^T。当限定L的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky）分解。 3．追赶法：

设系数矩阵为三对角矩阵

b1a20A00c1b2a3000c2b300000000an1bn10an000 cn1bn则方程组Ax=f称为三对角方程组。

设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记

b120 L000000002300300n1n000,0n110012001U0000000000000 n11可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角

方程组Ux=y。

4．雅克比迭代法：

】

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

之后确定迭代格式，X(k1) = 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）分解）求线性方程组

111x1 +x2 +x3= 9， 456》 111x1 +x2 +x3= 8， 3451x1 + x2 +2x3= 8 2的解。

2．(第五章习题9)用追赶法解三对角方程组Ax=b，其中

210001121000A=01210，b=0. 0012100001203．（第五章习题10）用改进的平方根法解线性方程组

211x1123x2 = 131x345 64．（第六章习题7）用SOR方法解线性方程组（分别取松弛因子ω=，ω=1，ω=）

4x1 - x2 = 1， -x1 +4x2- x3= 4，

—

-x2 +4x3= -3.

11，1，-）T.要求当x*x(k)<5×106时迭代终止，并且对每22一个ω值确定迭代次数.

5.（第六章习题8）用SOR方法解线性方程组（取ω=） 精确解x*=（

5x1 -2x2+ x3= -12， -x1 +4x2- 2x3= 20， 2x1 -3x2+10x3= 3.

要求当x(k1)x(k)<104时迭代终止.

6．（第六章习题9）设有线性方程组Ax=b，其中A为对称正定阵，迭代公式

x(k1)x(k)+ω(b- Ax(k))，k=0，1，2…，

试证明当0<ω<

:

2时上述迭代法收敛（其中0<(A)）. 

7．（第六章计算实习题1）给出线性方程组Hnx=b，其中系数矩阵Hn为希尔伯特矩阵：

Hnx=(hij)Rnn, hij=

1，i，j=1，2，…，n.

ij1假设x*=（1，1，…，1）TRn，b= Hnx*.若取n=6，8，10，分别雅克比迭代法及SOR迭代（ω=1，，）求解.比较计算结果. 三．实验手段：

指操作环境和平台:win7系统下MATLAB R2009a

程序语言：一种类似C语言的程序语言，但比C语言要宽松得多，非常方便。 四.程序

1.

①直接三角分解（文件） 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)*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.

①追赶法（文件） :

function x=ZG_SDJ(a,b,c,f) %aÊǶԽÇÏßÔªËØ

%bÊǶԽÇÏßÉÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö %cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö %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); y(j1)=temp1/a1(j1); end j1=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.

①改进的平方根法（文件） function GJPFG(A,b) n=length(b);% nΪÁÐά£» ￥

% 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; end for i=1:n l(i,i)=1; end

·

disp('µ¥Î»ÏÂÈý½Ç¾ØÕóLΪ£º'); %½â³öµ¥Î»ÏÂÈý½Ç¾ØÕóL£» l

disp('¶Ô½Ç¾ØÕóDΪ£º'); %½â³ö¶Ô½Ç¾ØÕóD£» d

%ÓÉLDL'x=bÇó½âx£» %ÓÉLy=b£¬Çóy£»

%ÓÉL'x=inv£¨D£©y£¬Çó½â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; end x(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; end

disp('ÓÉLy=bÇó½âyµÃ£º'); y

disp('Ax=bµÄ½âxΪ£º'); x

②控制台输入代码： >> A=[2 -1 1;-1 -2 3;1 3 1]; >> b=[4;5;6]; `

>> GJPFG(A,b)

4.

①SOR方法（文件）

function SOR_1(A,b,x0,x_a,w) %x_aΪ¾«È·½â if(w<=0 || w>=2) error('²ÎÊý·¶Î§´íÎó'); return; end

:

eps=;

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)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡'); return; end end

disp('Ax=bµÄ½âΪ£º');

—

x

disp('µü´ú´ÎÊýΪ£º'); n

②控制台输入代码： >> A=[4 -1 0;-1 4 -1;0 -1 4]; >> b=[1;4;-3]; >> x0=[0;0;0]; >> x_a=[;1;]; >> w=;

>> SOR_1(A,b,x0,x_a,w) 、

>> w=1;

>> SOR_1(A,b,x0,x_a,w)

>> w=;

>> SOR_1(A,b,x0,x_a,w) 5．

①SOR方法（文件） function SOR_2(A,b,x0,w,eps) if(w<=0 || w>=2) error('²ÎÊý·¶Î§´íÎó'); —

return; end

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-x0)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡'); return; end end

}

disp('Ax=bµÄ½âΪ£º'); x

disp('µü´ú´ÎÊýΪ£º'); n

②控制台输入代码： >> A=[5 2 1;-1 4 2;2 -3 10]; >> b=[-12;20;3]; >> x0=[0;0;0]; >> w=;

>> eps=10e-4; —

>> SOR_2(A,b,x0,w,eps)

6．此题为证明题，无程序代码。 7.

① 雅克比迭代法（文件） 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(sum)`

for i=1:n x(i)=y(i); end end end

② SOR方法（文件） function SOR_3(n,w) %x_aΪ¾«È·½â if(w<=0 || 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)); %Çó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)>=eps x0=x;

'

x =B*x0+f; n=n+1; if(n>=2000)

disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡'); return; end end

disp('Hx=bµÄ½âΪ£º'); x

disp('µü´ú´ÎÊýΪ£º');

）

n

③控制台输入代码： >> x=Jocabi(6) >> x=Jocabi(8) >> x=Jocabi(10) >> SOR_3(6,1) >> SOR_3(6, >> SOR_3(6, >> SOR_3(8,1) >> SOR_3(8, [

>> SOR_3(8,

>> SOR_3(10,1)

>> SOR_3(10, >> SOR_3(10,

五．实验结果比较与分析 1.

2.

￥

3.

4.

5.

] 9.证：

x(k1)(IA)x(k)+b，(k=0，1，2…)

故迭代矩阵B=I-A，其特征值=1-(A). 由||<1，|1-(A)|<1得 0<<

2(A) 故当0<<

2时，更有0<<2(A)，从而有||<1，7.

①雅克比迭代法：

(B)<1，迭代格式收敛。

可以看到用雅克比迭代法求希尔伯特阵方程组的解是病态的，这是因为希尔伯特阵的谱半径大于1，并不收敛。 #

②SOR迭代法：

n=6 ,<104时，

松弛因子 迭代次数 近似解 1 620 （ ，，，，，） 588 （ ，，，，，） 539 （ ，，，，，） n=8 ,<104时

松弛因子 迭代次数 近似解 1 426 （，，，，，，， ） 1157 （，，，，，，，） 1701 （，，，，，，，） n=10 ,<104时

松弛因子 迭代次数 近似解 1 1216 （，，，，，，，，，） 1379 （，，，，，，，，，） 1520 （，，，，，，，，，） 六.学习心得

直接法可以求得线性方程组的精确解，但实际计算中会有误差，一般求得近似解。适用于解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组。

迭代法是逐步逼近求方程组的近似解，但存在收敛性及收敛速度问题，适用于解大型稀疏矩阵方程组。