《北京科技大学应用计算方法作业与答案》由会员分享，可在线阅读，更多相关《北京科技大学应用计算方法作业与答案（15页珍藏版）》请在金锄头文库上搜索。
1、1 一、第一次作业（一） 2-6 计算下列向量的 1-范数、 -范数、 2- 范数。（1）x=(12,-4,-6,2)T A=12,-4,-6,2 A = 12 -4 -6 2 norm(A,1) ans = 24 norm(A,inf) ans = 12 norm(A,2) ans = 14.1421 (2) x=(1,3,-4)T A=1,3,-4 A = 1 3 -4 norm(A,1) ans = 8 norm(A,inf) ans = 4 norm(A,2) ans = 5.0990 （二） 2-9 计算下列矩阵的行范数、列范数、谱范数、F 范数。（1）112111113A A=3,
2、-1,1;1,1,1;2,1,-1 A = 3 -1 1 1 1 1 2 1 -1 norm(A,1) ans = 6 norm(A,inf) 2 ans = 5 norm(A,2) ans = 3.7888 norm(A,fro) ans = 4.4721 （2）Ra aaA,00 A=0,1;-1,0 A = 0 1 -1 0 norm(A,1) ans = 1 norm(A,inf) ans = 1 norm(A,2) ans = 1 norm(A,fro) ans = 1.4142 二、第二次作业用牛顿迭代法求方程0133xx在20x附近的根。要求：给成程序和运行结果. 1、牛顿法的基
3、本原理在求解非线性方程0)(xf时，它的困难在于)(xf是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为kx，在kx处的Taylor展开式为)()()( kkkxxxfxfxf)1.2(令0)(xf，解其方程得到), 1 ,0( ,)()(1kxfxfxxkk kk)2 .2(3 )2 .2(式为牛顿迭代公式，用牛顿迭代公式求方程0)(xf根的方法称为牛顿迭代法。此即牛顿迭代法的设计原理。 2、matlab程序代码function root=NewtonRoot(f,a,b,eps) if (nargin=3) eps=1.0e-4; endf1=subs(sym(f),findsy
4、m(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)
5、,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); endend 3、运行结果截图4 结论：通过计算可以看出0133xx在20x附近的根为1.8794. 三、第三次作业编写高斯顺序消元法求解下面方程组的程序并计算结果。2 .3525.82102 .610321321321xxx
6、xxxxxx)1.3(1、高斯顺序消元法的设计原理高斯顺序消元法的基本思想是将线性方程组nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111)2 .3(通过消元，逐步转化为等价的上（或下）三角形方程组，然后用回代法求解。 2、matlab程序代码function x,XA=GaussXQByOrder(A,b) % 高斯顺序消元法N = size(A); n = N(1); for i=1:(n-1) for j=(i+1):n if (A(i,i)=0) disp( 对角元素为 0！ ); % 防止对角元素为0return ; endl =
7、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; endend x=SolveUpTriangle(A,b); % 通用的求上三角系数矩阵线性方程组的函数XA = A; % 消元后的系数矩阵function x=SolveUpTriangle(A,b) 5 N=size(A); n=N(1); for i=n:-1:1 if (i=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endend3
8、、运行结果截图7 结论：)1.4(式方程组的解为2200.1,1800.1,8600.0321xxx，需要迭代 12 步。（二）用Seidel求解题设方程组 1、Seidel迭代原理在Jacobi迭代计算过程中可看出，Jacobi方法计算)1(k ix时，并未用到已算出的)1( 1)1( 2)1( 1,k ikkxxx，这时想到，如果Jacobi迭代收敛，)1, 1()1(ijxk j比)1, 1()(ijxk j更接近方程组的解，若能在迭代过程中尽快用新的信息)1(k jx去替换)(k jx，则可望收敛更快。由此，可将Jacobi迭代公式改写为,2, 1 ,0,2, 1),(11)(11)1
9、()1(knixaxabaxnijk jijijk jiji iik i)5.4(此式即为Seidel迭代法。 2、Seidel迭代法matlab程序代码function x,n=gauseidel(A,b,x0,eps,M) if nargin=3 eps=1.0e-3; M=200; elseif nargin=4 M=200; elseif nargin=eps x0=x; x=G*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endend3、运行结果截图结论：)1.4(式方程组的解为2200.1,1800.1,86
10、00.0321xxx，需要迭代 5 步。五、第五次作业用归一化算法（归一化幂法6-1）求矩阵 A的最大模特征值和特征向量，其中TVA)1 ,1 , 1(,322244461)0(1、幂法的基本思想任取一个非零初始向量nRv0且00v，由矩阵A的乘幂构造一迭代序列为9 0102 1201vAAvvvAAvvAvvk kk)1.5(假设矩阵A有n个线性无关的特征向量),2, 1(nkxk，于是给定的初始向量0v可以用这组特征向量线性表示，即niiinnxaxaxaxav122110并设0ia。把0v代入迭代序列的第一条式子，得niiiinnnnnxaxaxaxaxaxaxaAAvv12221112
11、21101)(同理可得),2, 1( , 1222111nkxaxaxaxaAvvniik iink nnkk kk2、matlab程序代码function l,v,s=pmethod(A,x0,eps) if nargin=2 eps = 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) x=1:12; y=80.9 67.2 67.1 50.5 32.0 33.6 36.6 46.8 52.3 6
12、2.0 64.1 71.2; xx=0.0:0.1:12.0; yy=spline(x,y,xx); plot(x,y,o,xx,yy) 图 1 通过图像可以看出日照随时间的大概走势，我们也可以推测出以后每个月的大概日照长度的大概走势。12 八、第八次作业在某化学反应中，由实验得分解物浓度与时间关系如下：表 2：分解物浓度与时间关系时间0 5 10 15 20 25 浓度0 1.27 2.16 2.86 3.44 3.87 时间30 35 40 45 50 55 浓度4.15 4.37 4.51 4.58 4.62 4.64 用最小二乘法求)(tfy。解：利用matlab 数值逼近可以得到：即
13、最小二乘拟合二次多项式为2453.01992.00022.0)(2tttf其中拟合值和实际值的对比图见下: 图 2：图中圆圈为已知数据的点，曲线为最小二乘拟合二次多项式九、第九次作业用龙贝格求积算法计算积分dxeIx10。计算过程中数值保留6 位有效数字。1、龙贝格算法原理13 龙贝格求积法的算法为其中)(hTm逼近I的阶达到m2次。2、matlab程序代码function R,quad,err,h=romber(f,a,b,n,tol) % f是被积函数。% a, b分别为是积分的上下限。% n1是T数表的列数。% tol 是允许误差。% R是T数表。% quad是所求积分值。M=1; h=
14、b-a; err=1; J=0; R=zeros(4,4); R(1,1)=h*(feval( f,a)+feval( f,b)/2; while(errtol) h=h/2; s=0; for p=1:M x=a+h*(2*p-1); s=s+feval(f,x); endR(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K)/(4K-1); enderr=abs(R(J,J)-R(J+1,K+1); %取绝对值。endquad=R(J+1,J+1); 再定义函数function y=f(x);
15、 y=exp(x); %exp(x)以e为底的指数。命令窗口中输入,2, 1,41)(4)2( )()()(11mhThT hThThTmmm mm14 clc clear romber(f, 0, 1, 7, 10(-6) 3、运行结果截图结论：求解得最终结果为1.7183 。十、第十次作业用经典 R-K方法计算初值问题 2)0(38yyy， 取步长2 .0h，计算)4 .0(y的近似值1、经典 R-K 方法设计原理在,1nnxx区间上取),(yxf在P个点处的数值作线性组合，即得到 P级的 R-K 方法， 在 P级 R-K 方法中，最常用的是标准的（或经典的）四阶R-K 方法。)22(6),() 2, 21()2,21(),(432113423121kkkkhyyhkyhxfkkhyhxfkkhyhxfkyxfknnnnnnnnnn)1 .10(2、matlab程序代码function fv = Funval(f,varvec,varval) var = findsym(f); varc = findsym(varvec); s1 = length(var); s2 = length(varc); m =floor(s1-1)/3+1); varv = zeros(1,m); if s1 = s2 for i=0: (
