计算方法上机答案.doc

《计算方法上机答案.doc》由会员分享，可在线阅读，更多相关《计算方法上机答案.doc（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、上海电力学院数值分析上机实验报告 题目： 数值分析上机实验报告 学生姓名： 学 号： 1111 专 业： 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=-. I = -2.0558e+009(2) 粗略估计

2、I20：Mathcad计算结果：for循环程序： While循环程序： I=0.; I=0.; 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

3、)./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)-

4、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.3596m =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

5、=-lambertw(1/10*exp(1/5)+1/5 format long -lambertw(1/10*exp(1/5)+1/5ans = 0.725小结：所以我们可以看到，在要求同一运算精度的情况下，采用二分法运算了11次，采用题设的迭代方法运算了6次，采用加速迭代法只运算了2次，采用牛顿迭代法需运算2次。也就是说加速迭代速度都是超线性收敛的，而简单迭代法是线性收敛的。而二分法收敛速度较慢,所以在工程上不经常使用。3钢水包使用次数多以后，钢包的容积增大，数据如下：x2345678910111213141516y6.428.29.589.59.7109.939.9910.4910.59

6、10.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

7、 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;pl

8、ot(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设，分析下列迭代法的收敛性，并求的近似解及相应的迭

9、代次数。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);

10、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

11、(j,5),tx(j,6);end运行结果：1 0. 0. 0. 0. 0. 0.2 0. 1. -0. 1. -0. 1.3 0. 1. 0. 1. 0. 1.4 0. 1. 0. 1. 0. 1.5 0. 1. 0. 1. 0. 1.6 0. 1. 0. 1. 0. 1.7 0. 1. 0. 1. 0. 1.8 0. 1. 0. 1. 0. 1.9 0. 1. 0. 1. 0. 1.10 0. 1. 0. 1. 0. 1.11 0. 1. 0. 1. 0. 1.12 0. 1. 0. 1. 0. 1.13 0. 1. 0. 1. 0. 1.14 0. 1. 0. 1. 0. 1.15 0.

12、 1. 0. 1. 0. 1.16 0. 1. 0. 1. 0. 1.17 0. 1. 0. 1. 0. 1.18 0. 1. 0. 1. 0. 1.19 0. 1. 0. 1. 0. 1.20 0. 1. 0. 1. 0. 1.21 0. 1. 0. 1. 0. 1.22 0. 1. 0. 1. 0. 1.23 0. 1. 0. 1. 0. 1.24 0. 1. 0. 1. 0. 1.25 0. 1. 0. 1. 0. 1.26 0. 1. 0. 1. 0. 1.27 0. 1. 0. 1. 0. 1.28 0. 1. 0. 1. 0. 1.29 0. 1. 0. 1. 0. 1.（2）GA

13、USS-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-x

14、0)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)e

15、nd（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 t

16、x=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*ex

17、p(-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.24

18、07 -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*

19、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.5

20、772 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

21、-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

22、/(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.

23、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,利用fm

24、insearch寻找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.拟合形如的函数的一种快速方法是将最小二乘法用于下列问题：，试用这一方法拟合下面的表格给出的中国人口数据。中国人口数据 次序 年份 人口（亿） 第一次 195

25、3 5.82 第二次 1964 6.95 第三次 1982 10.08 第四次 1990 11.34 第五次 2000 12.66x=1953 1964 1982 1990 2000; y=5.82 6.95 10.08 11.34 12.66; x1=ones(5,1); A=x1 x x.*y; z=Ayz = 2.9456 -0.00140.0005 因此，方程为： y=( 2.9456-0.0014*x)/(1+0.0005*x)10、已知20世纪美国人口的统计数字如表所示（单位：百万）表 美国人口统计数据年份 1900 1910 1920 1930 1940 1950 1960 19

26、70 1980 1990人口 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4解:由logisitc模型，。其中r为固有增长率，M为最大人口容量。解得。由程序：t=1900:10:1990;p=76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4;n=10;dp=p(2:n)-p(1:n-1);dt=t(2:n)-t(1:n-1);dpdt=(dp./dt)./p(2:n);f=polyfit(p(2:n),dpdt,1);disp(求人口增长率r：),r=f(2)disp(求最大人口容量M：),s=-f(1);M=r/sp1=M./(1+(M/p(1)-1).*exp(-r.*(t-1900);plot(t,p,*,t,p1)运行结果为：求人口增长率r：r = 0.0164求最大人口容量M：M =659.8734美国人口增长模型综上所诉：本题利用logisitc模型拟合，先拟合。人口的长期预测是一件复杂而又困难的事情,但中短期的预报却是现实可能的,。人口增长都不是线性关系,美国人口在局部时段是指数增长的, 但长期趋势却是Logistic 过程。