北科大——计算方法上机作业——卫鸿儒

计算方法作业姓名：学号：专业：学院：成绩：__________________任课教师：卫宏儒2012年11月作业一：1、计算下列向量的1-范数、-范数、2-范数。(1)(12,4,6,2)Tx(2)(1,3,4)Tx解：(1)(12,4,6,2)Tx程序：x=[12,-4,-6,2];norm(x,1)norm(x,inf)norm(x,2)运行结果：得到(12,4,6,2)Tx的1-范数、-范数、2-范数分别为：24、12、14.1421。(2)(1,3,4)Tx程序：x=[1,3,4];norm(x,1)norm(x,inf)norm(x,2)运行结果：得到(1,3,4)Tx的1-范数、-范数、2-范数分别为：8、4、5.0990。2、计算下列矩阵的行范数、列范数、谱范数、F范数。(1)311111211A(2)0-0aAa,aRaR解：(1)311111211A程序：clcx=[3-11;111;21-1];norm(x,inf)norm(x,1)norm(x,2)norm(x,'fro')运行结果：(2)0-0aAa,aRaR取a=1;程序：clcx=[01;-10];norm(x,inf)norm(x,1)norm(x,2)norm(x,'fro')得到：矩阵311111211A的行范数、列范数、谱范数、F范数分别为：5、6、3.7888、4.4721；0-0aAa，当a=1时，其值分别为：1、1、1、1.4142。作业二：1、用牛顿迭代法求方程0133xx在20x附近的根。要求:给出程序和运行结果，计算结果保留4位有效数字。解：设定初值p0，最大迭代步数N，误差限Tol，近似根p1，迭代步数k，则程序：symsx;f=x^3-3*x-1;%定义f以及df函数df=3*x^2-3;p0=2;N=1000;Tol=1e-5;fork=1:Np1=p0-subs(f/df,x,p0);%subs为置换函数，把f/df表达式中的x用p0替代ifabs(p1-p0)Tolbreakendp0=p1;enddisp(p1);disp(k);%迭代完毕后输出p1和k的值运行结果：有运行结果可知：当迭代4次时，方程的解为1.879，已经满足精度要求。作业三：1、编写高斯顺序消元法求解下面方程组的程序及并计算结果。123123123106.21028.5253.2xxxxxxxxx解：原理：高斯消去法的求解过程：分为两个阶段:首先，把原方程组化为上三角形方程组，称之为“消去”过程；然后，用逆次序逐一求出三角方程组(原方程组的等价方程组)的解，并称之为“回代”过程。程序：%A为系数矩阵，B为增广矩阵，n为方程中未知数的个数A=[10-1-1;-110-2;-2-15];b=[6.28.53.2]';n=3;B=[A,b];fork=1:n-1ifB(k,k)==0disp('Gaussfailstosolvethis');returnendfori=k+1:nm=B(i,k)/B(k,k);B(i,k+1:n+1)=B(i,k+1:n+1)-m*B(k,k+1:n+1);%顺序消元过程endendX(n)=B(n,n+1)/B(n,n);%回代过程fori=n-1:-1:1X(i)=(B(i,n+1)-B(i,i+1:n)*X(i+1:n)')/B(i,i);enddisp(X)运行结果：结果显示：经高斯顺序消元法得到线性方程组的解为：10.86x，21.18x，31.22x。作业四：1、编写Jacobi迭代法和Seidel迭代法求解上述方程组的程序，并计算出结果。精度要求：(1)()310kkxx解：(1)Jacobi迭代程序：A=[10-1-1;-110-2;-2-15];b=[6.28.53.2]';n=3;Tol=1e-3;N=1000;IX=[000]';X=IX;fork=1:Nfori=1:nX(i)=(b(i)-A(i,:)*IX+A(i,i)*IX(i))/A(i,i);%把X(i)项用其余各X示ifnorm(X-IX,inf)Toldisp(X);disp(k);return;endIX=X;endend运行结果：由运行结果显示，经5次Jacobi迭代之后，就能得到满足精度要求的线性方程组解。(2)高斯赛德尔迭代程序：A=[10-1-1;-110-2;-2-15];b=[6.28.53.2]';Tol=1e-3;N=1000;X0=[1;1;1];X=X0;k=0;BGS=-(tril(A))^(-1)*(triu(A)-diag(diag(A)));s=max(abs(eig(BGS)))ifs1,warning('发散');end;n=length(b);whilekNfori=1:nifi==1X(1)=(b(1)-A(1,2:n)*X0(2:n))/A(1,1);elseifi==nX(n)=(b(n)-A(n,1:n-1)*X(1:n-1))/A(n,n);elseX(i)=(b(i)-A(i,1:i-1)*X(1:i-1)-A(i,i+1:n)*X0(i+1:n))/A(i,i);endendendifabs(X-X0)Tolbreak;endX0=X;k=k+1;endifk==N,warning('已达到迭代次数上限');end;disp(['k=',num2str(k)]),X运行结果：由运行结果可知：高斯赛德尔迭代法收敛速度明显比Jacobi迭代法快，其只要迭代3次，就能满足精度要求。作业五：1、用归一化算法(归一化幂法)求矩阵A的最大摸特征值和特征向量，其中164442223A,(0)(111)TV，，计算过程中保留小数点后5位。解：程序：A=[164;442;223];%输入所需要的各个条件v=[111]';eps=1e-5;N=1000;lamda=0;erro=1;k=1;while(k=N&&erroeps)u=A*v;m=max(abs(u));%求出向量u中模最大的元素dc=abs(lamda-m);u=u/m;%进行归一化处理dv=norm(u-v);err=max(dc,dv);v=u;lamda=m;k=k+1;enddisp('Themaximumeigenvalueis:');%输出结果disp(lamda);disp('itscorrespondingeigenvectoris:')disp(u);运行结果：运行结果显示：矩阵A的最大摸特征值和特征向量分别为9.3579和(1.0000,0.9785,0.6222)TV。作业六：1、已知fx的函数值和导数值如下:(1)1f，(2)2f，(3)3f，'(1)1f，'(2)2f。求次数小于等于4的多项式()Px，使得：(1)1P，(2)2P，(3)3P，'(1)1P，'(2)2P。并给出余项公式。解：0201122012010210122021()()()()()()()()()()()()()xxxxxxxxxxxxLxyyyxxxxxxxxxxxx2(2)(3)(1)(3)(1)(2)()123(12)(13)(21)(23)(31)(32)xxxxxxLxx构造多项式：2()()()pxLxqx，1,1AB则有：(1)(2)(3)0qqq所以：()(1)(2)(3)()qxxxxrx因为()px为四次多项式，知()rx为一次多项式，所以设()rxABx,则()(1)(2)(3)()pxxxxxABx由导数值值可知：'(1)1()21pAB'(2)1(2)2pAB解得：1,1AB所以432()(1)(1)(2)(3)7176Pxxxxxxxxxx522()()(1)(2)(3)5!Rxxxxf作业七：1、一维插值问题的应用题：已知某地区在不同月份的平均日照时间的观测数据如下表所示(h/月)，试分析日照时间的变化规律。月份123456789101112日照80.967.267.150.532.033.636.646.852.362.064.171.2解：程序：x=1:1:12;y=[80.967.267.150.532.033.636.646.852.362.064.171.2];xi=1:0.01:12;y1i=interp1(x,y,xi);y2i=interp1(x,y,xi,'spline');plot(x,y,'o',xi,y1i,'-',xi,y2i,'-.')运行结果：运行结果中，圆圈是已知数据点，折线是线性插值曲线，点划线是三次样条曲线。作业八：1、在某化学反应中，由实验得分解物浓度和时间关系如下：时间t/s0510152025303540455055浓度y(*10-4)mol/m301.272.162.863.443.874.154.374.514.584.624.64用最小二乘法求()yft。解：程序：clct=[0:5:55];y=[01.272.162.863.443.874.154.374.514.584.624.64];p=polyfit(t,y,2)ti=0:5:55;yi=polyval(p,ti);plot(t,y,'o',ti,yi)运行结果：作业九：1、用龙贝格求积算法计算积分10xIedx。计算过程中数值保留6位有效数字。解：程序：clca=0;b=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(f(a)+f(b))/2;while(erreps)J=J+1;h=h/2;x=a+h:2*h:b-h;R(J+1,1)=R(J,1)/2+h*sum(f(x));forK=1:min(3,J)R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K))/(4^K-1);endif(J3)err=abs(R(J+1,4)-R(J,4));endendquad=R(J,4);vpa(quad,6)其中，10()xfxedx。运行结果：作业十：1、用经典R-K方法计算初值问题'83yy，(0)2y，取步长0.2h，计算(0.4)y的近似值，小数点后至少保留4位。解：程序：clcx=0;y=2;h=0.2;n=2;fori=1:nK1=h*ff(x,y);K2=h*ff(x+h/2,y+K1/2);K3=h*ff(x+h/2,y+K2/2);K4=h*ff(x+h,y+K3);y=y+(K1+2*K2+2*K3+K4)/6;x=x+h;enddisp(y)运行结果：