北京科技大学应用计算方法作业与答案 .doc
一、第一次作业（一）2-6计算下列向量的1-范数、-范数、2-范数。（1）x=(12,-4,-6,2)T>> A=12,-4,-6,2A = 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,-4A = 1 3 -4>> norm(A,1)ans = 8>> norm(A,inf)ans = 4>> norm(A,2)ans = 5.0990（二）2-9 计算下列矩阵的行范数、列范数、谱范数、范数。（1）>> A=3,-1,1;1,1,1;2,1,-1A = 3 -1 1 1 1 1 2 1 -1>> norm(A,1)ans = 6>> norm(A,inf)ans = 5>> norm(A,2)ans = 3.7888>> norm(A,'fro')ans = 4.4721（2）>> A=0,1;-1,0A = 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二、第二次作业用牛顿迭代法求方程在附近的根。要求：给成程序和运行结果.1、 牛顿法的基本原理在求解非线性方程时，它的困难在于是非线性函数，为克服这一困难，考虑它的线性展开。设当前点为，在处的展开式为 令，解其方程得到 式为牛顿迭代公式，用牛顿迭代公式求方程根的方法称为牛顿迭代法。此即牛顿迭代法的设计原理。2、程序代码function root=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*f2>0) disp('两端点函数值乘积大于0！') return;else tol=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(dfa>dfb) root=a-fa/dfa; else root=b-fb/dfb; end while(tol>eps) 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、运行结果截图结论：通过计算可以看出在附近的根为1.8794.三、第三次作业编写高斯顺序消元法求解下面方程组的程序并计算结果。 1、 高斯顺序消元法的设计原理高斯顺序消元法的基本思想是将线性方程组 通过消元，逐步转化为等价的上（或下）三角形方程组，然后用回代法求解。2、 程序代码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！'); %防止对角元素为0 return; end l = 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)N=size(A);n=N(1);for i=n:-1:1 if(i<n) s=A(i,(i+1):n)*x(i+1):n,1); else s=0; end x(i,1)=(b(i)-s)/A(i,i);end 3、 运行结果截图结论：高斯顺序消元法求解出的结果为。四、第四次作业编写迭代法和迭代法求解方程组的程序，并计算出结果。 精度要求：（一） 用求解题设方程组1、迭代原理设有一个元线性方程组 它的矩阵形式为，如果非奇异，且。由上式方程组可以得到 而其相应的迭代公式 此迭代公式即为迭代。2、迭代法程序代码functionx,n=jacobi(A,b,x0,eps,varargin)% 求解线性方程组的迭代法其中% A为方程组的系数矩阵% b为方程组的右端项% eps为精度要求，默认值为1e-5% varargin为最大迭代次数，值100% x为方程组的解% n为迭代次数if nargin=3 eps=1.0e-6; M=200;elseif nargin<3 error returnelseif nargin=5 M=varargin1;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=D(L+U);f=Db;x=B*x0+f;n=1;while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; endend3、运行结果截图结论：式方程组的解为，需要迭代12步。（二） 用求解题设方程组1、迭代原理在迭代计算过程中可看出，方法计算时，并未用到已算出的，这时想到，如果迭代收敛，比更接近方程组的解，若能在迭代过程中尽快用新的信息去替换，则可望收敛更快。由此，可将迭代公式改写为 此式即为迭代法。2、迭代法程序代码functionx,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps=1.0e-3; M=200;elseif nargin=4 M=200;elseif nargin<3 error return;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1;while norm(x-x0)>=eps x0=x; x=G*x0+f; n=n+1; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end3、运行结果截图结论：式方程组的解为，需要迭代5步。五、第五次作业用归一化算法（归一化幂法6-1）求矩阵A的最大模特征值和特征向量，其中1、 幂法的基本思想任取一个非零初始向量且，由矩阵的乘幂构造一迭代序列为 假设矩阵有个线性无关的特征向量，于是给定的初始向量可以用这组特征向量线性表示，即并设。把代入迭代序列的第一条式子，得同理可得2、程序代码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)<eps) l = m; %到所需精度，退出，l为主特征值 s = k; %s为迭代步数 return; else if(k=M) disp('迭代步数太多，收敛速度太慢！'); l = m; s = M; else l = m; end endend3、运行结果截图结论：由输出结果可知，经过15步迭代，求得矩阵特征值为9.3597，对应的迭代向量为（1.0000,0.9785，0.6222）T六、第六次作业已知的函数值和导数值如下：求次数小于等于4的多项式，使并给出余项公式。解：首先构造拉格朗日插值函数化简得构造函数对上式求导得根据导数信息可以解出所以余项构造辅助函数其中然后反复利用定理最好可以得到七、第七次作业已知某地区在不同月份的平均日照时间的观测数据如下表，试分析日照时间的变化规律。表1：日照时间月份123456日照h/月80.967.267.150.532.033.6月份789101112日照h/月36.646.852.362.064.171.2解：首先在matlab中输入语句>> x=1:12;>> y=80.9 67.2 67.1 50.5 32.0 33.6 36.6 46.8 52.3 62.0 64.1 71.2;>> xx=0.0:0.1:12.0;>> yy=spline(x,y,xx);>> plot(x,y,'o',xx,yy)图1通过图像可以看出日照随时间的大概走势，我们也可以推测出以后每个月的大概日照长度的大概走势。八、第八次作业在某化学反应中，由实验得分解物浓度与时间关系如下：表2：分解物浓度与时间关系时间0510152025浓度01.272.162.863.443.87时间303540455055浓度4.154.374.514.584.624.64用最小二乘法求。解：利用matlab数值逼近可以得到：即最小二乘拟合二次多项式为其中拟合值和实际值的对比图见下:图2：图中圆圈为已知数据的点，曲线为最小二乘拟合二次多项式九、第九次作业用龙贝格求积算法计算积分。计算过程中数值保留6位有效数字。1、 龙贝格算法原理龙贝格求积法的算法为其中逼近的阶达到次。2、 程序代码functionR,quad,err,h=romber(f,a,b,n,tol)% f是被积函数。% a, b分别为是积分的上下限。% n1是T数表的列数。% tol是允许误差。% R是T数表。% quad是所求积分值。M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(feval('f',a)+feval('f',b)/2;while(err>tol)&(J<n)|(J<4) J=J+1; h=h/2; s=0; for p=1:M x=a+h*(2*p-1); s=s+feval('f',x); end R(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); end err=abs(R(J,J)-R(J+1,K+1);%取绝对值。endquad=R(J+1,J+1);再定义函数function y=f(x);y=exp(x);%exp(x)以e为底的指数。命令窗口中输入clcclearromber('f', 0, 1, 7, 10(-6) 3、 运行结果截图结论：求解得最终结果为1.7183。十、第十次作业用经典R-K方法计算初值问题，取步长，计算的近似值 1、 经典R-K方法设计原理在区间上取在个点处的数值作线性组合，即得到P级的R-K方法，在P级R-K方法中，最常用的是标准的（或经典的）四阶R-K方法。 2、 程序代码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 = s2for i=0: (s1-1)/3)k = findstr(varc,var(3*i+1);index = (k-1)/3;varv(i+1) = varval(index+1);endfv = subs(f,var,varv);elsefv = subs(f,varvec,varval);endfunction y=DELGK4_lungkuta(f, h,a,b,y0,varvec)%f:一阶常微分方程的一般表达式的右端函数%h:积分步长%a:自变量取值下限%b:自变量取值上限%y0:函数初值%varvee：常微分方程的变量组format long;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);for i=2:N+1 K1=Funval(f,varvec,x(i-1) y(i-1); K2=Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3=Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2); K4=Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformat short;3、运行结果截图结论：最终求得的近似值为2.4654.
