资源描述:
《MATLAB数学实验》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、姚文20063100009106信计数学实验报告一、实验目的(1)复习用Gauss_Seidel,Jacobi,超松弛迭代法求解线性方程组。(2)比较用以上各种方法求解问题所要求满足的条件二、实验问题构造3阶Hilbert矩阵A,并使其精确解x为3维全一向量(111),进而构造出右端向量b。对于方程组Ax=b,从理论上分析Jacobi迭代法和Gauss_Seidel迭代法是否收敛,并求解各自的谱半径和收敛速度。逐步增加n的值,重复求解上述问题,并总结计算结果。三、问题的分析对于三阶Hilbert矩阵,以下分析其用Gauss_Seidel迭代法和用
2、Jacobi迭代法所得的迭代矩阵的谱半径,并由此分析其用相应方法求解的可行性。(1)在MATLAB的commandwindow中输入以下命令:>>formatrat>>%有理数格式>>a=hilb(3)a=11/21/31/21/31/41/31/41/5>>D=diag(diag(a)',0)D=10001/30001/5>>%取矩阵a的对角元素>>L=-tril(a,-1)L=000-1/200-1/3-1/40第12页共12页姚文20063100009106信计>>%取矩阵a的下三角矩阵的负值>>U=-triu(a,1)U=0-1/2-1/
3、300-1/4000>>%取矩阵a的上三角矩阵的负值>>B=inv(D-L)*UB=0-1/2-1/303/4-1/40-5/48125/144>>%矩阵B为矩阵a用Gauss_seidel法求解时得到的的迭代矩阵>>pr=max(abs(eig(B)))pr=2101/2142%迭代矩阵B的谱半径由上述结果可见,用Gauss_Seidel法得到的迭代矩阵谱半径小于一,用此方法求解,结果收敛。(2)在MATLAB的commandwindow中输入以下命令:>>formatrat>>a=hilb(3);>>D=diag(diag(a)',0)D=1
4、0001/30001/5>>%取矩阵a的对角元素>>B=inv(D)*(D-a)B=0-1/2-1/3第12页共12页姚文20063100009106信计-3/20-3/4-5/3-5/40>>%Jacobi方法的迭代矩阵>>pr=max(abs(eig(B)))pr=1449/841>>%迭代矩阵的谱半径由上述结果可见,用Jacobi法得到的迭代矩阵谱半径大于一,用此方法求解,结果将不收敛。下面编写程序,对于输入的正整数w,构造维数为w的Hilbert矩阵,并以全一向量为精确解构造右端向量,然后分别用Gauss_Seidel方法和Jacobi方
5、法求解方程组,并与精确解比较误差。(a)用Gauss_Seidel迭代法分析该问题,用MATLAB编辑程序如下:%Gauss_Seidel迭代法的程序function[res,k,pr,wc]=Gauss_Seidel(err,M)%输出结果为方程组的解res,迭代次数k和系数矩阵的谱半径pr%err为误差控制,M为限制的最大迭代次数formatlongw=input('pleaseinputanumber:');%输入一个正整数w,为hilbert矩阵的维数A=hilb(w);%根据输入的正整数构造hilbert矩阵Am=size(A);%求A
6、的维数n=m(1);xjq=ones(n,1);%方程组的精确解x0=zeros(n,1);b=A*xjq;%方程组的右端向量D=diag(diag(A)',0);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L)*U;g=inv(D-L)*b;pr=max(abs(eig(B)));%求迭代矩阵B的谱半径x1=B*x0+g;k=0;第12页共12页姚文20063100009106信计while(norm((x1-x0),2)>err)k=k+1;x0=x1;x1=B*x0+g;if(k>M)disp('迭代次数过大!')
7、,return;endendk;res=x1;wc=norm(res-xjq,2);%wc为所求结果和精确解之差的二范数(b)用jacobi迭代法分析该问题,用MATLAB编辑程序如下:%jacobi迭代法的程序function[res,k,pr,wc]=Jacobi(err,M)%输入:err为误差控制,M为最大迭代次数%输出:方程组的解res,迭代次数k和系数矩阵的谱半径pr,精确解与所得解的误差wcformatlongw=input('pleaseinputanumber:');%输入一个正整数w,为hilbert矩阵的维数A=hilb(w
8、);%根据输入的正整数构造hilbert矩阵Am=size(A);%求A的维数n=m(1);xjq=ones(n,1);%方程组的精确解