现代数值计算上机实验报告

1、太原科技大学现代数值计算方法上机报 告院系： 华科学院 专业年级： 计算机科学与技术 学生姓名： 张栩嘉 学 号： 201522030129 指导教师： 2017年5月12日数值计算方法上机实习题1 设，（1） 由递推公式，从的几个近似值出发，计算；（2） 粗糙估计，用，计算；（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。%上机习题1%（1）I=0.182;for n=1:20 I=(-5)*I+1/n;endfprintf(I20 的值是 %en,I);%(2)I=0.0087;for n=1:20 I=(-1/5)*I+1/(5*n);endfprintf(I0 的值是 %f

2、n,I);3） 现象产生的原因：假设Sn的真值为Sn*，误差为n，即n=Sn*-Sn 对于真值，我们也有关系式Sn*+5Sn-1*=1n 综合两个递推公式有n=-5n-1 这就意味着哪怕开始只有一点点误差，只要n足够大，按照这种每计算一步误差增长5倍的方式，所得结果总是不可信的。因此整个算法是不稳定的。对于第二种方法 n=-15n-1 误差会以每计算一步缩小到1/5的方式进行，所以以这种方法计算的结果是可靠的，整个算法是稳定的。2. 求方程的近似根，要求，并比较计算量。（1） 在0，1上用二分法；（2） 取初值，并用迭代；（3） 加速迭代的结果；（4） 取初值，并用牛顿迭代法；（5） 分析绝对

3、误差。%（1）在0，1上用二分法；a=0;b=1.0;ci=0;while abs(b-a)5*1e-4 c=(b+a)/2; ci=ci+1; if exp(c)+10*c-20 b=c; else a=c; end endfprintf(二分法求得 c 的值是 %ft,c);fprintf(迭代次数是 %dn,ci);%（2）不动点迭代法，x=0;a=1;ci=0;while abs(x-a)5*1e-4 a=x; x=(2-exp(x)/10; ci=ci+1;endfprintf(不动点迭代求得 x 的值是 %ft,x);fprintf(迭代次数是 %dn,ci);%（3）艾特肯加速迭

4、代x=0;a=0;b=1;ci=0;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; ci=ci+1;endfprintf(艾特肯加速迭代求得 x 的值是 %ft,x);fprintf(迭代次数是 %dn,ci);%（4）用牛顿迭代法x=0;a=0;b=1;ci=0;while abs(b-a)5*1e-4 a=x; x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x; ci=ci+1;endfprintf(牛顿迭代法求得 x 的值是 %ft,x);

5、fprintf(迭代次数是 %dn,ci);%（5）分析绝对误差solve(exp(x)+10*x-2=0);fprintf(x的真值是%.8f,x);运行结果：二分法求得 c 的值是 0.090332迭代次数是 11绝对为误差为0.000193不动点迭代求得 x 的值是 0.090513迭代次数是 4绝对为误差为0.000012艾特肯加速迭代求得 x 的值是 0.099488迭代次数是 3绝对为误差为0.008963牛顿迭代法求得 x 的值是 0.090525迭代次数是 2绝对为误差为0.000000x的真值是0.090525分析可知加速迭代绝对误差最大，二分法次之，牛顿法和不动点法迭代效果

6、最好。3钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系，计算均方差。（注：增速减少，用何种模型）解：设y=f(x)具有指数形式（a0,b1e-4 x0=y; y=B*x0+f;n=n+1;endynb) %高斯-赛德尔迭代function y=gs(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+

7、f;n=1;while norm(y-x0)10(-4) x0=y; y=G*x0+f;n=n+1;endync)%SOR迭代(omega=1.334,1.95,0.95)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;endynA=4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1

8、 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;fprintf(雅可比迭代的结果及迭代次数是：n);jacobi(A,b,x0);fprintf(高斯赛德尔迭代的结果及迭代次数是：n);gs(A,b,x0);c=1.334 1.95 0.95;for i=1:3 w=c(i); fprintf(omega的值为:%.3f 时，SOR迭代的结果及迭代次数是：n,w)sor(A,b,w,x0);end运行结果：雅可比迭代的结果及迭代次数是：y = 1.0000 2.0000 1.0000 2.0000

