计算方法上机答案.doc

键入文字键入文字键入文字上海电力学院数值分析上机实验报告 题目： 数值分析上机实验报告 学生姓名： 11111111111 学 号： 111111111111111 专 业： 1111 2013年12月30日数值计算方法上机实习题1 设，（1） 由递推公式，从的几个近似值出发，计算；（2） 粗糙估计，用，计算；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。(1) 解答：n=0，这里可以用for循环，while循环，根据个人喜好与习惯：for循环程序： While 循环程序：I=0.1823; I=0.1823;for n=1:20 i=1; I=(-5)*I+1/n; while i II20=-2055816073.851284 I = -2.0558e+009(2) 粗略估计I20：Mathcad计算结果：for循环程序： While循环程序： I=0.007998; I=0.007998; for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n I n=n+1;I =0.0083 end II =0.0083(3) 算法误差分析：计算在递推过程中传递截断误差和舍入误差第一种算法：（从120） 误差放大了5n倍，算法稳定性很不好；第二种算法：（从201）误差在逐步缩小，算法趋近稳定，收敛。2 求方程的近似根，要求，并比较计算量。（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.0906i = 11（2） 取初值，并用迭代；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) 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 m1000 break; end end format long y=x(m-1); m=m-1;运行结果y =0.090483741803596m =2（4） 取初值，并用牛顿迭代法；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) 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/5ans = 0.09052510130725小结：所以我们可以看到，在要求同一运算精度的情况下，采用二分法运算了11次，采用题设的迭代方法运算了6次，采用加速迭代法只运算了2次，采用牛顿迭代法需运算2次。也就是说加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。而二分法收敛速度较慢,所以在工程上不经常使用。3钢水包使用次数多以后，钢包的容积增大，数据如下：x2345678910111213141516y6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.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*)与指数函数趋势类似（但是趋势不同，一个趋于递减，一个趋于递增，使用倒数），故拟合模型为： 两边同时取对数： 用此模型进行数据拟合：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即是： 最终拟合曲线为：将此拟合曲线和源数据进行对比：主程序：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设，分析下列迭代法的收敛性，并求的近似解及相应的迭代次数。1） JACOBI迭代；2） GAUSS-SEIDEL迭代；3） SOR迭代（）。解答：（1） JACOBI迭代；定义函数。function tx=jacobi(A,b,imax,x0,tol) %jacobi迭代 %初始值x0,次数imax,精度toldel=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)del disp(对角元太小); return endendfor k=1:imax %jacobi循环 for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x0(j); end end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0)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 %fn,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:n dg=A(i,i); if abs(dg)del disp(对角元太小); return endend for k=1:imax %G-S循环 x=x0; for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); end tx=tx;x; if norm(x-x0)tol return else x0=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 %fn,j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)end（3）SOR迭代（）。function tx=sor(A,b,imax,x0,tol,w) %sor迭代 %初始值x0,次数imax,精度toldel=10-10;tx=x0;n=length(x0);for i=1:n dg=A(i,i); if abs(dg)del disp(对角元太小); return endend for k=1:imax %G-S循环 x=x0; for i=1:n sm=b(i); for j=1:n if j=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i); end tx=tx;x; if norm(x-x0)ep %norm(v-v0)ep k=k+1; q=v; u=v/norm(v,inf) v=B*u; v0=q;end mt=1/norm(v,inf)+pmy=u主程序：A=6 3 1;3 2 1;1 1 1; maxtr(A,11,0.001)6用经典R-K方法求解初值问题（1）， ；（2）， 。和精确解比较，分析结论。解答：（1）主程序： 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);主窗口输入：x=0:10;y1=2*exp(-x)+sin(x);y2=2*exp(-x)+cos(x);x,y=ode45(lorenzeq,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) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5775 1.5772 1.2758 1.2761 2 1.1802 1.1800 -0.1457 -0.1455 3 0.2406 0.2407 -0.8903 -0.8904 4 -0.7202 -0.7202 -0.6170 -0.6170 5 -0.9454 -0.9454 0.2971 0.2971 6 -0.2745 -0.2745 0.9652 0.9651 7 0.6589 0.6588 0.7557 0.7557 8 0.9901 0.9900 -0.1449 -0.1448 9 0.4124 0.4124 -0.9109 -0.9109 10 -0.5440 -0.5439 -0.8389 -0.8390（2）主程序：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);主窗口输入：x=0:10;y1=2*exp(-x)+sin(x);y2=2*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) y2 0 2.0000 2.0000 3.0000 3.0000 1 1.5772 1.5772 1.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.8390小结：四阶RungeKutta方法得到的结果已很接近精确解，证明这种迭代方法精确度很好，是一种有效的算法。7用有限差分法求解边值问题（h=0.1）：.解答：主程序：h=0.1;n=(1-(-1)/h+1;x(1)=-1;x(n-1)=1;y(1)=1;y(n-1)=1;for i=1:n-1 x(i)=x(1)+(i-1)*h; q(i)=(1+x(i)2); B(i)=2/(h2)+q(i);endfor i=1:n-2 C(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-2 g(i)=0;endy=Hg运行结果：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.83080.9072小结：有限差分法是将微分方程离散化为差分方程进行求解，而求解差分方程就是解一个三对角矩阵线性方程组，从以上程序可以看出编程很容易实现用追赶法求解三对角矩阵线性方程组，且具有较高的精度。8用函数拟合数据x0.10.20.30.40.50.60.70.8y0.61.11.61.82.01.91.71.5解：此为一个非线性最小二乘问题，按照最小二乘原理，应选取参数使得表达式达到极小值。用matlab仿真需要借助fminsearch函数。完整的matlab程序如下：解答：非线性最小二乘问题：x=0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8;y=0.6 1.1 1.6 1.8 2 1.9 1.7 1.4;第一步：创建 fitfun函数 fitfun.mfunction err=fitfun(c,x,y)a=c(1);b=c(2);err=y-a*sin(b*x);err=err*err;第二步：建立nlfit.m,利用fminsearch寻找fitfun极小值：function err,a,b=nlfit(x,y)if nargin nlfitThe nonlinear least square fitting y=a*sin(b*x) for data 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3 is y= 1.9751*sin( 3.0250*x)最终输出图如图所示：9.拟合形如的函数的一种快速方法是将