解线性方程组的基本迭代法实验(ca参考).doc

Lab．解线性方程组的基本迭代法实验
【实验目的和要求】
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．Matlab 6.0及以上版本。
实验讲评：

实验成绩：
评阅教师：
年 月 日

解线性方程组的基本迭代法实验报告
一、算法描述

得到

则有：
第一步
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);
if nargin<3
error('输入错误，最少要输入三个参数');
return;
end
if nargin==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;
if max(abs(eig(B)))>=1
disp('谱半径大于等于1，迭代不收敛，无法进行');
return;
end
n=1;
for i=1:1:n
if sum(A(i,i)~=n)~=n
error('输入的A矩阵的对角线元素不能为0');
return;
end
end
y=B*x0+f;
while norm(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);
if nargin<3%针对这个nargin我还有一个疑问，过一段时间在来处理他!
error('输入错误，最少要输入三个参数');
return;
end
if nargin==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;
if max(abs(eig(B)))>=1
disp('谱半径大于等于1，迭代不收敛，无法进行');
return;
end
n=1;
for i=1:1:n
if sum(A(i,i)~=n)~=n
error('输入的A矩阵的对角线元素不能为0');
return;
end
end
y=B*x0+f;
while norm(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);
if nargin<3%针对这个nargin我还有一个疑问，过一段时间在来处理他!
error('输入错误，最少要输入三个参数');
return;
end
if nargin==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;
for i=1:1:n
if sum(A(i,i)~=n)~=n
error('输入的A矩阵的对角线元素不能为0');
return;
end
end
y=B*x0+f;
while norm(y-x0)>=eps&n<100
x0=y;
y=B*x0+f;
n=n+1;
end
三、数值计算
1）首先编写如下程序实现输入大矩阵A：
A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j1:1:20
if abs(i-j)==1
A(i,j)=-1/2;
end
end
end
for i=1:1:20
for j=1:1:20
if abs(i-j)==2
A(i,j)=-1/4;
end
end
end
第一次给定初始向量为20行一列的0，
右端面项向量b=20行一列的1
迭代误差要求0.005
Jacobi迭代法求解：
A=zeros(20,20);
for i=1:1:20
A(i,i)3;
end
for i=1:1:20
for j=1:1:20
if abs(i-j)==1
A(i,j)=-1/2;
end
end
end
for i=1:1:20
for j=1:1:20
if abs(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,bx0,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

>>
在Command Window中输入：
Gauss-Seidel迭代法求解：
在Command Window中输入：
A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j=1:1:20
if abs(i-j)==1
A(i,j)=-1/2;
end
end
end
for i=1:1:20
for j=1:1:20
if abs(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迭代法求解：
在Command Window中输入
A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j=1:1:20
if abs(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
05003
0.5002
0.4999
0.4978
0.4856
0.4146

n =

7

>>
Gauss-Seidel迭代法求解：
在Command Window中输入
A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j=1:1:20
if abs(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
在Command Window中输入

A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j=1:1:20
if abs(i-j)==1
A(i,j)=-1/2;
end
end
end
for i=1:1:20
for j=1:1:20
if abs(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
-15400
-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
在Command Window中输入
A=zeros(20,20);
for i=1:1:20
A(i,i)=3;
end
for i=1:1:20
for j=1:1:20
if abs(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迭代法的原理及其内在含义的认识、了解、掌握，同时也使我体会到了一些数值计算理论的研究规律，由浅入深，由表及里，有特殊到一般。