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

数值计算方法上机作业热能工程1．设105dxxxInn，（1）由递推公式nIInn151，从0I的几个近似值出发，计算20I；解：易得：0Iln6-ln5=0.1823,程序为：I=0.182;forn=1:20I=(-5)*I+1/n;endI输出结果为：20I=-3.0666e+010（2）粗糙估计20I，用nIInn515111，计算0I；因为0095.0560079.01020201020dxxIdxx所以取0087.0)0095.00079.0(2120I程序为：I=0.0087;forn=1:20I=(-1/5)*I+1/(5*n);endI0I=0.0083（3）分析结果的可靠性及产生此现象的原因(重点分析原因)。首先分析两种递推式的误差；设第一递推式中开始时的误差为000IIE，递推过程的舍入误差不计。并记nnnIIE，则有01)5(5EEEnnn。因为20E20020)5(IE，所此递推式不可靠。而在第二种递推式中nnEEE)51(5110，误差在缩小，所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。2．求方程0210xex的近似根，要求41105kkxx，并比较计算量。（1）在[0，1]上用二分法；程序：a=0;b=1.0;whileabs(b-a)5*1e-4c=(b+a)/2;数值计算方法上机作业热能工程ifexp(c)+10*c-20b=c;elsea=c;endendc结果：c=0.0903（2）取初值00x，并用迭代1021xkex；程序：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;whileabs(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）取初值00x，并用牛顿迭代法；程序：x=0;a=0;b=1;whileabs(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.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）设y=f(x)具有指数形式xbaey（a0,b0）。对此式两边取对数，得xbay1lnln。记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下x23456789t=1/x0.50000.33330.25000.20000.16670.14290.12500.1111z=lny1.85942.10412.25972.25132.27212.30262.29562.3016101112131415160.10000.09090.08330.07690.07140.06670.06252.35042.35992.36092.37952.36092.38882.3758程序：t=[0.50000.33330.25000.20000.16670.14290.12500.11110.10000.09090.08330.07690.07140.06670.0625];z=[1.85942.10412.25972.25132.27212.30262.29562.30162.35042.35992.36092.37952.36092.38882.3758];polyfit(t,z,1)结果：ans=-1.11072.4578由此可得A=2.4578，B=-1.1107，6791.11Aea，b=B=-1.1107方程即为xey1107.16791.11数值计算方法上机作业热能工程计算均方差编程：x=[2:16];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];f(x)=11.6791*exp(-1.1107./x);c=0;fori=1:15a=y(i);b=x(i);c=c+(a-f(b))^2;endaverge=c/15结果：averge=0.05944．设410100141010014101101410010141001014A，625250b，bxA分析下列迭代法的收敛性，并求42110kkxx的近似解及相应的迭代次数。（1）JACOBI迭代；程序：functiony=jacobi(a,b,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;whilenorm(y-x0)1e-4x0=y;y=B*x0+f;n=n+1;endyn以文件名jacobi.m保存。程序：a=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[05-25-26]';x0=[000000]';jacobi(a,b,x0);数值计算方法上机作业热能工程运行结果为:y=1.00002.00001.00002.00001.00002.0000n=28（2）GAUSS-SEIDEL迭代；程序：functiony=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;whilenorm(y-x0)10^(-4)x0=y;y=G*x0+f;n=n+1;endyn以文件名deisel.m保存。程序：a=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[05-25-26]';x0=[000000]';jacobi(a,b,x0);运行结果为：y=1.00002.00001.00002.00001.00002.0000数值计算方法上机作业热能工程n=15（3）SOR迭代（95.0,95.1,334.1）。程序：functiony=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;whilenorm(y-x0)10^(-4)x0=y;y=lw*x0+f;n=n+1;endyn以文件名sor.m保存。程序：a=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[05-25-26]';x0=[000000]';c=[1.3341.950.95];fori=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．用逆幂迭代法求111123136A最接近于11的特征值和特征向量，准确到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;whileabs(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;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方法求解初值问题（1）xxyyyxyyysin2cos22sin22212211，]10,0[x，3)0(2)0(21yy；程序：functionydot=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.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.8389（2）xxyyyxyyysin999cos999999998sin22212211，]10,0[x，3)0(2)0(21yy。和精确解xexyxexyxxcos2)(sin2)(21比较，分析结论。程序：functionydot=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)y2\n')forj=1:length(y)fprintf('%4d%.4f%.4f%.4f%.4f\n',x(j),y(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.9106-0.910910-0.5439-0.5439-0.8389-0.8390结论：R-K方法求解的结果精度较高。7．用有限差分法求解边值问题（h=0.1）：数值计算方法上机作业热能工程1)1()1(0)1(2yyyxy.程序为：h=0.1;n=(1-(-1))/h+1;x(1)=-1;x(n-1