欢迎来到天天文库
浏览记录
ID:21850737
大小:547.82 KB
页数:17页
时间:2018-10-25
《矩阵与数值分析》由会员上传分享,免费在线阅读,更多相关内容在学术论文-天天文库。
1、2013级工科硕士研究生《矩阵与数值分析》课程数值实验题目N1Q6一、设分别编制从小到大和从大到小的顺序程序分别计算7=2J-1o10000,J1000000并指山两种方法计算结果的有效位数。Matlab程序如下:function[$i,sd]=S(N)formatlong;si=0;sd=0;forj=N:-l:2si=1.0e6/(jA2-l)+si;endforj=2:Nsd=1.0e6/(jA2-l)+sd;endend在matlab命令窗口中输入:[si,sd]=S(10000)运行结果:si=7.499000049995000e+005sd=7.49900004
2、9994994e+005在matlab命令窗口中输入:[si,sd]=S(1000000)运行结果:si=7.499990000005000e+005sd=7.499990000005200e+0051结果分析:si为从大到小的顺序求和的值,sd为从小到大的顺序求和的值。当N分别为10000和1000000时,si分别为7.499000049995000e+005和7.499990000005000e+005,可以看出这两个数的有效值均为13位;而sd分别为7.499000049994994e+005和7.499990000005200e+005,这两个数的有效值均为16位
3、。这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。为了使结果更为精确我们必须避免在四则运算屮出现“大数吃小数”的情况,应该按从小到大的顺序进行求和。二、解线性方程组1.分別利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组/U=/>,其屮常向fi为维随机生成的列向y:,系数矩阵A具有如下形式其中1-、z2-1Tn^+2/^,w-lH-l为《-1阶矩阵,L为n-1阶单位矩阵,迭代法计算停止的条件-12为:k+r'l4、=jz(n)fori=l:n-l%此处n赋值即可,如n=100forj=l:n-lif(i==j)T(ij)=2;endif(i==(j+D)T(i,j)=-1;endif(i==(j_l))T(ij)=-1;endendendI=eye(n-1);k=l;formm=l:(n-l)A(k:(k+n-2),k:(k+n-2))=T+2*I;k=k+n-l;endk=l;forxx=l:(n-2)A(k:(k+n-2),(k+n-l):(k+2*n-3))=-I;A((k+n-l):(k+2*n-3),k:(k+n-2))=-I;k=k+n-l;endb=randn((n-15、)A2,1);xO=zeros((n-l)A2,1);Jacobi迭代法Matlab程序如下:functionn=jacobi(A,b,xO)eps=1.0e-10;M=10000;D=diag(diag(A));L=-tril(A,-l);U=-triu(A,l);B=D(L+U);f=Db;x=B*xO+f;n=l;whilenorm(x-xO)>=epsx0=x;x=B*xO+f;n=n+l;if(n〉=M)dispCWarning:迭代次数太多,可能不收敛return;endendGauss-Seidel迭代法Matlab程序如下:functionn=gause6、idel(A,b,xO)eps=1.0e-10;%求A的对角矩阵%求A的下三角阵%求A的上三角阵%迭代次数M=10000;D=diag(diag(A));L=-tril(A,-l);U=-triu(A,l);G=(D-L)U;f=(D-L)b;x=G*xO+f;n=l;whilenorm(x-xO)〉=epsx0=x;x=G*xO+f;n=n+l;if(n>=M)dispCWarning:return;end%求A的对角矩阵%求A的下三角阵%求A的上三角阵%迭代次数迭代次数太多,可能不收敛end在matlab命令窗口中输入:[A,b,x0J=jz(10);J10=jac7、obi(A,b,x0)G1O=gauseidel(A,b,xO)[A,b,x0]=jz(20);J20=jacobi(A,b,x0)G20=gauseidel(A,b,x0)[A,b,xO]=jz(3O);J30=jacobi(A,b,x0)G3O=gauseidel(A,b,xO)运行结果:J10=436;GIO=214J20=1810;G20=913J30=3990;G30=2001结果分析:N=10且M=1000吋,Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为436和214;N=20且M
4、=jz(n)fori=l:n-l%此处n赋值即可,如n=100forj=l:n-lif(i==j)T(ij)=2;endif(i==(j+D)T(i,j)=-1;endif(i==(j_l))T(ij)=-1;endendendI=eye(n-1);k=l;formm=l:(n-l)A(k:(k+n-2),k:(k+n-2))=T+2*I;k=k+n-l;endk=l;forxx=l:(n-2)A(k:(k+n-2),(k+n-l):(k+2*n-3))=-I;A((k+n-l):(k+2*n-3),k:(k+n-2))=-I;k=k+n-l;endb=randn((n-1
5、)A2,1);xO=zeros((n-l)A2,1);Jacobi迭代法Matlab程序如下:functionn=jacobi(A,b,xO)eps=1.0e-10;M=10000;D=diag(diag(A));L=-tril(A,-l);U=-triu(A,l);B=D(L+U);f=Db;x=B*xO+f;n=l;whilenorm(x-xO)>=epsx0=x;x=B*xO+f;n=n+l;if(n〉=M)dispCWarning:迭代次数太多,可能不收敛return;endendGauss-Seidel迭代法Matlab程序如下:functionn=gause
6、idel(A,b,xO)eps=1.0e-10;%求A的对角矩阵%求A的下三角阵%求A的上三角阵%迭代次数M=10000;D=diag(diag(A));L=-tril(A,-l);U=-triu(A,l);G=(D-L)U;f=(D-L)b;x=G*xO+f;n=l;whilenorm(x-xO)〉=epsx0=x;x=G*xO+f;n=n+l;if(n>=M)dispCWarning:return;end%求A的对角矩阵%求A的下三角阵%求A的上三角阵%迭代次数迭代次数太多,可能不收敛end在matlab命令窗口中输入:[A,b,x0J=jz(10);J10=jac
7、obi(A,b,x0)G1O=gauseidel(A,b,xO)[A,b,x0]=jz(20);J20=jacobi(A,b,x0)G20=gauseidel(A,b,x0)[A,b,xO]=jz(3O);J30=jacobi(A,b,x0)G3O=gauseidel(A,b,xO)运行结果:J10=436;GIO=214J20=1810;G20=913J30=3990;G30=2001结果分析:N=10且M=1000吋,Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为436和214;N=20且M
此文档下载收益归作者所有