(完整版)数值计算方法上机实习题答案

1、数值计算方法上机作业 热能工程 1 (1) 由递推公式In 5In 1 ，从I。的几个近似值出发，计算 120 ； n 解：易得：I。ln 6-1 n5=0.1823, 程序为： 1=0.182; for n=1:20 l=(-5)*l+1/n; end I 输出结果为：I 20 = -3.0666e+010 1 所以取 I20 丄(0.0079 0.0095) 0.0087 2 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5* n); end I I 0= 0.0083 (3) 分析结果的可靠性及产生此现象的原因 (重点分析原因)。 首先分析两种递推式的误差

2、；设第一递推式中开始时的误差为 并记 En In I n，则有 E* 5E.1 (5)0。因为 E20 20 1 (5)20E 120，所此递推式不可靠。而在第二种递推式中 E0 -E1 5 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中, 即算法是否数值稳定。 x 4 2.求方程e 10 x 2 0的近似根，要求 Xk1 Xk 5 10 ，并比较计算量。 (1) 在0 , 1上用二分法； 设In 05 xdX (2) 粗糙估计I 20,用ln1 5In 1 丄，计算i 5n 因为 0.0079 20 1 x dx 0 6 20 20 dx 0 5 0.0095 E。 I

3、。 I。，递推过程的舍入误差不计。 1 (-)nEn，误差在缩小, 5 考虑误差是否得到控制， 数值计算方法上机作业 热能工程 程序：a=0;b=1.0; while abs(b-a)5*1e-4 c=(b+a)/2;数值计算方法上机作业 热能工程 if exp(c)+10*c-20 b=c; else a=c; end end c 结果：c = 0.0903 程序：x=0; a=1; while abs(x-a)5*1e-4 a=x; x=(2-exp(x)/10; end x 结果：x = 0.0905 (3) 加速迭代的结果； 程序：x=0; a=0;b=1; while abs(b-a

4、)5*1e-4 a=x; y=exp(x)+10*x-2; z=exp(y)+10*y-2; x=x-(y-x)A2/(z-2*y+x); b=x; end x 结果：x = 0.0995 (4) 取初值X0 0，并用牛顿迭代法; 程序：x=0; a=0;b=1; while abs(b-a)5*1e-4 a=x;(2) 取初值X。 0，并用迭代Xk 1 X 2 e 10 数值计算方法上机作业 热能工程 x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; end x 结果： x = 0.0905 (5) 分析绝对误差。 solve(exp(x)+10*x-2=0) 3.

5、钢水包使用次数多以后，钢包的容积增大，数据如下: x 2 3 4 5 6 7 8 9 y 6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10 11 12 13 14 15 16 10.49 10.59 10.60 10.8 10.6 10.9 10.76 试从中找出使用次数和容积之间的关系，计算均方差。 (注：增速减少，用何种模型) b - 1 设y=f(x)具有指数形式 y ae (a0,b1e-4 x0=y; y=B*x0+f; n=n+1; end y n 以文件名jacobi.m保存。 程序： a=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -

6、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); 运行结果为 y = 1.0000 end averge=c/15 结果：averge = 0.0594 4 1 0 4设A 1 0 0 1 0 1 4 1 0 1 4 1 0 1 4 1 0 1 0 1 0 0 0 1 0 0 1 ， 1 0 4 1 1 4 0 5 2 ，Ax b 5 2 6 分析下列迭代法的收敛性，并求 数值计算方法上机作业 热能工程 2.0000 1.0000 2

7、.0000 1.0000 2.0000 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)10A(-4) x0=y; y=G*x0+f; n=n+1; end y n 以文件名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 -1

8、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.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n = 15 数值计算方法上机作业 热能工程 (3) SOR 迭代( 1.334, 1.95, 0.95 )。 程序： function y=sor(a,b,w,xO) 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;

9、 n=1; while norm(y-x0)10A(-4) x0=y; y=lw*x0+f; n=n +1; end y n 以文件名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:3 w=c(i); sor(a,b,w,x0); end 运行结果分别为： y = 1.0000 2.0000 1.0000 2.0000

