当前位置:文档之家› (完整版)数值计算方法上机实习题答案

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

1. 设⎰+=105dx xx I nn , (1) 由递推公式nI I n n 151+-=-,从0I 的几个近似值出发,计算20I ; 解:易得:0I =ln6-ln5=0.1823, 程序为:I=0.182; for n=1:20I=(-5)*I+1/n; end I输出结果为:20I = -3.0666e+010 (2) 粗糙估计20I ,用nI I n n 515111+-=--,计算0I ; 因为 0095.056 0079.01020201020≈<<≈⎰⎰dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为:I=0.0087; for n=1:20I=(-1/5)*I+1/(5*n); end I0I = 0.0083(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。

首先分析两种递推式的误差;设第一递推式中开始时的误差为000I I E '-=,递推过程的舍入误差不计。

并记nn n I I E '-=,则有01)5(5E E E n n n -==-=- 。

因为=20E 20020)5(I E >>-,所此递推式不可靠。

而在第二种递推式中n nE E E )51(5110-==-= ,误差在缩小,所以此递推式是可靠的。

出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。

2. 求方程0210=-+x e x的近似根,要求41105-+⨯<-k k x x ,并比较计算量。

(1) 在[0,1]上用二分法; 程序:a=0;b=1.0;while abs(b-a)>5*1e-4 c=(b+a)/2;if exp(c)+10*c-2>0 b=c; else a=c; end end c结果:c =0.0903(2) 取初值00=x ,并用迭代1021x k e x -=+;程序: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)>5*1e-4 a=x;y=exp(x)+10*x-2; z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x); b=x; end x结果:x =0.0995(4) 取初值00=x ,并用牛顿迭代法; 程序:x=0; a=0;b=1;while abs(b-a)>5*1e-4 a=x;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.钢水包使用次数多以后,钢包的容积增大,数据如下:试从中找出使用次数和容积之间的关系,计算均方差。

(注:增速减少,用何种模型) 设y=f(x)具有指数形式xbaey =(a>0,b<0)。

对此式两边取对数,得xba y 1ln ln +=。

记A=lna ,B=b ,程序:t=[0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769 0.0714 0.0667 0.0625];z=[1.8594 2.1041 2.2597 2.2513 2.2721 2.3026 2.2956 2.3016 2.3504 2.3599 2.3609 2.3795 2.3609 2.3888 2.3758]; polyfit(t,z,1) 结果:ans = -1.1107 2.4578由此可得 A=2.4578,B=-1.1107,6791.11==Ae a ,b=B=-1.1107方程即为xey 1107.16791.11-=计算均方差编程:x=[2:16]; y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76]; f(x)=11.6791*exp( -1.1107./x); c=0;for i=1:15 a=y(i); b=x(i);c=c+(a-f(b))^2; endaverge=c/15 结果:averge =0.05944.设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛----------------=410100141010014101101410010141001014A ,⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=625250b ,b x =A 分析下列迭代法的收敛性,并求42110-+≤-kk x x 的近似解及相应的迭代次数。

(1) JACOBI 迭代; 程序:function y=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;while norm (y-x0)>1e-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 -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=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 -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 =15(3) SOR 迭代(95.0,95.1,334.1=ω)。

程序: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; 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.00002.0000 1.0000 2.0000 1.0000 2.0000 n =13 y =1.00002.0000 1.0000 2.0000 1.0000 2.0000 n =241 y =1.00002.0000 1.0000 2.0000 1.0000 2.0000 n =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;while abs(norm(v,inf)-norm(v0,inf))>ep %norm(v-v0)>ep k=k+1;q=v;u=v/norm(v,inf) v=B*u; v0=q; endmt=1/norm(v,inf)+p my=u主界面中输入:A=[1 -2 -3]; maxtr(A,11,0.001) 结果为: 特征值: mt =11.0919特征向量: my =0.3845 -1.0000 0.73066.用经典R-K 方法求解初值问题(1)⎩⎨⎧-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211,]10,0[∈x ,⎩⎨⎧==3)0(2)0(21y y ; 程序: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 =0 1 2 3 4 5 6 7 8 9 10 y =2.00003.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 (2)⎩⎨⎧-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212211,]10,0[∈x ,⎩⎨⎧==3)0(2)0(21y y 。

和精确解⎩⎨⎧+=+=--xe x y xe x y xx cos 2)(sin 2)(21比较,分析结论。

程序: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 保存。

相关主题