数值分析上机题报告

1、数值计算方法上机实习作业上海电力学院数值计算方法上机实习（报告） 院系： 电力与自动化工程学院 专业年级： 电力系统及其自动化一班 学生姓名： 学 号： 指导教师： 年 月 日1 设， （1） 由递推公式，从的几个近似值出发，计算；（2） 粗糙估计，用，计算；（3） 分析结果的可靠性及产生此现象的原因。（1）因为由积分可以得到=1/6，所以取的两个近似值I1=0.1667和I2=0.167程序：I1=0.1667;I2=0.167;for n=1:20 I1=-5*I1+1/n; I2=-5*I2+1/n; fprintf(%.1d %.4f %.4fn,n,I1,I2)end 运行结果：1

2、0.1665 0.16502 -0.3325 -0.32503 1.9958 1.95834 -9.7292 -9.54175 48.8458 47.90836 -244.0625 -239.37507 1220.4554 1197.01798 -6102.1518 -5984.96439 30510.8700 29924.932510 -152554.2502 -149624.562711 762771.3419 748122.904412 -3813856.6262 -3740614.438713 19069283.2078 18703072.270314 -95346415.9675 -

3、93515361.280015 476732079.9042 467576806.466716 -2383660399.4583 -2337884032.270817 11918301997.3502 11689420161.412718 -59591509986.6956 -58447100807.008019 297957549933.5305 292235504035.092420 -1489787749667.6023 -1461177520175.4116（2） 因为 所以取程序：I=0.0095;for n=20:-1:1 I=-1/5*I+1/(5*n); fprintf(%.1

4、d %.4fn,n,I)end运行结果：20 0.008119 0.008918 0.009317 0.009916 0.010515 0.011214 0.012013 0.013012 0.014111 0.015410 0.01699 0.01888 0.02127 0.02436 0.02855 0.03434 0.04313 0.05802 0.08841 0.1823（3）首先分析两种递推式的误差；设第一递推式中开始时的误差为，递推过程的舍入误差不计。并记，则有。因为，所此递推式不可靠。而在第二种递推式中，误差在缩小，所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过

5、程中，应考虑误差是否得到控制，即算法是否数值稳定。2.求方程的近似根，要求误差不超过，并比较计算量。（1） 在0，1上用二分法；（2） 取初值，并用迭代；（3） 取初值，并用牛顿迭代法。x2345678910111213141516y6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76先在命令窗口中输入命令： fplot(exp(x)+10*x-2,0 1);grid得到如下图形，可以看到曲线与X轴的交点，并且能够从图中大致估算到近似根的位置。（1）在0，1上用二分法； 编写Matlab函数M文件bisect如下：functio

6、n x=bisect(fname,a,b,e)if nargin0,error(函数在两端点值必须异号 );endx=(a+b)/2while(b-a)(2*e), fx=feval(fname,x); if fa*fx fname=inline(exp(x)+10*x-2); bisect(fname,0,1,5*10-4)x = 0.5000x = 0.2500x = 0.1250x = 0.0625x = 0.0938x = 0.0781x = 0.0859x = 0.0898x = 0.0918x = 0.0908x = 0.0903ans = 0.0903（2）取初值，并用迭代； 编

