LAB07-解线性方程组的基本迭代法实验.doc

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．Matlab 6.0及以上版本。
实验讲评：

实验成绩：
评阅教师：
200 年 月 日
Lab07．解线性方程组的基本迭代法实验
一、算法描述
1、雅可比迭代法描述如下：
将线性方程组中的系数矩阵分为三部分
设选取M为A的对角元素部分，即选取M=D（对角矩阵），A=D-N
，由得到解Ax=b的雅可比迭代法
其中，称J为解Ax=b的雅可比迭代法的迭代矩阵。
2、高斯-塞德尔迭代法描述如下：
选取分裂矩阵M为A的下三角部分，即选取M=D-L（下三角矩阵），A=M-N，于是由得到解Ax=b的高斯-塞德尔迭代法
，其中
为解Ax=b的高斯-塞德尔迭代法的迭代矩阵。
3、逐次超松弛迭代法描述如下：
选取分裂矩阵M为带参数的下三角矩阵，其中为可选择的松弛因子，于是由可构造一个迭代法，其迭代矩阵为
，从而得到解Ax=b的逐次超松弛迭代法，简称SOR法。解Ax=b的SOR方法为，其中

二、算法程序
1、编写Jacobi迭代法的M文件如下：
function [x,n]=Jacobi(A,b,x0,r)
format long
n=length(A);
D=diag(diag(A));
L=(-1)*tril(A,-1);
U=(-1)*triu(A,1);
B=inv(D)*(L+U);
f=inv(D)*b;
x=B*x0+f;
n=1;
while norm(x-x0)>=r
x0=x;
x=B*x0+f;
n=n+1;
end
2、编写Gauss-Seidel迭代法的M文件如下：
function [x,n]=GaussSeidel(A,b,x0,r)
n=length(A);
D=diag(diag(A));
L=(-1)*tril(A,-1);
U=(-1)*triu(A,1);
B=inv(D-L)*U;
f=inv(D-L)*b;
x=B*x0+f;
n=1;
while norm(x-x0)>=r
x0=x;
x=B*x0+f;
n=n+1;
end
3、编写SOR迭代法的M文件如下：
function [x,n]=SOR(A,b,x0,w,r)
format long
n=length(A);
D=diag(diag(A));
L=(-1)*tril(A,-1);
U=(-1)*triu(A,1);
Lw=inv(D-w*L)*((1-w)*D+w*U);
f=w*inv(D-w*L)*b;
x=Lw*x0+f;
n=1;
while norm(x-x0)>=r
x0=x;
x=Lw*x0+f;
n=n+1;
end
三、数值实验
矩阵A的程序表示如下：
function A=lucius()
A=[ 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
-1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
-1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2-1/4;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3 -1/2;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1/4 -1/2 3];
(1)选取不同的初始向量及右端面项向量b，给定迭代误差要求，分别用编写的Jacobi迭代法和Gauss-Seidel迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。
1、用Jacobi迭代法程序求解：
clear all;
clc;
r=1.0e-6;
x0=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]';
A=lucius();
b=[7 1 7 1 1 1 5 1 1 7 1 7 1 1 3 1 7 1 2 1]';
[x,n]=Jacobi(A,b,x0,r)

改变数值：
clear all;
clc;
r=1.0e-6;
x0=[5 1 5 1 1 1 5 1 1 5 1 1 5 5 5 5 5 5 5 5]';
A=lucius();
b=[8 7 1 7 7 9 4 7 7 3 1 7 8 7 2 7 5 1 7 7]';
[x,n]=Jacobi(A,b,x0,r)

2、用Gauss-Seidel迭代法程序求解：
clear all;
clc;
r=1.0e-6;
x0=[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]';
A=lucius();
b=[1 8 7 9 8 4 0 6 3 9 4 1 3 4 1 8 1 3 1 4 ]';
[x,n]=GaussSeidel(A,b,x0,r)

改变数值：
clear all;
clc;
r=1.0e-6;
x0=[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]';
A=lucius();
b=[1 3 3 1 4 4 0 8 7 1 0 0 8 5 2 5 8 2 1 3 ]';
[x,n]=GaussSeidel(A,b,x0,r)

根据以上结果可知得到的序列是收敛的
(2)用编写的SOR迭代法程序，对于(1)所选取的初始向量及右端面项向量b进行求解，松驰系数ω取1<ω<2的不同值，在时停止迭代，通过迭代次数分析计算结果并得出你的结论。
SOR迭代法：
clear all;
clc;
w=1.2;
r=1.0e-5;
x0=[7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ]';
A=lucius();
b=[2 0 1 2 0 6 1 5 1 3 1 4 2 6 9 7 1 5 2 3 ]';
[x,n]=SOR(A,b,x0,w,r)

改变数值：
clear all;
clc;
w=1.2;
r=1.0e-5;
x0=[8 5 1 2 5 7 2 4 8 3 2 4 5 8 7 3 2 6 5 3]';
A=lucius();
b=[5 2 2 1 2 5 1 9 9 2 0 1 0 4 0 0 1 5 1 3 ]';
[x,n]=SOR(A,b,x0,w,r)

四、总结
通过以上实验了解了三种迭代法的算法结构，掌握了三种基本迭代法的使用，通过对结果的分析了解了每一种迭代法的特点，并且可以得出结论：相对于雅可比迭代法，高斯-塞得尔迭代法加快了收敛速度，而SOR迭代法的收敛速度与松弛因子ω有关。