数值分析上机试题对应参考答案.doc

一、 问答题1、什么是近似值x* 有效数字？若近似值x*的误差限是某一位的半个单位，该位到x*的第一位非零数字共有n位，就说x*有n位有效数字。它可表示为X=10m(a1+a210-1+an10-(n-1),其中ai(i=1,2,n)是0到9中的一个数字，a10,m为整数，且xx*10m-n+12、数值计算应该避免采用不稳定的算法，防止有效数字的损失. 因此，在进行数值运算算法设计过程中主要注意什么？（1）简化计算过程，减少运算次数；（2）避免两个相近的数相减；（3）避免除数的绝对值远小于被除数的绝对值；（4）防止大数“吃掉”小数的现象；（5）使用数值稳定的算法，设法控制误差的传播。3、写出“n 阶阵A 具有n 个不相等的特征值”的等价条件（至少写3 个）(1)|A|不为零(2)n阶矩阵A的列或行向量组线性无关(3)矩阵A为满秩矩阵(4)n阶矩阵A与n阶可逆矩阵B等价4、迭代法的基本思想是什么？就是用某种极限过程去逐步逼近线性方程组精确解得方法。其基本思想为：先任取一组近似解初值X0，然后按照某种迭代原则，由X0计算新的近似解X1，以此类推，可计算出X2,X3，XK，。,如果X收敛，则取为原方程组的解。5、病态线性方程组的主要判断方法有哪些？（1）系数矩阵的某两行（列）几乎近似相关（2）系数矩阵的行列式的值很小（3）用主元消去法解线性方程组时出现小主元（4）近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题要求。6、Lagrange 插值的前提条件是什么？并写出二次Lagrange 插值的基函数。 前提条件是： 二次Lagrange 插值的基函数：7、什么是数值积分的代数精度？如果某一个求积公式对于次数不超过m的多项式均能准确地成立，但对于m+1次多项式就不准确成立，则称该求积公式具有m次代数精度（或代数精确度）。8、写出解线性方程组的全主元消去法、列主元消去法和标度化列主元消去法各自的优缺点。一般来说列主元法能保证算法的稳定，所谓算法的稳定是指在运算过程中计算误差（对消去法这种直接法来说主要指由于计算机字长有限带来的舍入误差）能得到控制；全主元是较列主元更稳定的算法，但是它的计算量要比列主元算法大得多，列主元算法每做一次消元仅与同列的元素做比较，比较的次数与线性方程组的阶n是同阶的量，而全主元每做一次要与系数矩阵所有元素进行比较，计算量是与n2同阶的量。9、求矩阵特征值的 QR 方法的算法步骤是什么？对于给定的n阶实对称矩阵A与迭代次数M；（1）令A1=A，对于k=1,2，M；（2）迭代计算下一个矩阵：Ak=QkRk（对Ak作QR分解）， （3）Ak=RQ(交换乘法次序)，令k=k+1,Ak+1=QkTAKQK（4）返回到2，直到k=M,输出A的主对角元素。10、什么是高斯型求积公式与高斯点？一般可设求积公式=In()，为具有 n个节点的插值求积公式，且具有2n+1次最高代数精度，则称其求积点x0,x1,xn称为高斯点具有，相应的公式称为高斯型求积公式。11、 什么是插值计算的龙格现象？增加插值节点，提高插值多项式的次数，可以使插值函数在更多的点与所逼近的函数取相同的值，但会使插值函数在两端发生激烈的振荡，这就是插值计算的龙格现象。12、龙格库塔方法的基本思想是什么？基本思想是，在每一步内多预报几个点的斜率值，将其加权平均作为平均斜率，从而构造出更高的计算格式。具体做法是，用函数f(x,y)在若干点上的函数值的线性组合来构造近似公式，构造时让近似公式在（xi,yi）处的泰勒展开式与解y(x)在xi处的泰勒展开式的前面几项重合，从而使近似公式达到所需要的阶数，这样既避免求偏导数，又提高了计算方法的精度。二、计算题1、随机生成4*4的矩阵，做列主元的三角分解，并验证等式PA=LU 。解： A=rand(4) %随机生成的四位矩阵A = 0.8147 0.6324 0.9575 0.9572 0.9058 0.0975 0.9649 0.4854 0.1270 0.2785 0.1576 0.8003 0.9134 0.5469 0.9706 0.1419 L,U,P=lu(A)L = 1.0000 0 0 0 0.9917 1.0000 0 0 0.8920 -0.3250 1.0000 0 0.1390 -0.4552 0.2567 1.0000U = 0.9134 0.5469 0.9706 0.1419 0 -0.4448 0.0024 0.3447 0 0 0.0925 0.9426 0 0 0 0.6955P = 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 P*A %验证PA=LUans = 0.9134 0.5469 0.9706 0.1419 0.9058 0.0975 0.9649 0.4854 0.8147 0.6324 0.9575 0.9572 0.1270 0.2785 0.1576 0.8003 L*Uans = 0.9134 0.5469 0.9706 0.1419 0.9058 0.0975 0.9649 0.4854 0.8147 0.6324 0.9575 0.9572 0.1270 0.2785 0.1576 0.80032、用Lagrange 插值和Newton 插值拟合0, 2*pi上的sin 函数，并画图比较。解：建立拉格朗日函插值函数function yt,L = LagInterp1(x,y,xt)%拉格朗日差值%x，y：差值条件%xt：用拉格朗日插值函数要计算的自变量，可以是多个%yt：用拉格朗日插值函数计算出xt对应的函数值数组%L：拉格朗日插值多项式表达式syms t;n=length(x);ny=length(y);if n=ny error(差值节点x与函数值y不一致);end L = 0.0;for k = 1:n lk =1;for j = 1:n if j=k lk = lk*(t-x(j)/(x(k)-x(j); end end; L = L + y(k)*lk; endsimplify(L); %简化拉格朗日插值多项式表达式L=collect(L); % 将拉格朗日插值多项式展开yt = subs(L,t,xt); %计算插值点处的函数值 建立牛顿函数：function yt,N = NewtInterp(x,y,xt)%已知数据点的牛顿插值%x，y：差值条件%xt：要计算的插值点，可以是多个%yt：用牛顿插值函数算出xt对应的函数值数组%L：牛顿插值多项式表达式syms t;n=length(x);ny=length(y);if n=ny error(差值节点x与函数值y不一致);end a=zeros(1,n);N = y(1);w = 1;for k=1:n-1 yy=zeros(1,n); % 记录差商 for j=k+1:n yy(j) = (y(j)-y(k)/(x(j)-x(k); end a(k) = yy(k+1); w = w*(t-x(k); N = N + a(k)*w; y = yy;endyt = subs(N,t,xt);simplify(N);N = collect(N); % 将插值多项式展开N = vpa(N, 6); % 系数转化为6位精度命令： x=0:pi/10:2*pi; y=sin(x); z=0:pi/20:2*pi; y1=LagInterp1(x,y,z); %拉格朗日拟合法 y2=NewtInterp(x,y,z); %牛顿拟合法 figure; plot(z,sin(z),-r,z,y1,-b,z,y2,-.k) %绘制函数图象 hold on plot(x,y,+) xlabel(x) ylabel(y)3、对不同初值用Jacobi 迭代法解方程组Ax=b，其中解：建立雅可比函数function tx=jacobi(A,b,imax,x0,acc)%利用Jacobi迭代法解线性方程组AX=b，迭代初值为x0，迭代次数为imax，精度为acc%利用Jacobi迭代法解线性方程组AX=b，%迭代初值为x0，%迭代次数为imax，%精度为accdel=10(-10); %主对角元素不能太小 tx=x0;n=length(x0);for i=1:ndg=A(i,i);if abs(dg)del disp(主对角元素太小); returnendendfor k=1:imax for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x0(j); end end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0)A=2,-1,0,0,0,0;-1,2,-1,0,0,0;0,-1,2,-1,0,0;0,0,-1,2,-1,0;0,0,0,-1,2,-1,;0,0,0,0,-1,2;b=1,0,1,0,0,1;acc=1.0*10(-6); imax=20; （10,20,30都可取） x0=zeros(1,6); tx=jacobi(A,b,imax,x0,acc) （此命令一行一行的复制）tx = 0 0 0 0 0 0 0.5000 0 0.5000 0 0 0.5000 0.5000 0.5000 0.5000 0.2500 0.2500 0.5000 0.7500 0.5000 0.8750 0.3750 0.3750 0.6250 0.7500 0.8125 0.9375 0.6250 0.5000 0.6875 0.9063 0.8438 1.2188 0.7188 0.6563 0.7500 0.9219 1.0625 1.2813 0.9375 0.7344 0.8281 1.0313 1.1016 1.5000 1.0078 0.8828 0.8672 1.0508 1.2656 1.5547 1.1914 0.9375 0.9414 1.1328 1.3027 1.7285 1.2461 1.0664 0.9688 1.1514 1.4307 1.7744 1.3975 1.1074 1.0332 1.2153 1.4629 1.9141 1.4409 1.2153 1.0537 1.2314 1.5647 1.9519 1.5647 1.2473 1.1077 1.2823 1.5917 2.0647 1.5996 1.3362 1.1237 1.2958 1.6735 2.0956 1.7004 1.3616 1.1681 1.3368 1.6957 2.1870 1.7286 1.4343 1.1808 1.3479 1.7619 2.2122 1.8106 1.4547 1.2171 1.3809 1.7800 2.2862 1.8335 1.5139 1.2274 1.3900 1.8336 2.3067 1.9001 1.5304 1.2569 1.4168 1.8484 2.3668 1.9186 1.5785 1.2652 1.4242 1.8918 2.3835 1.9727 1.5919 1.28934、用最小二乘法求一个多项式，使它与下列数据相拟合x-1.0-0.50.00.51.01.52.0y-4.447-0.4520.5510.048-0.4470.5494.552(1) 求拟合曲线y=Pn(x), (n=1, 2, 3)(2) 求拟合误差 (3) 画出y=Pn(x)的图像解：（1）命令 x=-1.0 -0.5 0.0 0.5 1.0 1.5 2.0; y=-4.447 -0.452 0.551 0.048 -0.447 0.549 4.552; p1=polyfit(x,y,1) %一次多项式拟合p1 = 2.0001 -0.9495 p2=polyfit(x,y,2) %二次多项式拟合p2 = 0.0010 1.9991 -0.9502 p3=polyfit(x,y,3) %三次多项式拟合p3 = 1.9991 -2.9977 -0.0000 0.5491 z=-1.0:0.1:2.0; y1=polyval(p1,z); y2=polyval(p2,z); y3=polyval(p3,z);即拟合的多项式为：y1 = 28001/14000*x-5317/5600y2 = 1152921504607243/1152921504606846976*x2+27987/14000*x-13303/14000y3 = 2249/1125*x3-8993/3000*x2-45750853357891/1152921504606846976*x+23063/42000（2）（3） plot(x,y,o,z,y1,-r,z,y2,-b,z,y3,-k) %绘制拟合函数图象 5、分别用欧拉公式、梯形公式和改进欧拉方法求解：取步长h=0.1. 要求列表输出解的结果，并比较不同方法的精度。解：建立函数 欧拉公式functionx,y=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n); end x=x;y=y; 梯形公式的Matlab实现程序functionx,y=tixing(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x;y=y; function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1;if kK,error(迭代发散);end end 改进欧拉公式的Matlab实现程序: Function x,y=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(k1+k2)/2; end x=x;y=y; dyfun=inline(y+2*x/y2); x,y=naeuler(dyfun,0,1,1,0.1); %调用欧拉公式 x,yans = 0 1.0000 0.1000 1.1000 0.2000 1.2265 0.3000 1.3758 0.4000 1.5450 0.5000 1.7331 0.6000 1.9397 0.7000 2.1655 0.8000 2.4119 0.9000 2.6806 1.0000 2.9737 dyfun=inline(y+2*x/y2); x,y=tixing(dyfun,0,1,1,0.1); %调用梯形公式 x,yans = 0 1.0000 0.1000 1.1286 0.2000 1.2810 0.3000 1.4549 0.4000 1.6492 0.5000 1.8644 0.6000 2.1017 0.7000 2.3631 0.8000 2.6510 0.9000 2.9682 1.0000 3.3182 dyfun=inline(y+2*x/y2); x,y=naeuler2(dyfun,0,1,1,0.1); %调用改进欧拉公式 x,yans = 0 1.0000 0.1000 1.1133 0.2000 1.2520 0.3000 1.4128 0.4000 1.5936 0.5000 1.7939 0.6000 2.0143 0.7000 2.2560 0.8000 2.5207 0.9000 2.8107 1.0000 3.1287 x=0:0.1:1; y=sqrt(1+2.*x); %精确解 x,yans = 0 1.0000 0.1000 1.0954 0.2000 1.1832 0.3000 1.2649 0.4000 1.3416 0.5000 1.4142 0.6000 1.4832 0.7000 1.5492 0.8000 1.6125 0.9000 1.6733 1.0000 1.73216、分别用复化梯形公式、复化辛普森公式计算（区间10 等分），并与积分精确值作比较。解：建立函数复化梯形公式function I=T_quad(x,y)n=length(x);m=length(y);if n=m error(wrong); return;endh=(x(n)-x(1)/(n-1);a=1 2*ones(1,n-2) 1;I=h/2*sum(a.*y);建立复化辛普森公式function I=S_quad(x,y)n=length(x);m=length(y);if n=m error(hh); return; endif rem(n-1,2)=0 I=T_quad(x,y); return;endN=(n-1)/2;h=(x(n)-x(1)/N;a=zeros(1,n);for k=1:N a(2*k-1)=a(2*k-1)+1;a(2*k-1)=a(2*k)+4; a(2*k+1)=a(2*k+1)+1;endI=h/6*sum(a.*y);命令 x=0:0.1:1; y=sqrt(1-exp(-x)./exp(x); T=T_quad(x,y) %复化梯形公式求积分T = 0.3284 S=S_quad(x,y) %复化辛普公式求积分S = 0.2012 syms x; %定义符号变量 a=int(sqrt(1-exp(-x)/exp(x),0,1) %使用int函数求函数积分a =(2*(1 - 1/exp(1)(3/2)/3 %有理说解 vpa(a) %转化为数值解ans =0.33504922214016864454877581388462综合题第一题命令：v0=100 %初速度L=200 %击球处离本垒的水平距离g=9.8/0.3048 %重力加速度，转为英尺/平方秒 syms aa=30 %速度与水平夹角x=a*pi/180 %转为弧度制 vx=v0*cos(x) %水平速度vy=v0*sin(x) %垂直方向速度t=L/vx LH=vy*t-0.5*g*t2 if(LH+335) disp(可以飞过围墙)else disp(不可以飞过围墙)end syms x vx=v0*cos(x)vy=v0*sin(x) t=L/vxsyms hh=vy*t-0.5*g*t2-32 fzero(inline(h),0)v0 = 100L = 200g = 32.1522a = 30x = 0.5236vx = 86.6025vy = 50.0000t = 2.3094LH = 29.7308不可以飞过围墙 vx = 100*cos(x) vy = 100*sin(x) t = 2/cos(x) h = 200*sin(x)/cos(x)-24500/381/cos(x)2-32 ans = 0.5372验证性的题目：2.2.2 第1题（1） fun=inline(x2-0.8*x+0.15);x_star,inedx,it=bisect(fun,0.1,0.6)x_star = 0.0800 0.0300inedx = 0it = 0当index=0时，表明初始区间不是有根区间。（2）（3）3.2.1 第1题A=1 2 -1;0 1 2;-3 6 4;B=-1 0 1;0 2 2;3 5 1C=BD=A+BE=A*BF=A2G=inv(A)*BH=A*BB = -1 0 1 0 2 2 3 5 1C = -1 0 3 0 2 5 1 2 1D = 0 2 0 0 3 4 0 11 5E = -4 -1 4 6 12 4 15 32 13F = 4 -2 -1 -6 13 10 -15 24 31G = -1.0000 0.1304 1.3478 0 0.3478 0.2609 0 0.8261 0.8696H = -4 -1 4 6 12 4 15 32 13第2题 A=1 1 1 1 1;2 1 1 1 1;3 2 1 1 1;4 3 2 1 1;5 4 3 2 1;b=5 6 8 11 15;x,index=Gauss(A,b)x = 1 1 1 1 1index = 1第3题 A=-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4;b=1 1 1 1;ep=1e-15;it_max=100;x0=0;x,k,index=Gau_Seid(A,b,ep,it_max)x = -1.0000 -1.0000 -1.0000 -1.0000k = 61index = 1第4题 A=2 -1 0 0;-1 2 -1 0;-1 0 2 -1;0 0 -1 1;b=1 0 0 1x=zhuiganfa(A,b)b = 1 0 0 1x = 1 1 1 2第5题 A=1.5 3 -0.8 4;2 0 9 10;-7 4.8 -0.6 1;14 12.3 -4 5;b=4 0 1 -2;L,U,x,y,index=Doolittle(A,b)L = 1.0000 0 0 0 1.3333 1.0000 0 0 -4.6667 -4.7000 1.0000 0 9.3333 3.9250 -0.8386 1.0000U = 1.5000