数值计算方法上机实习题答案

1、数值计算方法上机作业热能工程1 设 I n1x ndx ，0 5x（ 1）由递推公式 I n5I n11 ，从 I 0 的几个近似值出发，计算I 20 ；n解：易得： I 0ln6-ln5=0.1823,程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： I 20 = -3.0666e+010（ 2）粗糙估计 I 20 ，用 I n 11 I n 11，计算 I 0 ；55n0.00791x201 x200.0095因为0dx I200dx65所以取 I 201 (0.00790.0095)0.00872程序为： I=0.0087;for n=1:20I

2、=(-1/5)*I+1/(5*n);endII 0 = 0.0083（ 3）分析结果的可靠性及产生此现象的原因(重点分析原因 )。首先分析两种递推式的误差；设第一递推式中开始时的误差为E0I 0 I 0 ，递推过程的舍入误差不计。并记 EnI nI n ，则有 En5En 1( 5) n E0 。因为 E20( 5) 20 E0I 20 ，所此递推式不可靠。 而在第二种递推式中E01 E1( 1 )n En ，误差在缩小，55所以此递推式是可靠的。 出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。2 求方程 ex10x20 的近似根，要求xk 1xk5

3、 10 4 ，并比较计算量。（ 1） 在 0， 1上用二分法；程序： a=0;b=1.0;while abs(b-a)5*1e-4c=(b+a)/2;数值计算方法上机作业热能工程if exp(c)+10*c-20b=c;else a=c;endendc结果： c =0.0903（ 2） 取初值 x00 ，并用迭代 xk 12 ex；10程序： x=0;a=1;whileabs(x-a)5*1e-4a=x;x=(2-exp(x)/10;endx结果： x =0.0905（ 3） 加速迭代的结果；程序： x=0;a=0;b=1;while abs(b-a)5*1e-4a=x;y=exp(x)+10

4、*x-2;z=exp(y)+10*y-2;x=x-(y-x)2/(z-2*y+x);b=x;endx结果： x =0.0995（ 4）取初值 x00 ，并用牛顿迭代法；程序： x=0;a=0;b=1;while abs(b-a)5*1e-4a=x;数值计算方法上机作业热能工程x=x-(exp(x)+10*x-2)/(exp(x)+10);b=x;endx结果：x =0.0905（ 5）分析绝对误差。solve(exp(x)+10*x-2=0)3钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.99101112131415161

5、0.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）b设 y=f(x) 具有指数形式yaex （a0,b1e-4x0=y;y=B*x0+f;n=n+1;endyn以文件名jacobi.m 保存。程序：a=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4; b=0 5 -2 5 -2 6;x0=0 0 0 0 0 0;jacobi(a,b,x0);数值计算方法上机作业热能工程运行结果为

6、:y =1.00002.00001.00002.00001.00002.0000n =28（ 2） GAUSS-SEIDEL 迭代；程序：function y=seidel(a,b,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);G=(D-L)U;f=(D-L)b;y=G*x0+f;n=1;while norm(y-x0)10(-4)x0=y;y=G*x0+f;n=n+1;endyn以文件名deisel.m 保存。程序：a=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0

7、 -1 4 -1;0 0 -1 0 -1 4; b=0 5 -2 5 -2 6;x0=0 0 0 0 0 0;jacobi(a,b,x0);运行结果为：y =1.00002.00001.00002.00001.00002.0000数值计算方法上机作业热能工程n =15（ 3） SOR 迭代（1.334, 1.95, 0.95 ）。程序：function y=sor(a,b,w,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);lw=(D-w*L)(1-w)*D+w*U);f=(D-w*L)b*w;y=lw*x0+f;n=1;while norm(y-x0

8、)10(-4)x0=y;y=lw*x0+f;n=n+1;endyn以文件名sor.m 保存。程序：a=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4; b=0 5 -2 5 -2 6;x0=0 0 0 0 0 0;c=1.334 1.95 0.95;for i=1:3w=c(i);sor(a,b,w,x0);end运行结果分别为：y =1.00002.00001.00002.00001.00002.0000n =数值计算方法上机作业热能工程13y =1.00002.0