10、 1.0000 2.0000数值计算方法上机作业 热能工程 function mt,my=maxtr(A,p,ep) n=len gth(A); B=A-p*eye( n); v0=ones(n ,1); k=1; v=B*v0; while abs (n orm(v,i nf)-norm(v0,i nf)ep %no rm(v-v0)ep k=k+1;13 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 n = 241 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000 6 3 1 5 用逆幕迭代法求 A 3 2 1 取接近

11、于 1 1 1 17 程序: 11的特征值和特征向量，准确到 10 3。 数值计算方法上机作业 热能工程 q=v; u=v/ norm(v,i nf) v=B*u; v0=q; end mt=1/norm(v,i nf)+p my=u 主界面中输入： A=1 -2 -3; maxtr(A,11,0.001) 结果为： 特征值： mt = 11.0919 特征向量: my = 0.3845 -1.0000 0.7306 6.用经典R-K方法求解初值问题 程序：function ydot=lore nzeq(x,y) ydot=_2*y(1)+y (2)+2*si n(x);y(1)-2*y(2)

12、+2*cos(x)-2*si n(x) 以文件民 lorenzeq.m 保存。 主窗口输入： x,y=ode45(lore nzeq,0:10,2;3) 运行结果为: 0 1 2 3 4 5 6 7 8 9 10(1) Y1 Y2 2y1 y2 2sinx 0,10, %（o） 数值计算方法上机作业 热能工程 和精确解yi (X) 2e x SinX比较，分析结论。 y2(x) 2e x cosx 程序：function ydot=lore nzeq1(x,y) ydot=-2*y(1)+y (2)+2*si n(x);998*y(1)-999*y (2)+999*cos(x)-999*si

13、n(x); 以文件名lorenzeq1.m保存。 程序：x=0:10; y1=2*exp(-x)+si n(x); y2=2*exp(-x)+cos(x); x,y=ode45(lore nzeq1,0:10,2;3); fprintf( X y(1) y1 y(2) y2n) for j=1:le ngth(y) fprin tf(%4d %.4f %.4f %.4f %.4fn,x(j),y(j,1),y1(j),y(j,2),y2(j) end 运行结果为： X y(1) y1 y(2) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5772 1.5772 1

14、.2759 1.2761 2 1.1800 1.1800 -0.1455 -0.1455 3 0.2407 0.2407 -0.8904 -0.8904 4 -0.7202 -0.7202 -0.6169 -0.6170 5 -0.9454 -0.9454 0.2972 0.2971 6 -0.2745 -0.2745 0.9648 0.9651 7 0.6588 0.6588 0.7554 0.7557 8 0.9900 0.9900 -0.1448 -0.1448 9 0.4124 0.4124 -0.9106 -0.9109 10 -0.5439 -0.5439 -0.8389 -0.8

15、390 结论：R-K方法求解的结果精度较高。 7.用有限差分法求解边值问题(h=0.1 ):2.000 3.0000 1.5775 1.2758 1.1802 -0.1457 0.2406 -0.8903 -0.7202 -0.6170 -0.9454 0.2971 -0.2745 0.9652 0.6589 0.7557 0.9901 -0.1449 0.4124 -0.9109 -0.5440 -0.8389 y1 2y1 y2 2si nx (2) y2 998y1 999y2 999cosx ,X 0,10 , y1(0) 2 999si nx y2(0) 3 数值计算方法上机作业 热

16、能工程 y (1 x2)y 0 y( 1) y(i) 1 程序为： h=0.1; n=(1-(-1)/h+1; x(1)=-1;x( n-1)=1; y(i)=i;y( n-i)=i; for i=1: n-1 x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)A2); B(i)=2/(hA2)+q(i); end for i=1: n-2 C(i)=-1/(hA2); end H=diag(B)+diag(C,1)+diag(C,-1); g(1)=0+1/(hA2); g( n-1)=0+1/(hA2); for i=2: n-2 g(i)=0; end y=Hg 运行结果为：

17、 y = 0.9027 0.8235 0.7592 0.7074 0.6661 0.6338 0.6095 0.5922 0.5814 0.5767 0.5778 0.5846 0.5974 0.6163 0.6420 0.6752 0.7167 0.7680 0.8308 0.9072 数值计算方法上机作业 热能工程 &拟合形如f( x)(a+bx) /(1+cx)的函数的一种快速方法是将最小二乘法用于下列问题： f (x) (1+cx) a+bx)，试用这一方法拟合表 4-4给出的中国人口数据。 表4-4 次序 年份 人口(亿) a) 1953 5.82 b) 1964 6.95 c) 1982 10.08 d) 1900 11.3