北京科技大学应用计算方法作业与答案

1一、第一次作业（一）2-6计算下列向量的1-范数、∞-范数、2-范数。（1）x=(12,-4,-6,2)TA=[12,-4,-6,2]A=12-4-62norm(A,1)ans=24norm(A,inf)ans=12norm(A,2)ans=14.1421(2)x=(1,3,-4)TA=[1,3,-4]A=13-4norm(A,1)ans=8norm(A,inf)ans=4norm(A,2)ans=5.0990（二）2-9计算下列矩阵的行范数、列范数、谱范数、F范数。（1）112111113AA=[3,-1,1;1,1,1;2,1,-1]A=3-1111121-1norm(A,1)ans=6norm(A,inf)2ans=5norm(A,2)ans=3.7888norm(A,'fro')ans=4.4721（2）RaaaA,00A=[0,1;-1,0]A=01-10norm(A,1)ans=1norm(A,inf)ans=1norm(A,2)ans=1norm(A,'fro')ans=1.4142二、第二次作业用牛顿迭代法求方程0133xx在20x附近的根。要求：给成程序和运行结果.1、牛顿法的基本原理在求解非线性方程0)(xf时，它的困难在于)(xf是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为kx，在kx处的Taylor展开式为))(()()('kkkxxxfxfxf)1.2(令0)(xf，解其方程得到),1,0(,)()('1kxfxfxxkkkk)2.2(3)2.2(式为牛顿迭代公式，用牛顿迭代公式求方程0)(xf根的方法称为牛顿迭代法。此即牛顿迭代法的设计原理。2、matlab程序代码functionroot=NewtonRoot(f,a,b,eps)if(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)roor=b;endif(f1*f20)disp('两端点函数值乘积大于0！')return;elsetol=1;fun=diff(sym(f));fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfadfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(toleps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);root=r1-fx/dfx;tol=abs(root-r1);endend3、运行结果截图4结论：通过计算可以看出0133xx在20x附近的根为1.8794.三、第三次作业编写高斯顺序消元法求解下面方程组的程序并计算结果。2.3525.82102.610321321321xxxxxxxxx)1.3(1、高斯顺序消元法的设计原理高斯顺序消元法的基本思想是将线性方程组nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111)2.3(通过消元，逐步转化为等价的上（或下）三角形方程组，然后用回代法求解。2、matlab程序代码function[x,XA]=GaussXQByOrder(A,b)%高斯顺序消元法N=size(A);n=N(1);fori=1:(n-1)forj=(i+1):nif(A(i,i)==0)disp('对角元素为0！');%防止对角元素为0return;endl=A(j,i);m=A(i,i);A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;%消元方程b(j)=b(j)-l*b(i)/m;endendx=SolveUpTriangle(A,b);%通用的求上三角系数矩阵线性方程组的函数XA=A;%消元后的系数矩阵functionx=SolveUpTriangle(A,b)5N=size(A);n=N(1);fori=n:-1:1if(in)s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end3、运行结果截图结论：高斯顺序消元法求解出)1.3(的结果为2200.1,1800.1,8600.0321xxx。四、第四次作业编写Jacobi迭代法和Seidel迭代法求解方程组的程序，并计算出结果。2.3525.82102.610321321321xxxxxxxxx)1.4(精度要求：3)()1(10kkxx（一）用Jacobi求解题设方程组1、Jacobi迭代原理设有一个n元线性方程组nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111)2.4(它的矩阵形式为BAX，如果nnijaA)(非奇异，且niaii,2,1,0。由上式方程组可以得到6),,2,1(),(111nixabaxnjjjijiiii)3.4(而其相应的迭代公式),,2,1(),(111)(nixabaxnjjkjijiiii)4.4(此迭代公式即为Jacobi迭代。2、Jacobi迭代法matlab程序代码function[x,n]=jacobi(A,b,x0,eps,varargin)%求解线性方程组的迭代法其中%A为方程组的系数矩阵%b为方程组的右端项%eps为精度要求，默认值为1e-5%varargin为最大迭代次数，值100%x为方程组的解%n为迭代次数ifnargin==3eps=1.0e-6;M=200;elseifnargin3errorreturnelseifnargin==5M=varargin{1};endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend3、运行结果截图7结论：)1.4(式方程组的解为2200.1,1800.1,8600.0321xxx，需要迭代12步。（二）用Seidel求解题设方程组1、Seidel迭代原理在Jacobi迭代计算过程中可看出，Jacobi方法计算)1(kix时，并未用到已算出的)1(1)1(2)1(1,,kikkxxx，这时想到，如果Jacobi迭代收敛，)1,,1()1(ijxkj比)1,,1()(ijxkj更接近方程组的解，若能在迭代过程中尽快用新的信息)1(kjx去替换)(kjx，则可望收敛更快。由此，可将Jacobi迭代公式改写为,2,1,0,,,2,1),(11)(11)1()1(knixaxabaxnijkjijijkjijiiiki)5.4(此式即为Seidel迭代法。2、Seidel迭代法matlab程序代码function[x,n]=gauseidel(A,b,x0,eps,M)ifnargin==3eps=1.0e-3;M=200;elseifnargin==4M=200;elseifnargin3errorreturn;endD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);8G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;whilenorm(x-x0)=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend3、运行结果截图结论：)1.4(式方程组的解为2200.1,1800.1,8600.0321xxx，需要迭代5步。五、第五次作业用归一化算法（归一化幂法6-1）求矩阵A的最大模特征值和特征向量，其中TVA)1,1,1(,322244461)0(1、幂法的基本思想任取一个非零初始向量nRv0且00v，由矩阵A的乘幂构造一迭代序列为901021201vAAvvvAAvvAvvkkk)1.5(假设矩阵A有n个线性无关的特征向量),,2,1(nkxk，于是给定的初始向量0v可以用这组特征向量线性表示，即niiinnxaxaxaxav122110并设0ia。把0v代入迭代序列的第一条式子，得niiiinnnnnxaxaxaxaxaxaxaAAvv1222111221101)(同理可得),,2,1(,1222111nkxaxaxaxaAvvniikiinknnkkkk2、matlab程序代码function[l,v,s]=pmethod(A,x0,eps)ifnargin==2eps=1.0e-6;endv=x0;%v为主特征向量M=5000;%迭代步数限制m=0;l=0;for(k=1:M)y=A*v;m=max(y);%m为按模最大的分量v=y/m;if(abs(m-l)eps)l=m;%到所需精度，退出，l为主特征值s=k;%s为迭代步数return;elseif(k==M)disp('迭代步数太多，收敛速度太慢！');l=m;s=M;elsel=m;end10endend3、运行结果截图结论：由输出结果可知，经过15步迭代，求得矩阵特征值为9.3597，对应的迭代向量为（1.0000,0.9785，0.6222）T六、第六次作业已知)(xf的函数值和导数值如下：2)2(,1)1(,3)3(,2)2(,1)1(''fffff求次数小于等于4的多项式)(xP，使2)2(,1)1(,3)3(,2)2(,1)1(''fffff并给出余项公式。解：首先构造拉格朗日插值函数3)23()13()2)(1(2)32()12()3)(1(1)31()21()3)(2(2xxxxxxP化简得xP2构造函数)3)(2)(1)((4xxxbaxxP对上式求导得161122121834223'4abaxbxaxbxaxP根据导数信息可以解出1,1ba11所以)3)(2)(1)(1(4xxxxxP余项)3()2()1)(()()()(2244xxxxkxPxfxR构造辅助函数)3()2()1)(()()()(224tttxktPtft其中)3()2()1()()(224xxxxRxk然后反复利用Rolle定理最好可以得到)3()2()1(!5)()(2254xxxfxR七、第七次作业已知某地区在不同月份的平均日照时间的观测数据如下表，试分析日照时间的变化规律。表1：日照时间月份123456日照h/月80.967.267.150.532.033.6月份789101112日照h/月36.646.852.362.064.