上海电力学院数值分析上机实验报告题目：数值分析上机实验报告学生姓名：11111111111学号：111111*********专业：11112013年12月30日数值计算方法上机实习题1． 设⎰+=105dx xx I nn ， （1） 由递推公式n I I n n 151+-=-，从0I 的几个近似值出发，计算20I ； （2） 粗糙估计20I ，用nI I n n 51511+-=-，计算0I ；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。

(1) 解答：n=0，0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=⎰⎰⎰x d xdx x dx x x I nn这里可以用for 循环，while 循环，根据个人喜好与习惯：for 循环程序： While 循环程序： I=0.1823; I=0.1823; for n=1:20 i=1;I=(-5)*I+1/n; while i<21 End I=(-5)*I+1/i; I i=i+1; fprintf('I20=%f',I) end I = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009 (2) 粗略估计I 20： Mathcad 计算结果： for 循环程序： While 循环程序： >> I=0.007998; I=0.007998; >> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n); >> I n=n+1; I =0.0083 end >> II =0.0083(3) 算法误差分析：计算在递推过程中传递截断误差和舍入误差 第一种算法：（从1——>20）*000e I I=-***21111120115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===1x x 205x +⎛⎜⎜⎜⎠d 7.998103-⨯=误差放大了5n 倍，算法稳定性很不好； 第二种算法：（从20——>1）*n n ne I I =-***111111111()()555555n n n n n n nn e I I I I I I e n n ---=-=-+--+=-=0111...()55nne e e ===误差在逐步缩小，算法趋近稳定，收敛。

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

（1） 在[0，1]上用二分法；function [t i]=erfenfa(a,b) f=@(x)( exp(x)+10*x-2) t=(a+b)./2; i=0;while abs(f(t))>0.001 if f(a)*f(t)<0 b=t;t=(a+b)/2; elseif f(b)*f(t)<0 a=t;t=(b+a)/2; end i=i+1; end 结果： t = 0.0906 i = 11（2） 取初值00=x ，并用迭代1021x k e x -=+；function x=diedai(x0) %x0初值 x=x0;for i=1:10000y=(2-exp(x))./10;x=y;y=(2-exp(x))./10; if abs(x-y)<5*10^(-4) disp('迭代次数'); 2*i break; end end结果：ans =6 x =0.090639135859584（3） 加速迭代的结果(艾特肯Aitken 加速方法)；function [y m]=aitken(a) func=@(x)(( 2-exp(x))./10) x(1)=a; wucha=1;m=1;while wucha> 5*10^(-4) p(m+1)=func(x(m));q(m+1)=func(p(m+1));x(m+1)=q(m+1)-((q(m+1)-p(m+1))^2)./(q(m+1)-2*p(m+1)+x(m));wucha=abs(x(m+1)-x(m)); m=m+1; if m>1000 break; end endformat long y=x(m-1); m=m-1;运行结果y =0.090483741803596m =2（4） 取初值00 x ，并用牛顿迭代法；function x=newtondiedai(x0) x=x0; for i=1:100y=x-(exp(x)+10*x-2)./(exp(x)+10);x=y; y=x-(exp(x)+10*x-2)./(exp(x)+10); if abs(x-y)<0.0001 disp('迭代次数'); i break; end end 运行结果 迭代次数i = 2 x =0.0905（5） 分析绝对误差。

通过指令求得方程精确解的近似解：>> solve('exp(x)+10*x-2=0') ans =-lambertw(1/10*exp(1/5))+1/5 >> format long>> -lambertw(1/10*exp(1/5))+1/5 ans =0.09052510130725 小结：所以我们可以看到，在要求同一运算精度的情况下，采用二分法运算了11次，采用题设的迭代方法运算了6次，采用加速迭代法只运算了2次，采用牛顿迭代法需运算2次。

也就是说加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。

而二分法收敛速度较慢,所以在工程上不经常使用。

3．钢水包使用次数多以后，钢包的容积增大，数据如下： x 2345678910111213141516y6.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试从中找出使用次数和容积之间的关系，计算均方差。

（注：增速减少，用何种模型） 解答：因为不知道其函数形式，可先plot 而后确定使用哪种模型比较合适。

