数值分析实验7.docx
《数值分析实验7.docx》由会员分享，可在线阅读，更多相关《数值分析实验7.docx（19页珍藏版）》请在冰豆网上搜索。
数值分析实验7
Lab07．解线性方程组的基本迭代法实验
【实验目的和要求】
1．使学生深入理解Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法；
2．通过对Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的程序设计，以提高学生程序设计的能力；
3．应用编写的程序解决具体问题，掌握三种基本迭代法的使用，通过结果的分析了解每一种迭代法的特点。
【实验内容】
1．根据Matlab语言特点，描述Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法。
2．编写Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的M文件。
3．给定
为五对角矩阵
（1）选取不同的初始向量
及右端面项向量b，给定迭代误差要求，分别用编写的Jacobi迭代法和Gauss-Seidel迭代法程序求解，观察得到的序列是否收敛？
若收敛，通过迭代次数分析计算结果并得出你的结论。
（2）用编写的SOR迭代法程序，对于
（1）所选取的初始向量
及右端面项向量b进行求解，松驰系数ω取1<ω<2的不同值，在
时停止迭代，通过迭代次数分析计算结果并得出你的结论。
【实验仪器与软件】
1．CPU主频在1GHz以上，内存在128Mb以上的PC；
2．Matlab6.0及以上版本。
实验讲评：
实验成绩：
评阅教师：
2011年7月1日
Lab07．解线性方程组的基本迭代法实验
一、算法描述
得到
则有：
第一步
Jacobi迭代法
令
则称
为雅克比迭代矩阵
由此可得雅克比迭代的迭代格式如下：
第二步
Gauss-Seidel迭代法
令
，则称
为Gauss-Seidel迭代矩阵
由此可得Gauss-Seidel迭代的迭代格式如下：
第三步
SOR迭代法
令
则有：
令
带入
的值可有
称
为SOR迭代矩阵
由此可得SOR迭代的迭代格式如下：
二、算法程序
Jacobi迭代法的M文件：
>>function[y,n]=Jacobi（A,b,x0,eps）
%*************************************************
%函数名称Jacobi雅克比迭代函数
%参数解释A系数矩阵
%b常数项
%x0估计解向量
%eps误差范围
%返回值
%y解向量
%n迭代次数
%函数功能实现线性方程组的Jacobi迭代求解
%*************************************************
n=length（A）;
ifnargin<3
error（'输入错误，最少要输入三个参数'）;
return;
end
ifnargin==3
eps=1e-6;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
M=D;
N=L+U;
B=M\N;
f=M\b;
ifmax（abs（eig（B）））>=1
disp（'谱半径大于等于1，迭代不收敛，无法进行'）;
return;
end
n=1;
fori=1:
1:
n
ifsum（A（i,i）~=n）~=n
error（'输入的A矩阵的对角线元素不能为0'）;
return;
end
end
y=B*x0+f;
whilenorm（y-x0）>=eps&n<100
x0=y;
y=B*x0+f;
n=n+1;
end
Gauss-Seidel迭代法的M文件和
>>function[y,n]=GaussSeidel（A,b,x0,eps）
%*************************************************
%函数名称GaussSeidel高斯赛德尔迭代函数
%参数解释A系数矩阵
%b常数项
%x0估计解向量
%eps误差范围
%返回值
%y解向量
%n迭代次数
%函数功能实现线性方程组的Jacobi迭代求解
%*************************************************
n=length（A）;
ifnargin<3%针对这个nargin我还有一个疑问，过一段时间在来处理他!
error（'输入错误，最少要输入三个参数'）;
return;
end
ifnargin==3
eps=1e-6;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
M=D-L;
N=U;
B=M\N;
f=M\b;
ifmax（abs（eig（B）））>=1
disp（'谱半径大于等于1，迭代不收敛，无法进行'）;
return;
end
n=1;
fori=1:
1:
n
ifsum（A（i,i）~=n）~=n
error（'输入的A矩阵的对角线元素不能为0'）;
return;
end
end
y=B*x0+f;
whilenorm（y-x0）>=eps&n<100
x0=y;
y=B*x0+f;
n=n+1;
end
SOR迭代法的M文件
>>function[y,n]=SOR（A,b,x0,w,eps）
%*************************************************
%函数名称SOR松弛迭代函数
%参数解释A系数矩阵
%w松弛因子
%b常数项
%x0估计解向量
%eps误差范围
%返回值
%y解向量
%n迭代次数
%函数功能实现线性方程组的Jacobi迭代求解
%*************************************************
n=length（A）;
ifnargin<3%针对这个nargin我还有一个疑问，过一段时间在来处理他!
error（'输入错误，最少要输入三个参数'）;
return;
end
ifnargin==3
eps=1e-6;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B1=D\（L+U）;
M=1/w*（D-w*L）;
N=1/w*（（1-w）*D-w*U）;
B=M\N;
f=M\b;
n=1;
fori=1:
1:
n
ifsum（A（i,i）~=n）~=n
error（'输入的A矩阵的对角线元素不能为0'）;
return;
end
end
y=B*x0+f;
whilenorm（y-x0）>=eps&n<100
x0=y;
y=B*x0+f;
n=n+1;end
三、数值计算
1）首先编写如下程序实现输入大矩阵A：
A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==2
A（i,j）=-1/4;
end
end
end
第一次给定初始向量
为20行一列的0，
右端面项向量b=20行一列的1
迭代误差要求0.005
Jacobi迭代法求解：
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==2
A（i,j）=-1/4;
end
end
end
>>b=ones（20,1）;
>>x0=zeros（20,1）;
>>eps=0.005;
>>[y,n]=Jacobi（A,b,x0,eps）
y=
0.4813
0.5729
0.6321
0.6513
0.6600
0.6632
0.6646
0.6651
0.6653
0.6653
0.6653
0.6653
0.6651
0.6646
0.6632
0.6600
0.6513
0.6321
0.5729
0.4813
n=
9
>>
在CommandWindow中输入：
Gauss-Seidel迭代法求解：
在CommandWindow中输入：
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==2
A（i,j）=-1/4;
end
end
end
>>b=ones（20,1）;
>>x0=zeros（20,1）;
>>eps=0.005;
>>[y,n]=GaussSeidel（A,b,x0,eps）
y=
0.4814
0.5732
0.6325
0.6518
0.6606
0.6640
0.6654
0.6660
0.6662
0.6663
0.6663
0.6663
0.6661
0.6656
0.6642
0.6609
0.6521
0.6328
0.5734
0.4816
n=
7
>>
第二次给定初始向量
为20行一列的0
右端面项向量b=20行一列的1.001
迭代误差要求0.005
Jacobi迭代法求解：
在CommandWindow中输入
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
>>
>>b=1.001*ones（20,1）;
>>x0=zeros（20,1）;
>>eps=0.005;
>>[y,n]=Jacobi（A,b,x0,eps）
y=
0.4146
0.4856
0.4978
0.4999
0.5002
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5002
0.4999
0.4978
0.4856
0.4146
n=
7
>>
Gauss-Seidel迭代法求解：
在CommandWindow中输入
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
>>b=1.001*ones（20,1）;
>>x0=zeros（20,1）;
>>eps=0.005;
>>[y,n]=GaussSeidel（A,b,x0,eps）
y=
0.4145
0.4856
0.4978
0.4999
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5003
0.5000
0.4980
0.4858
0.4146
n=5
观察计算结果得到的序列可以看出其是收敛，在较少的迭代次数下即可的到满足误差要求的解。
（2）第一次给定初始向量
为20行一列的0，
右端面项向量b=20行一列的1
迭代误差要求0.00005
松弛因子为1.5
在CommandWindow中输入
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==2
A（i,j）=-1/4;
end
end
end
>>b=ones（20,1）;
>>x0=zeros（20,1）;
>>w=1.5;
>>eps=1e-5;
>>[y,n]=SOR（A,b,x0,w,eps）
y=
1.0e+012*
-0.5082
-0.9690
-1.5400
-2.1738
-2.8767
-3.6356
-4.4375
-5.2635
-6.0901
-6.8885
-7.6243
-8.2578
-8.7437
-9.0319
-9.0675
-8.7940
-8.1452
-7.0831
-5.4598
-3.5651
n=
100
>>
第二次给定初始向量
为20行一列的0，
右端面项向量b=20行一列的1
迭代误差要求0.00005
松弛因子为1.2
在CommandWindow中输入
>>A=zeros（20,20）;
fori=1:
1:
20
A（i,i）=3;
end
fori=1:
1:
20
forj=1:
1:
20
ifabs（i-j）==1
A（i,j）=-1/2;
end
end
end
>>b=ones（20,1）;
>>x0=zeros（20,1）;
>>w=1.2;
>>eps=1e-5;
>>[y,n]=SOR（A,b,x0,w,eps）
y=
0.2792
0.3246
0.3319
0.3331
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3333
0.3334
0.3331
0.3348
0.3246
0.3874
n=
19
通过对迭代次数及其迭代结果的分析，我的得出的结论是松驰系数ω在SOR迭代中起着相当重要的作用，不同的松驰系数ω，可能对迭代结果带来很大的影响，恰当的松驰系数ω可以加速收敛，得到较为良好的迭代结果，而不恰当的松驰系数ω选取，则可能会得导致无法获得理想的结果，甚至还可能影响到迭代的收敛性。
四、总结
通过对本次实验，使我加深了对Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的原理及其内在含义的认识、了解、掌握，同时也使我体会到了一些数值计算理论的研究规律，由浅入深，由表及里，有特殊到一般。
升级会员 