9、0001.00002.00001.00002.0000n =241y =1.00002.00001.00002.00001.00002.0000n =176315用逆幂迭代法求A321 最接近于 11 的特征值和特征向量，准确到10 3 。111程序：function mt,my=maxtr(A,p,ep)n=length(A);B=A-p*eye(n);v0=ones(n,1);k=1;v=B*v0;while abs(norm(v,inf)-norm(v0,inf)ep%norm(v-v0)epk=k+1;数值计算方法上机作业热能工程q=v;u=v/norm(v,inf)v=B*u;v0=

10、q;endmt=1/norm(v,inf)+pmy=u主界面中输入：A=1 -2 -3; maxtr(A,11,0.001)结果为：特征值：mt =11.0919特征向量：my =0.3845-1.00000.73066用经典R-K 方法求解初值问题y12 y1y22sin x， xy1 (0)2（ 1）y12 y22 cosx 2 sin x0,10 ，；y2y2 (0)3程序： function ydot=lorenzeq(x,y)ydot=-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)以文件民lorenzeq.m 保存。主窗口输入：

11、x,y=ode45(lorenzeq,0:10,2;3)运行结果为：x =012345678910y =数值计算方法上机作业热能工程2.00003.00001.57751.27581.1802-0.14570.2406-0.8903-0.7202-0.6170-0.94540.2971-0.27450.96520.65890.75570.9901-0.14490.4124-0.9109-0.5440-0.8389y12y1y22 sin x 0,10 ，y1 (0)2（ 2）998y1， xy2 (0)。y2999y2 999cosx 999 sin x3和精确解y1 ( x)2e xsin

12、x2e x比较，分析结论。y2 ( x)cos x程序： function ydot=lorenzeq1(x,y)ydot=-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x);以文件名lorenzeq1.m 保存。程序： x=0:10;y1=2*exp(-x)+sin(x);y2=2*exp(-x)+cos(x);x,y=ode45(lorenzeq1,0:10,2;3);fprintf(xy(1)y1y(2)y2n)for j=1:length(y)fprintf(%4d%.4f%.4f%.4f%.4fn,x(j),y(

13、j,1),y1(j),y(j,2),y2(j)end运行结果为 :xy(1)y1y(2)y202.00002.00003.00003.000011.57721.57721.27591.276121.18001.1800-0.1455-0.145530.24070.2407-0.8904-0.89044-0.7202-0.7202-0.6169-0.61705-0.9454-0.94540.29720.29716-0.2745-0.27450.96480.965170.65880.65880.75540.755780.99000.9900-0.1448-0.144890.41240.4124-0

14、.9106-0.910910-0.5439-0.5439-0.8389-0.8390结论： R-K 方法求解的结果精度较高。7 用有限差分法求解边值问题（h=0.1 ）：数值计算方法上机作业热能工程y(1x2 ) y0.y( 1)y(1)1程序为：h=0.1;n=(1-(-1)/h+1;x(1)=-1;x(n-1)=1;y(1)=1;y(n-1)=1;for i=1:n-1x(i)=x(1)+(i-1)*h;q(i)=(1+x(i)2);B(i)=2/(h2)+q(i);endfor i=1:n-2C(i)=-1/(h2);endH=diag(B)+diag(C,1)+diag(C,-1);g

15、(1)=0+1/(h2);g(n-1)=0+1/(h2);for i=2:n-2g(i)=0;endy=Hg运行结果为：y =0.90270.82350.75920.70740.66610.63380.60950.59220.58140.57670.57780.58460.59740.61630.64200.67520.71670.76800.8308数值计算方法上机作业热能工程0.90728拟合形如f（ x）（ a+bx）/(1+cx) 的函数的一种快速方法是将最小二乘法用于下列问题：（ a+bx），试用这一方法拟合表4-4 给出的中国人口数据。f（ x）(1+cx)表 4-4次序年份人口（亿）a)19535.82b)19646.95c)198210.08d)190011.34e)200012