1、 数值计算方法上机作业 热能工程1 设 ，105dxInn（1） 由递推公式 ，从 的几个近似值出发，计算 ；nIn10I 20I解：易得： ln6-ln5=0.1823,0I程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： = -3.0666e+01020I（2） 粗糙估计 ，用 ，计算 ；nIIn5110I因为 9. 6 079. 022102dxdx所以取 87.)95.7.(2 I程序为：I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endI= 0.00830（3） 分析结果的可靠性及产生此现象的原因(重点分析原因
2、) 。首先分析两种递推式的误差；设第一递推式中开始时的误差为 ，递推过程的舍入误差不计。00IE并记 ，则有 。因为nnIE01)5(5Ennn 2，所此递推式不可靠。而在第二种递推式中 ，误差在缩2020)5( nE)51(510小，所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。2 求方程 的近似根，要求 ，并比较计算量。021xe 4105kx（1） 在0，1上用二分法；程序：a=0;b=1.0;while abs(b-a)5*1e-4 c=(b+a)/2;数值计算方法上机作业 热能工程if exp(c)+10*c-20b=
3、c;else a=c;end endc结果：c =0.0903（2） 取初值 ，并用迭代 ；0x102xke程序：x=0;a=1;while abs(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*x-2;z=exp(y)+10*y-2;x=x-(y-x)2/(z-2*y+x);b=x;endx结果：x =0.0995（4） 取初值 ，并用牛顿迭代法；0x程序：x=0;a=0;b=1;while abs(b-a)5*1e-4a=
4、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钢水包使用次数多以后，钢包的容积增大，数据如下：x 2 3 4 5 6 7 8 9y 6.42 8.2 9.58 9.5 9.7 10 9.93 9.9910 11 12 13 14 15 1610.49 10.59 10.60 10.8 10.6 10.9 10.76试从中找出使用次数和容积之间的关系，计算均方差。 （注：增速减少，用何种模型）设 y=f(x)具有指数形式 （a0,b1e-
5、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);数值计算方法上机作业 热能工程运行结果为:y =1.00002.00001.00002.00001.00002.0000n =28（2） GAUSS-SEIDEL 迭代；程序：function y=seidel(a,b,x0)D=dia
6、g(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 -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.
7、00001.00002.0000数值计算方法上机作业 热能工程n =15（3） SOR 迭代（ ） 。95.0,.1,34.程序：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)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
8、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.00001.00002.00001.00002.0000n =241y =1.00002.00001.00002.00001.00002.0000n =175用逆幂迭代法求 最接近于 11 的特征值和特征向量
9、，准确到 。1236A 310程序：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=q;end mt=1/norm(v,inf)+pmy=u主界面中输入：A=1 -2 -3; maxtr(A,11,0.001)结果为：特征值：mt =11.0919特征向量：my = 0.3845-1.00000
10、.73066用经典 R-K 方法求解初值问题（1） ， ， ； xysin2co2i12 10,3)(21y程序：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 保存。主窗口输入：x,y=ode45(lorenzeq,0:10,2;3)运行结果为：x =012345678910y =数值计算方法上机作业 热能工程2.0000 3.00001.5775 1.27581.1802 -0.14570.2406 -0.8903-0.7202 -0.61
11、70-0.9454 0.2971-0.2745 0.96520.6589 0.75570.9901 -0.14490.4124 -0.9109-0.5440 -0.8389（2） ， ， 。 xyyxsin9co98sin2121 10,3)(21y和精确解 比较，分析结论。xexxs)(i2程序：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*
12、exp(-x)+cos(x);x,y=ode45(lorenzeq1,0:10,2;3);fprintf( x y(1) y1 y(2) y2n)for j=1:length(y)fprintf(%4d %.4f %.4f %.4f %.4fn,x(j),y(j,1),y1(j),y(j,2),y2(j)end运行结果为:x y(1) y1 y(2) y20 2.0000 2.0000 3.0000 3.00001 1.5772 1.5772 1.2759 1.27612 1.1800 1.1800 -0.1455 -0.14553 0.2407 0.2407 -0.8904 -0.89044
13、 -0.7202 -0.7202 -0.6169 -0.61705 -0.9454 -0.9454 0.2972 0.29716 -0.2745 -0.2745 0.9648 0.96517 0.6588 0.6588 0.7554 0.75578 0.9900 0.9900 -0.1448 -0.14489 0.4124 0.4124 -0.9106 -0.910910 -0.5439 -0.5439 -0.8389 -0.8390结论：R-K 方法求解的结果精度较高。7用有限差分法求解边值问题（h=0.1）：数值计算方法上机作业 热能工程.1)(02yx程序为：h=0.1;n=(1-(-1
14、)/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(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.5
15、9740.61630.64200.67520.71670.76800.8308数值计算方法上机作业 热能工程0.90728拟合形如 f（x）（a+bx ）/(1+cx)的函数的一种快速方法是将最小二乘法用于下列问题： f（x）(1+cx)（a+bx） ，试用这一方法拟合表 4-4 给出的中国人口数据。表 4-4次序 年份 人口（亿）a) 1953 5.82b) 1964 6.95c) 1982 10.08d) 1900 11.34e) 2000 12.66解：把 f（x）(1+cx)（a+bx ）变成 f（x）a+bx-cx f（x）则近似看成基函数是 1，x，-x*f（x）而数据是（x i，f（x i） ）的最小二乘拟合问题，程序如下：x=1953 1964 1982 1900 2000;y=5.82 6.95 10.08 11.34 12.66;A=ones(5,1) x -x.*y;Z=Ay;a=Z(1)b=Z(2)c=Z(3)结果：a =11.5250b =-0.0059c =-5.0979e-004
