实验解线性方程组的基本迭代法实验

1/22数值分析实验报告姓名忘川学号1205025106系别数学系班级12级主讲教师指导教师实验日期2014/6/25专业信息与计算科学专业课程名称数值分析同组实验者无一、实验名称：实验八、解线性方程组的基本迭代法实验二、实验目的：1．深入理解Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法；2．通过对Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的程序设计，提高程序设计的能力；3．应用编写的程序解决具体问题，掌握三种基本迭代法的使用，通过结果的分析了解每一种迭代法的特点。三、实验内容及要求：1．根据Matlab语言特点，描述Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法。2．编写Jacobi迭代法、Gauss-Seidel迭代法和SOR迭代法的M文件。3．给定2020RA为五对角矩阵321412132141412132141412132141412132141213(1)选取不同的初始向量)0(x及右端面项向量b，给定迭代误差要求，分别用编写的Jacobi迭代法和Gauss-Seidel迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。(2)用编写的SOR迭代法程序，对于(1)所选取的初始向量)0(x及右端面项向量b进行求解，松驰系数ω取1ω2的不同值，在5)1()(10kkxx时停止迭代，通过迭代次数分析计算结果并得出你的结论。四、实验步骤（或记录）一、算法描述2/22由121,111112111212,1221222221,11,21,1212,100000000nnnnnnnnnnnnnnnnnnnnaaaaaaaaaaaaaaAaaaaaaaaaa得到1122nnaaDa211,11,212,10000nnnnnnaLaaaaa121,112,121,0000nnnnnnaaaaaUa,则有：ADLU则有：第一步：Jacobi迭代法()()\()\AxbDLUxbDxLUxbxDLUxDbADLU令\()\JDLUfDb则称J为雅克比迭代矩阵由此可得雅克比迭代的迭代格式如下：(0)(1)(),,0,1,2,kkxxJxfk初始向量第二步Gauss-Seidel迭代法()()()\()\AxbDLUxbDLxUxbxDLUxDLbADLU3/22令()\()\GDLUfDLb，则称G为Gauss-Seidel迭代矩阵由此可得Gauss-Seidel迭代的迭代格式如下：(0)(1)(),,0,1,2,kkxxGxfk初始向量第三步SOR迭代法0111(((1)))()((1))wADLUDwLwDwUDwLwDwU令1()MDwLw,1((1))NwDwUw则有：AMN()\\AxbMNxbMxNxbxMNxMbAMN令\\WLMNfMb带入,MN的值可有111((1))()((1))()()1()WwDwULDwLwDwUDwLbfwDwLbDwLw称WL为SOR迭代矩阵，由此可得SOR迭代的迭代格式如下：(0)(1)(),,0,1,2,kkWxxLxfk初始向量二、算法程序Jacobi迭代法的M文件：function[y,n]=Jacobi(A,b,x0,eps)%*************************************************4/22%函数名称Jacobi雅克比迭代函数%参数解释A系数矩阵%b常数项%x0估计解向量%eps误差范围%返回值%y解向量%n迭代次数%函数功能实现线性方程组的Jacobi迭代求解%*************************************************n=length(A);ifnargin3error('输入错误，最少要输入三个参数');return;endifnargin==3eps=1e-6;endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);M=D;N=L+U;B=M\N;f=M\b;5/22ifmax(abs(eig(B)))=1disp('谱半径大于等于1，迭代不收敛，无法进行');return;endn=1;fori=1:1:nifsum(A(i,i)~=n)~=nerror('输入的A矩阵的对角线元素不能为0');return;endendy=B*x0+f;whilenorm(y-x0)=eps&n100x0=y;y=B*x0+f;n=n+1;endGauss-Seidel迭代法的M文件和function[y,n]=GaussSeidel(A,b,x0,eps)%*************************************************%函数名称GaussSeidel高斯赛德尔迭代函数%参数解释A系数矩阵%b常数项%x0估计解向量6/22%eps误差范围%返回值%y解向量%n迭代次数%函数功能实现线性方程组的Jacobi迭代求解%*************************************************n=length(A);ifnargin3%针对这个nargin我还有一个疑问，过一段时间在来处理他!error('输入错误，最少要输入三个参数');return;endifnargin==3eps=1e-6;endD=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)))=1disp('谱半径大于等于1，迭代不收敛，无法进行');return;end7/22n=1;fori=1:1:nifsum(A(i,i)~=n)~=nerror('输入的A矩阵的对角线元素不能为0');return;endendy=B*x0+f;whilenorm(y-x0)=eps&n100x0=y;y=B*x0+f;n=n+1;endSOR迭代法的M文件function[y,n]=SOR(A,b,x0,w,eps)%*************************************************%函数名称SOR松弛迭代函数%参数解释A系数矩阵%w松弛因子%b常数项%x0估计解向量%eps误差范围%返回值%y解向量8/22%n迭代次数%函数功能实现线性方程组的Jacobi迭代求解%*************************************************n=length(A);ifnargin3%针对这个nargin我还有一个疑问，过一段时间在来处理他!error('输入错误，最少要输入三个参数');return;endifnargin==3eps=1e-6;endD=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:nifsum(A(i,i)~=n)~=nerror('输入的A矩阵的对角线元素不能为0');return;end9/22endy=B*x0+f;whilenorm(y-x0)=eps&n100x0=y;y=B*x0+f;n=n+1;end三、数值计算1）首先编写如下程序实现输入大矩阵A：A=zeros(20,20);fori=1:1:20A(i,i)=3;endfori=1:1:20forj=1:1:20ifabs(i-j)==1A(i,j)=-1/2;endendendfori=1:1:20forj=1:1:20ifabs(i-j)==2A(i,j)=-1/4;endend10/22end第一次给定初始向量为20行一列的0，右端面项向量b=20行一列的1迭代误差要求0.005Jacobi迭代法求解：A=zeros(20,20);fori=1:1:20A(i,i)=3;endfori=1:1:20forj=1:1:20ifabs(i-j)==1A(i,j)=-1/2;endendendfori=1:1:20forj=1:1:20ifabs(i-j)==2A(i,j)=-1/4;endendendb=ones(20,1);x0=zeros(20,1);eps=0.005;11/22[y,n]=Jacobi(A,b,x0,eps)y=0.48130.57290.63210.65130.66000.66320.66460.66510.66530.66530.66530.66530.66510.66460.66320.66000.65130.63210.57290.481312/22n=9在CommandWindow中输入：Gauss-Seidel迭代法求解：在CommandWindow中输入：A=zeros(20,20);fori=1:1:20A(i,i)=3;endfori=1:1:20forj=1:1:20ifabs(i-j)==1A(i,j)=-1/2;endendendfori=1:1:20forj=1:1:20ifabs(i-j)==2A(i,j)=-1/4;endendend13/22b=ones(20,1);x0=zeros(20,1);eps=0.005;[y,n]=GaussSeidel(A,b,x0,eps)y=0.48140.57320.63250.65180.66060.66400.66540.66600.66620.66630.66630.66630.66610.66560.66420.66090.65210.63280.573414/220.4816n=7第二次给定初始向量为20行一列的0右端面项向量b=20行一列的1.001迭代误差要求0.005Jacobi迭代法求解：在CommandWindow中输入A=zeros(20,20);fori=1:1:20A(i,i)=3;endfori=1:1:20forj=1:1:20ifabs(i-j)==1A(i,j)=-1/2;endendendb=1.001*ones(20,1);15/22x0=zeros(20,1);eps=0.005;[y,n]=Jacobi(A,b,x0,eps)y=0.41460.48560.49780.49990.50020.50030.50030.50030.50030.50030.50030.50030.50030.50030.50030.50020.49990.49780.48560.414616/22n=7Gauss-Seidel迭代法求解：在CommandWindow中输入A=zeros(20,20);fori=1:1:20A(i,i)=3;endfori=1:1:20forj=1:1:20ifabs(i-j)==1A(i,j)=-1/2;endendendb=1.001*ones(20,1);x0=zeros(20,1);ep