函数图形程序：x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; plot(x,y ，‘b*’)与指数函数趋势类似（但是趋势不同，一个趋于递减，一个趋于递增，使用倒数），故拟合模型为： bxy ae=两边同时取对数：ln ln by ax=+用此模型进行数据拟合：1,ln n n y x ⎛⎫ ⎪⎝⎭x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; t=1./x;s=polyfit(t,log(y),1) s =-1.1107 2.4578即是： 2.457811.679a e == 1.1107b =- 最终拟合曲线为：1.110711.679xy e-=将此拟合曲线和源数据进行对比： 主程序：x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76]; plot(x,y,'ro') hold onlims1=[2,16];fplot('11.679*exp(-1.1107/x)',lims1) hold off均方差的求解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.679*exp(-1.1107./x); c=0;for i=1:15 a=y(i); b=x(i);c=c+(a-f(b))^2;endaverge=c/15averge = sqrt(averge) 0.24384．设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛----------------=410100141010014101101410010141001014A ，⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=625250b ，b x =A 分析下列迭代法的收敛性，并求42110-+≤-kk x x 的近似解及相应的迭代次数。

1） JACOBI 迭代；2） GAUSS-SEIDEL 迭代； 3） SOR 迭代（95.0,95.1,334.1=ω）。

解答：（1） JACOBI 迭代； 定义函数。

function tx=jacobi(A,b,imax,x0,tol) %jacobi 迭代 %初始值x0,次数imax,精度tol del=10^-10;tx=[x0];n=length(x0); for i=1:ndg=A(i,i); if abs(dg)<deldisp('对角元太小'); return end endfor k=1:imax %jacobi 循环 for i=1:n sm=b(i); for j=1:n if j~=ism=sm-A(i,j)*x0(j); end endx(i)=sm/A(i,i); endtx=[tx;x];if norm(x-x0)<tol return elsex0=x;endend保存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];imax=100;tol=10^-4;tx=jacobi(A,b,imax,x0,tol);for j=1:size(tx,1)fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6));end运行结果：1 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000002 0.000000 1.250000 -0.500000 1.250000 -0.500000 1.5000003 0.625000 1.000000 0.500000 1.000000 0.500000 1.2500004 0.500000 1.656250 0.312500 1.656250 0.312500 1.7500005 0.828125 1.531250 0.765625 1.531250 0.765625 1.6562506 0.765625 1.839844 0.679688 1.839844 0.679688 1.8828137 0.919922 1.781250 0.890625 1.781250 0.890625 1.8398448 0.890625 1.925293 0.850586 1.925293 0.850586 1.9453139 0.962646 1.897949 0.948975 1.897949 0.948975 1.92529310 0.948975 1.965149 0.930298 1.965149 0.930298 1.97448711 0.982574 1.952393 0.976196 1.952393 0.976196 1.96514912 0.976196 1.983742 0.967484 1.983742 0.967484 1.98809813 0.991871 1.977791 0.988895 1.977791 0.988895 1.98374214 0.988895 1.992415 0.984831 1.992415 0.984831 1.99444815 0.996208 1.989639 0.994820 1.989639 0.994820 1.99241516 0.994820 1.996462 0.992923 1.996462 0.992923 1.99741017 0.998231 1.995167 0.997583 1.995167 0.997583 1.99646218 0.997583 1.998349 0.996699 1.998349 0.996699 1.99879219 0.999175 1.997745 0.998873 1.997745 0.998873 1.99834920 0.998873 1.999230 0.998460 1.999230 0.998460 1.99943621 0.999615 1.998948 0.999474 1.998948 0.999474 1.99923022 0.999474 1.999641 0.999282 1.999641 0.999282 1.99973723 0.999820 1.999509 0.999755 1.999509 0.999755 1.99964124 0.999755 1.999832 0.999665 1.999832 0.999665 1.99987725 0.999916 1.999771 0.999886 1.999771 0.999886 1.99983226 0.999886 1.999922 0.999844 1.999922 0.999844 1.99994327 0.999961 1.999893 0.999947 1.999893 0.999947 1.99992228 0.999947 1.999964 0.999927 1.999964 0.999927 1.99997329 0.999982 1.999950 0.999975 1.999950 0.999975 1.999964（2）GAUSS-SEIDEL迭代；function tx=gseidel(A,b,imax,x0,tol) %gseidel迭代 %初始值x0,次数imax,精度toldel=10^-10;tx=[x0];n=length(x0);for i=1:ndg=A(i,i);if abs(dg)<deldisp('对角元太小');returnendendfor k=1:imax %G-S循环x=x0;for i=1:nsm=b(i);for j=1:nif j~=ism=sm-A(i,j)*x(j);endendx(i)=sm/A(i,i);endtx=[tx;x];if norm(x-x0)<tolreturnelsex0=x;endend主程序：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];imax=100;tol=10^-4;tx=gseidel(A,b,imax,x0,tol);for j=1:size(tx,1)fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5 ),tx(j,6))endω）。