9、 1.0000 2.0000n = 28高斯赛德尔迭代的结果及迭代次数是：y = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000n = 15omega的值为:1.334 时，SOR迭代的结果及迭代次数是：y = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000n = 13omega的值为:1.950 时，SOR迭代的结果及迭代次数是：y = 1.0000 2.0000 1.0000 2.0000 1.0000 2.0000n = 241omega的值为:0.950 时，SOR迭代的结果及迭代次数是：y = 1.0000 2.

10、0000 1.0000 2.0000 1.0000 2.0000n = 175用逆幂迭代法求最接近于11的特征值和特征向量，准确到。程序：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;end mt=1/norm(v,inf)+pmy=uA=6 3 1;3 2 1;1 1 1; maxtr(A,11,0.001

11、);运行结果：mt = 11.0919my =0.3845 -1.0000 0.73066用经典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) 以文件名lorenzeq.m保存。主窗口输入：x,y=ode45(lorenzeq,0:10,2;3)运行结果为：x = 0 1 2 3 4 5 6 7 8 9 10y = 2.0000 3.0000 1.5775 1.2758 1.1802 -0.14

12、57 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）：程序：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保存。程序：x=0:10;y1=2*exp(-x)+sin(x);y2=2*exp(-x)+cos(x);x,y=ode

13、45(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.720

14、2 -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结论：R-K方法求解的结果精度较高。7用有限差分法求解边值问题（h=0.1）：.程序为：h=0.1;n=(1-(-1)/h+1;x(1)=-1;x(n-1)=1;y(1)=1;y(

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

16、420 0.6752 0.7167 0.7680 0.83080.90728P100，例4.4.6用函数y=asinbx拟合数据x0.10.20.30.40.50.60.70.8y0.61.11.61.82.01.91.71.3解：使用非线性拟合函数nlinfit，编写MATLAB程序如下：functionerr,a,b=nlfitb(x,y) if nargin2, x=1:8/10; y=0.6 1.1 1.6 1.8 2.0 1.9 1.7 1.3; end beta0=1 1; beta=nlinfit(x,y,mymodel,beta0); fprintf(The nonlinear

17、 least square fitting y=a*sin(b*x) for datann); fprintf(%6.1f,x); fprintf(n); fprintf(%6.1f,y); fprintf(nn isnt y=%7.4f*sin(%7.4f*x)nn,beta); z=linspace(x(1),x(end),100); plot(x,y,ro,z,beta(1)*sin(beta(2)*z),b-.) grid on; function yb=mymodel(beta,xb) yb=beta(1)*sin(beta(2)*xb);拟合结果为：The nonlinear le

18、ast 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.9750*sin( 3.0249*x)拟合曲线为：9拟合形如f（x）（a+bx）/(1+cx)的函数的一种快速方法是将最小二乘法用于下列问题：f（x）(1+cx)（a+bx），试用这一方法拟合表4-4给出的中国人口数据。表4-4次序 年份 人口（亿）第一次 1953 5.82第二次 1964 6.95第三次 1982 10.08第四次 1900 11.34第五次 2000

19、 12.66解：把f（x）(1+cx)a+bx进一步变形为f（x）a+bx-cx f（x），这个方程可以看成基函数是1，x，-xf(x)而数据为（xi,f(xi)）的最小二乘拟合问题。编程如下：x=1953 1964 1982 1900 2000;y=5.82 6.95 10.08 11.34 12.66;A=ones(5,1) x -x.*y;Z=Ay;a=Z(1)b=Z(2)c=Z(3)v=linspace(1953,2000,100);plot(x,y,b-+,v,(a+b*v)./(1+c*v),k-);求得的解为a =11.5250 b =-0.0059 c =-5.0979e-0410. 已知20世纪美国人口的统计数字如表（单位：百万）：年份1900191019201930194019501960197019801990人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4试拟合美国人口20世纪的增长率。 t=1900:10:1990;x=76.092.0106