7、写Matlab函数M文件diedai如下：k=0;x=0;fprintf(n=%2d x%2d=%f n,k,k,x)for k=1:10 if abs(2-exp(x)/10)5*10-4 break else x=(2-exp(x)/10; fprintf(n=%2d x%2d=%f n,k,k,x) endend在命令窗口中编写内嵌函数表达式，并调用函数M文件diedai：n= 0 x 0=0.000000 n= 1 x 1=0.100000 n= 2 x 2=0.089483 n= 3 x 3=0.090639 n= 4 x 4=0.090513 n= 5 x 5=0.090526 n

8、= 6 x 6=0.090525 n= 7 x 7=0.090525 n= 8 x 8=0.090525 n= 9 x 9=0.090525 n=10 x10=0.090525ans = 0.0905（3）取初值，并用牛顿迭代法。 编写Matlab函数M文件newton如下：function x = newton(fname,dfname,x0,e,N)if nargin 5,N=500;endif nargin e&k fun=inline(exp(x)+10*x-2); dfun=inline(exp(x)+10); newton(fun,dfun,0,0.5e-4)n= 0 x 0= 0

9、.000000000n= 1 x 1= 0.090909091n= 2 x 2= 0.090525109n= 3 x 3= 0.090525101ans = 0.0905小结：我们可以看到，在运算要求到同一精度的情况下，采用（1）的二分法运算了11次，采用（2）的方法运算了6次，采用（3）的加速迭代法运算了3次，采用（4）的牛顿迭代法也需运算3次。3钢水包使用次数多以后，钢包的容积增大，数据如下：x23456789y6.428.29.589.59.7109.939.991011121314151610.4910.5910.6010.810.610.910.76试从中找出使用次数和容积之间的关系

10、，计算均方差。（注：增速减少，用何种模型）（1）由数据点分布图可知，拟合曲线y=f(x)随着x的增加而上升，但上升速度由快到慢，当x趋于无穷大时，y趋于某个常数，故曲线有一水平渐进线。根据上述特征很容易想到用Logistic人口模型来拟合该曲线。设y=f(x)具有指数形式（a0,b0）。对此式两边取对数，得。记A=lna，B=b，并引入新变量z=lny，t=1/x。引入新变量后的数据表如下x23456789t=1/x0.50000.33330.25000.20000.16670.14290.12500.1111z=lny1.85942.10412.25972.25132.27212.30262

11、.29562.3016X10111213141516t=1/x0.10000.09090.08330.07690.07140.06670.0625z=lny2.35042.35992.36092.37952.36092.38882.3758程序：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.359

12、9 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，b=B=-1.1107方程即为 （2）程序：x=2 3 4 5 6 7 8 9 10 11 12 13 14 15 16;y=6.42 8.20 9.58 9.50 9.70 10.00 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76;sum=(11.6791*exp(-1.1107/x(1)-y(1)2;for i=1:14 sum=sum

13、+(11.6791*exp(-1.1107/x(i+1)-y(i+1)2 ; d=sum0.5; endd运行结果：d =0.9441所以均方差为0.9441结论与分析：本题也可用表达式来拟合，也可以得到较小的均方差，但与Logistic相比均方差稍大些。在编译程序中用到polyfit函数拟合，它其实就是利用了最小二乘法的原理。4设，分析下列迭代法的收敛性，并求的近似解及相应的迭代次数。（1） JACOBI迭代；（2） GAUSS-SEIDEL迭代；（3） SOR迭代（）。因为系数矩阵是严格对角占优，所以用JACOBI迭代和GAUSS-SEIDEL迭代都是收敛的；又因系数矩阵对称正定，且(0,

14、2)，所以SOR迭代也收敛。（1）程序：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

15、-x0)tol return else x0=x; end endA=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

16、运行结果：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.83984

17、4 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.9

18、76196 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.9

19、98231 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.999

20、282 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.99994

21、7 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）程序：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

22、%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; endendA=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,

23、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.187500 1.203125 0.113281 1.4814453 0.613281 1.384766 0.517334 1.560974 0.606796 1.7810334 0.736435 1.715

24、141 0.764287 1.776880 0.818263 1.8956385 0.873005 1.863889 0.884102 1.893843 0.913342 1.9493616 0.939433 1.934219 0.944356 1.949283 0.958216 1.9756437 0.970875 1.968362 0.973322 1.975603 0.979902 1.9883068 0.985991 1.984804 0.987178 1.988268 0.990344 1.9943819 0.993268 1.992698 0.993837 1.994362 0.9

25、95360 1.99729910 0.996765 1.996490 0.997038 1.997291 0.997770 1.99870211 0.998445 1.998313 0.998577 1.998698 0.998928 1.99937612 0.999253 1.999189 0.999316 1.999374 0.999485 1.99970013 0.999641 1.999610 0.999671 1.999699 0.999752 1.99985614 0.999827 1.999813 0.999842 1.999855 0.999881 1.99993115 0.9

26、99917 1.999910 0.999924 1.999931 0.999943 1.99996716 0.999960 1.999957 0.999964 1.999967 0.999973 1.999984（3）W=1.334程序：function tx=sor(A,b,imax,x0,tol,w) %sor迭代 %初始值x0,次数imax,精度tol，松弛因子wdel=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 % SOR循环 x=

27、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)tol return else x0=x; end endA=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-

28、4;w=1.334;tx=sor(A,b,imax,x0,tol,w);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.667500 -0.110889 1.630519 0.432889 2.1083873 1.099889 1.584755 1.145478 2.016105 1.092449 2

29、.0431474 0.833524 2.162521 1.025372 1.978394 1.030507 2.0042245 1.102598 1.998570 0.985252 2.046688 1.006313 1.9957766 0.980826 1.991270 1.016176 1.985512 0.988740 2.0030507 0.998661 2.004109 0.992153 1.998020 1.005488 1.9981948 1.001157 1.998227 1.000767 2.003133 0.998018 2.0001989 1.000067 2.00021

30、0 1.000925 1.998623 1.000339 2.00035510 0.999588 2.000214 0.999422 2.000243 1.000158 1.99974111 1.000290 1.999885 1.000149 2.000118 0.999862 2.00009012 0.999904 2.000010 1.000023 1.999890 1.000043 1.99999213 0.999999 2.000018 0.999959 2.000037 1.000001 1.99999014 1.000019 1.999987 1.000018 2.000000

31、0.999992 2.000007W=0.95程序：function tx=sor(A,b,imax,x0,tol,w) %sor迭代 %初始值x0,次数imax,精度tol，松弛因子wdel=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 % SOR循环 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

32、)+(1-w)*x0(i); end tx=tx;x; if norm(x-x0)tol return else x0=x; end endA=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=1000;tol=10-4;w=0.95;tx=sor(A,b,imax,x0,tol,w);for j=1:size(tx,1) fprintf(%d %f %f %f %f %f %fn,j,tx(j,1

33、),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.187500 -0.192969 1.141670 0.078178 1.3977373 0.553178 1.350992 0.439321 1.498869 0.537714 1.7269334 0.704501 1.654414 0.706017 1.725149 0.764678 1.8606375 0.837871 1.818505 0.843820 1

34、.854770 0.877538 1.9268546 0.914297 1.904393 0.917620 1.923734 0.935685 1.9615037 0.954895 1.949667 0.956671 1.959909 0.966165 1.9797498 0.976269 1.973521 0.977213 1.978912 0.982201 1.9893489 0.987516 1.986072 0.988015 1.988907 0.990638 1.99439710 0.993433 1.992674 0.993696 1.994165 0.995076 1.99705

35、311 0.996546 1.996146 0.996684 1.996931 0.997410 1.99845012 0.998183 1.997973 0.998256 1.998386 0.998637 1.99918513 0.999044 1.998934 0.999082 1.999151 0.999283 1.99957114 0.999497 1.999439 0.999517 1.999553 0.999623 1.99977415 0.999736 1.999705 0.999746 1.999765 0.999802 1.99988116 0.999861 1.99984

36、5 0.999866 1.999876 0.999896 1.99993817 0.999927 1.999918 0.999930 1.999935 0.999945 1.99996718 0.999962 1.999957 0.999963 1.999966 0.999971 1.999983 结论分析：比较三种迭代方法，jaccobi迭代法迭代次数是29次，gauss-seidel迭代法迭代次数是16，速度要比jaccobi迭代法快。而在sor迭代中选取不同的松弛因子w迭代的次数各不相同且可能会相差很大，在w=1.95时迭代次数有242次。所以改变松弛因子w的值，不仅会影响迭代过程的收敛速度，也会影响其收敛性。5用逆幂迭代法求最接近于11的特征值和特征向量，准确到。程序：A=6 3 1;3 2 1;1 1 1;x,lamda=eig(A), v=1 1 1;q=11,B=inv(A-q*eye(3);for k=1:20 z=B*v; m=max1(z); v=z/m; fprintf(%6d %f %f %f %fn,k,v(1),v(2),v(3),m-1+q), if abs(m-1+q-lamda(3,3)0.001 break endend运行结果：x = -0.3065 -0.4082 -0.8599 0.7812 0.4082 -0.