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

实 验 报 告

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

哈尔滨工程大学教务处 制

实验四 解线性方程组

一．解线性方程组的基本思想 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?c1b2a3?000?c2?b3??00?000?0000?an?an?1bn?10??0?0?? ??cn?1??bn??则方程组Ax=f称为三对角方程组。

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

?b1???2?0 L?????0??0?000?00??000??2?3?00?3???n?1??n0??0?0??,??0???n???1?10??01?2?001U???????000??000??????0??00?00??

???0?n?1???1??0可先依次求出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

1．(第五章习题8)用直接三角分解（杜利特尔（Doolittle）分解）求线性方程组

111x1 +x2 +x3= 9， 456111x1 +x2 +x3= 8，

4531x1 + 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??4????????1?23x = ???2??5? ?1???6?31????x3???4．（第六章习题7）用SOR方法解线性方程组（分别取松弛因子ω=1.03，ω=1，

ω=1.1）

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

-x2 +4x3= -3.

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

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

精确解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

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，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)*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ê????????a??

%bê???????é?·?μ??a??￡???êy±èaéùò??? %cê?????????·?μ??a??￡???êy±èaéùò??? %fê?3￡êy??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.

①改进的平方根法（文件GJPFG.m） function GJPFG(A,b) n=length(b);% n?aáD??￡? % LDL'·??a￡? 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('μ￥????èy?????óL?a￡o'); %?a3?μ￥????èy?????óL￡? l

disp('???????óD?a￡o'); %?a3????????óD￡? d

%óéLDL'x=b?ó?ax￡? %óéLy=b￡??óy￡?

%óéL'x=inv￡¨D￡?y￡??ó?ax￡? 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?ó?ayμ?￡o'); y

disp('Ax=bμ??ax?a￡o'); 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?a??è·?a

if(w<=0 || w>=2)

error('2?êy·??§′í?ó'); return; end

eps=5.0e-6;

D=diag(diag(A)); %?óAμ????????ó L=-tril(A,-1); %?óAμ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x_a-x)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?');

return; end end

disp('Ax=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); 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('2?êy·??§′í?ó'); return; end

D=diag(diag(A)); %?óAμ????????ó L=-tril(A,-1); %?óAμ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x-x0)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=200)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?'); return; end end

disp('Ax=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); 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．此题为证明题，无程序代码。 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(sum)

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

end

② SOR方法（文件SOR_3.m） function SOR_3(n,w) %x_a?a??è·?a

if(w<=0 || w>=2)

error('2?êy·??§′í?ó'); 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μ???èy???ó U=-triu(A,1); %?óAμ?é?èy???ó B=inv(D-L*w)*((1-w)*D+w*U); f=w*inv((D-L*w))*b; x=B*x0+f;

n=1; %μü′ú′?êy

while norm(x-x0)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=2000)

disp('Warning: μü′ú′?êyì??à￡??é?ü2?ê?á2￡?'); return; end end

disp('Hx=bμ??a?a￡o'); x

disp('μü′ú′?êy?a￡o'); 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.证：

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

①雅克比迭代法：

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

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

n=6 ,?<10?4时，

松弛因子? 迭代次数 近似解 1 620 （1.0005 ，1.0045，0.9626，1.0441，1.0285，0.9583） 1.25 588 （ 0.9997，1.0134，0.9362，1.0706，1.0230，0.9555） 1.5 539 （ 0.9991，1.0211，0.9082，1.1146，0.9899，0.9656） n=8 ,?<10?4时

松弛因子? 迭代次数 近似解 1 426 （0.9970，1.0417，0.8967，1.0155，1.0654，1.0505，0.9991， 0.9309） 1.25 1157 （0.9971，1.0353，0.9174，1.0167，1.0408，1.0378，1.0022，0.9508） 1.5 1701 （0.9980，1.0232，0.9484，1.0045，1.0260，1.0324，1.0037，0.9623） n=10 ,?<10?4时

松弛因子? 迭代次数 近似解 1 1216 （0.9977，1.0203，0.9797，0.9654，1.0010，1.0300，1.0367，1.0223，0.9924，0.9525） 1.25 1379 （0.9985，1.0103，1.0049，0.9503，1.0001，1.0257，1.0347，1.0221，0.9946，0.9570） 1.5 1520 （0.9993，1.0003，1.0354，0.9171，1.0187，1.0116，1.0405，1.0185，0.9982，0.9588） 六.学习心得

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

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

本文来源：https://www.bwwdw.com/article/02v6.html