1、笫八章 常用算法程序举例,本章介绍的一些例子，这是计算机解题中所常遇到的，通过它们可以学习程序设计的方法与技巧。 切实掌握基本的算法，并在此基础上举一反三。,8.1 数值积分,求一个函数f(x)在a,b上的定积分 ，其几何意义是求 f(x)曲线和直线 x=a，y=0，x=b 所围成的曲边梯形面积。 为了近似求出此面积，可将a,b区间分成若干个小区间，每个区间的宽度为(ba)n，n为区间个数。 近似求出每个小的曲边梯形面积，然后将n个小面积加起来，就近似得到总的面积。即定积分的近似值，当n愈大(即区间分得愈小，近似程度愈高。,y f(b)f(a)a a+h ba+(i-1)h a+ih图 8.1
2、,8.1 数值积分,近似求小曲边梯形面积的方法，常用的有以下三种： (l)用小矩形代替小曲边梯形，求出各小矩形面积，然后累加 (2)用小梯形代替小曲边梯形。 (3)在小区间范围内，用一条抛物线代替该区间内的f(X)，然后求出由该抛物线与 x=a+(i-1) h， x=0， x=a+ih 形成的小曲边梯形面积。,图 8.2 图 8.3,8.1.l 矩形法,矩形的面积为底高。第一个小矩形的面积，底的值就是 (b-a)/n，高为 f(a)，当然也可以用 f(a+h) 为高。第i个小矩形的面积为：si= hf(a+(i-l)h)用NS图描述求定积分的算法，如右图,输入 A,B,N H=(B-A)/N
3、S=0 Do I=1,nX= A+(I-1)*HS=S+H*F(X) 打印面积 S图 8.4,8.1.l 矩形法,例 8.1 求 按图8.4写出程序。READ (*,*) A,B,N H=(B-A)/NS=0.0DO 10 I=1,NX=A+(I-1)*HS=S+H*EXP(X) CONTINUEWRITE (*,100)A,B,NWRITE(*,200)S FORMAT(1X,A=,F10.3,3X,B=,F10.3,3X,N=,I4) FORMAT(1X,S=,F15.8)END,8.1.2 梯形法,这也是近似地用小梯形表示小曲边梯形，见右图 。 笫 i 个小梯形面积为那么,yf(a+(i
4、-1)h)f(a+ih)0 xa+(i-1)h a+ih,8.1.2 梯形法,例 8.2 求 READ (*,*) A,B,N H=(B-A)/NS=0.0DO 10 I=1,NS=S+H*(SIN(0.0+(I-1)*H)+SIN(0.0+I*H)/2.0 CONTINUEWRITE (*,100) A,B,NWRITE (*,200) S FORMAT(1X,A=,F10.3,3X,B=,F10.3,3X,N=,I4) FORMAT(1X,S=,F15.8)END,8.1.3 辛普生(Sinpson)法,其基本方法是：在一小区间内用一抛物线 f1(x)代替原来的曲线 f(x)，见图85。抛
5、物线是如何决定的呢？取a，b的中点c，c的坐标为求出f(c)，过f(a),f(b),f(c)三点可以作一条唯一的抛物线：抛物线定积分公式其中,yf1(x) f(c) f(b)f(a) f(x)a c b x,8.1.3 辛普生(Sinpson)法,如果把a,b分成n个小区间，则其中,8.1.3 辛普生(Sinpson)法,用辛普生法求READ(*,*) A,B,NH=(B-A)/(2.*N)S=0.0FA=1.0/(1.0+A)FB=1.0/(1.0+B)X=A+HF2=0.0F4=1.0/(1.0+X)DO 10 I=1,N-1X=X+HF2=F2+1.0/(1.0+X),X=X+HF4=F
6、4+1.0/(1.0+X) 10 CONTINUE S=H/3.0*(FA+FB+4.0*F4+2.0*F2)WRITE(*,100)A,B,NWIRTE(*,200) S FORMAT(1X,A=,F8.2,2X, B=,F8.2,2X,N=,I4) FORMAT(1X,S=,F16.7)END,8.2 解一元方程 8.2.1 迭代法,用迭代法求一元方程f(X)的根，其基本方法如下：(l)将f(X)改写成求X的式子：Xg(X)形式。(2)大致估计出一个根X的范围，给X定一个初值X0，把它代入上式等号的右边，求出 X的第一次近似值X1。(3)再将X1代入g(X)得X2。这样一次一次地将求出的新
7、值又作为下一次的初值代入g(X)。 直到前后两次求出的X值很接近，即 为止。是一个给定的很小的数。这时 xn+1 就是所求的近似值。 例8.4 用迭代法求解x32x22x10先找出xg(x)式子，可得到： x=(-x3-2x2-1)/2,8.2.1 迭代法,估计在0附近有一个根。设x0=0。 在设计算法时，应考虑一个问题,有可能经过许多次迭代后仍不收敛。为防止无休止的迭代下去，应规定一个最高的循环次数，如达到此次数仍不能满足 就不再进行下去，应打印出“经xx次迭次后仍未收敛”。应当说明：对同一个f(X)，可以写出不同的X=g(X)形式。例如，对X32X22X+10，可以写出： (l) X=(-
8、X3-2X2-2)2(2) X=(-X3-2X-2)/2)1/2(3) X(-2X2-2X-2)1/3 有的g(X)是收敛的，而有的g(X)则不收敛。同一个g(X)，对某些X0来说是收敛的，对有的X0则不收敛。 如果g(x)具有一阶导数连续，且对于所有的x，若有 (q为一个定数)，那么 x=g(x) 对于任意的x0均收敛，且q愈小，收敛速度愈快。如果不满足对所有的X存在 则可能对有的x0收敛，对有的则不收敛。,8.2.1 迭代法,READ (*,*) X,MDO 10 I=1,MX1=(-X*3-2.0*X*X-2.0)/2.0WRITE(*,100) I,X1IF(ABS(X-X1).GT.
9、1E-6) THENX=X1ELSESTOPEND IF CONTINUEWRITE(*,200) M FORMAT(1X,I=,I3,5X,X1=,F15.7) FORMAT(1X,COMPUTATION HAS NOT CONVERGED AFTER,I4,ITERATION)END,8.2.1 牛顿迭代法,用牛顿选代法求 f(X)=0 在 X。附近的一个实根的方法是(l)选一个接近于真实根X的近似根X1；(2)通过X1求出 f(X1)。在几何上就是作X=X1，交f(X)于f(X1)。见下页图8.9 。(3)过f(X1)作f(X)的切线，交X轴于X2。可以用公式求出X2。(4) 过f(X2
10、)作f(X)的切线交X轴于X3。当 就认为XN+1已足够精确。牛顿迭代公式是：,8.2.1 牛顿迭代法,yf(x)f(x1)f(x2)f(x3) a 0 x3 x2 x1图 8.9,8.2.1 牛顿迭代法,例8.5 用牛顿迭代法求f(X)=X3-2X24X十1=0在X0附近的一个实根。已知：f(X)X3-2X24X+l求出： f(X)3X2-4X+4可以作出NS流程图(图810)。其思路是：X1 第一次X的初始值， N 迭代次数。F f(Xl)，FI f(X1)，X 通过迭代求出的X的下一个近似值。,读X1N=1F=X3-2X24X十1=0F13X2-4X+4 X=X1-F/F1N=N+1 直
11、到 X-X1 10-6图 8.10,8.2.1 牛顿迭代法,程序READ (*,*) X1N=1 F=X1*3-2.0*X1*2+4.0*X1+1F1=3.0*X1*2-4.0*X1+4.0X=X1-F/F1WRITE(*,100) N,X1 N=N+1IF(ABS(X-X1).GT.1E-6) thenX1=XGOTO 10END IF FORMAT(1X,N=,I3,3X,X1=,F15.7)END,8.2.3 二分法,1. 其基本思路2. 怎样做到“舍弃”一半区间呢？ (1)如果f(x)与f(x1)不同符号则用x作为新的x2，这就舍掉了原(x,x2)这个区间。否则，就舍弃了原(x1,x)
12、区间。(2)再根据新的 x1，x2，找中点x，重复上述步骤。,y f(x2)f(x)x1 x x2 xf(x1) x1 x x2图 8.11,8.2.3 二分法,3. 例8.6 用二分法求 f(X)=X3-6X-l=0 在X=2附近的一个实根。取X1=0，X2=5。READ(*,*) X1, X2 F1=X1*3-6.0*X1-1.0 F2=X2*3-6.0*X2-1.0 10 X=(X1+X2)/2.0F=X*3-6.0*X-1.0IF(SIGN(F,F1).EQ.F) THEN / if( f1.ne.0.0 .and. f/f1 .gt. 0.0) thenX1=XF1=FELSEX2=
13、XF2=FIF( (ABS(X1-X2).GT.1E-5) .AND.(ABS(F).GT.1E-6) GOTO 10IF (ABS(F).GT.1E-6) X=(X1+X2)/2.0WRITE(*,100) XFORMAT(1X,X=,F15.7)END,8.2.4 弦截法,弦截法的基本思路与“二分法相似，有一点不同，二分法每次取区间的中点，然后从中舍弃一半的区问。而弦截法则是取 f(X1)与f(X2)联线与 X轴的交点 X，从(X1,X)与(X,X2) 二个区问中舍去一个。取舍的方法与二分法相同。见图8.12。如果f(X1)与f(X)异号，则它们之间的联 y f(x2)线必然交X轴，设交点
14、为X 可用下式求出X判别 f(X1)与 f(X2)是否同符号的方法与 x1 x x2二分法采用的方法相同。 x例8.7 用弦裁法求根。X3-2X27X40 f(x1) f(x)要求进行到先后二次求出的X(弦截点)的值之差不超过10-6为止 x1 x2 图 8.12,8.2.4 弦截法,10 READ( *,*) X1, X2F1=X1*3-2.0*X1*2+7.0*X1+4.0F2=X2*3-2.0*X2*2+7.0*X2+4.0IF(SIGN(F1, F2).EQ.F1) GOTO 10F=1.0 20 IF(ABS(X1-X2).GT.1E-5).AND.(ABS(F).GT.1E-6)
15、THENX=X2-(X2X1)(F2F1)*F2F=X*3-2.0*X*27.0*X+4.0IF(SIGN(F,F1).EQ.F) THENX1=X,8.2.4 弦截法,F1=FELSEX2XF2=FEND IFGOTO 20END IFIF(ABS(F).GT. 1E6) X=(X1+X2)/2.0WRITE(*,100) XFORMAT(1X,X=，F15.7)END,8.3 求函数的极小值,0.618法求函数 F(X) 在给定区间a1，a2上的 极小值： 一个试验点a4 取在a1,a 2的 0.618位置上，令一个试验点a3 取在区间 a1，a2的 0. 382的位置上，即a3与a4 为
16、区间a1，a2上的两个对称点。a3a10.382(a2a1) a4a10.618(a2a1) 若函数F(a3) F(a4） 说明a3是好点，极 小点必在a1,a 4之内。否则在a3，a2 之内。然后在保留区间内继续取点，连 续 N1次缩短后，最后的区间长度为： 0.618(N1)(a2a1)应满足精度要求，若相对精度为， 则0.618(N1) ,F(X)F(a1) f(a2)f(a3) f(a4)xa1 a3 a4 a2a1 a3 a4 a2图 8.13,0.618法黄金分割法,二. 具体步骤： (1)确定试点个数N根据精度要求，0.618(N1) Nlnln0.618; N为整数NINT(l
17、n()ln(0.618)1;(2)取试点(令t0.618)a3a2t(a1a2) a4a1t(a2a1）计算F(a3),F(a4),并令F3F(a3), F4F(a4),0.618法黄金分割法,(3) 比较F3和F4如果 F3 F4 令a2a4; a4a3; F4F3a3a2t(a1a2）;F3F(a3)否则。令a1a3； a3a4; F3F4a4a1t(a2a1); F4F(a4);(4) 上述计算从N1开始。以1为步长，直到1为止。(5) 比较F3与F4,如果F3 F4,则a3为极小点，F3为极小值。否则a4为极小点，F4 为极小值。,0.618法黄金分割法,三.黄金分割法（0.618法）
18、的FORTRAN程序 1.子程序段头形式SUBROUTINE GOLD(X,H,EPS) 2.哑元说明X 实型量，开始时存放由使用者给定的初始点，计算完毕时，存放极小点；H 实型量，开始时存放由使用者给定的初始步长， 用以寻找搜索区间a,b,离开子程序时，存放函数最优值；EPS 实型量， 当区间长度与初始区间长度之比小于EPS 时，迭代终止。 3.子程序功能,0.618法黄金分割法,子程序GOLD首先从初始点出发，确定搜索区间，然后，用 0.618法，缩小区间，求出近似极小点。作为例子，计算了F(X)0.5(exe-x) 的极小值。取X1.0 , h0.1 , EPS10-5 ,计算结果输出为
19、Xmin ，Fmin ，及 NF其中 NF 为计算函数的次数。该例题的计算程序(见D:FORGOLD.FOR)。4.黄金分割法（0.618法）子程序（见gold.for）。,GOLD.FOR,PROGRAM JYZCOMMON NFNF=0X=1.0H=0.1CALL GOLD(X,H,1.E-5)WRITE(*,10) X,H,NFWRITE(2,10) X,H,NF 10 FORMAT(1X,XMIN=,E15.6,5X,FMIN=,E15.6,5X,NF=,I4)STOPEND FUNCTION F(X)COMMON NFF=0.5*(EXP(X)+EXP(-X)NF=NF+1,GOLD
20、.FOR,RETURN ENDSUBROUTINE GOLD(X,H,EPS)F0=F(X)IF(F(X+H)-F0) 1,1,2 2 X2=X+HH=-H 1 X1=X+HF1=F(X1)IF(F1.LE.F0) THEN 3 H=2*HX2=XX=X1F0=F1,GOLD.FOR,GO TO 1 ELSE 4 A=X2 B=X2IF(X2.GT.X1) A=X1IF(X2.LT.X1) B=X1END IFS=B-AT=(SQRT(5.0)-1.0)*0.5X1=A+T*(B-A)F1=F(X1)X2=B+T*(A-B)F2=F(X2) 55 IF(F1.GE.F2) GO TO 90,G
21、OLD.FOR,A=X2X2=X1F2=F1X1=A+T*(B-A) F1=F(X1)GO TO 115 90 B=X1X1=X2 F1=F2X2=B+T*(A-B)F2=F(X2) 115 IF(ABS(B-A)/S).GT.EPS) GO TO 55,GOLD.FOR,IF(F1.GE.F2) GO TO 135X=X1 H=F1GO TO 140 135 X=X2H=F2 140 RETURNEND,8.4 打印图案,更重要的绘图技术将在 第14章中介绍,8.5 计算机模拟,计算机模拟也称为“仿真” 例10. 有一球以10m/s的速度从水平线以450向斜上方抛出，该球着地后仍保持原方向向
22、前弹跳。但速度是抛出时的80%，以后每着地一次速度都比上次减小20%。求此球运动轨迹。 YX=X0+(V0cos450)TY=Y0+(V0sin450)T-0.5gT2 g=9.8m/s2 当Y0=0球上抛越过顶点开始下落着地，着地时刻Ti=2(V0*0.8i)sin450/g 0 X I=0 T0=2V0sin450/g=2*10*sin450/9.8=1.443075064S I=1 T1= 2V0*0.8*sin450/g =2*10*0.8*sin450/9.8=1.154460051S,程序,PROGRAM INCLUDE FGRAPH.FIINCLUDE FGRAPH.FDPARA
23、METER (N=100,M=5)DIMENSION X(N),Y(N)INTEGER*4 BCOLINTEGER*2 du,MAXX,MAXY,NEWX,NEWYCOMMON MAXX,MAXYRECORD /xycoord/xyMAXX=10MAXY=10WRITE(*,*) ENTER V0,DLREAD(*,*) V0,DLdu=setvideomode($vres16color),BCOL=setbkcolor($WHITE) du=setcolor(int2(1)X0=0.Y0=0.DO 100 I=0,15TI=2.0*V0*DL*I*SIN(3.1415926/4.0)/9.8
24、DT=TI/REAL(N)DO 200 J=1,NX(J)=X0+V0*DL*I*COS(3.1415926/4.)*DT*JY(J)=Y0+V0*DL*I*SIN(3.1415926/4.0)*DT*J-0.5*9.8*(DT*J)*2 200 CONTINUECALL setvieworg(20,240,xy)IF(I.EQ.0) CALL MOVETO(NEWX(X0),NEWY(-Y0),XY) IF(I.EQ.0) DU=LINETO(580,NEWY(-Y0) IF(I.EQ.0) CALL MOVETO(NEWX(X0),NEWY(-Y0),XY),DO 300 J=1,N 30
25、0 du=lineto(NEWX(X(J),NEWY(-Y(J)X0=X(N)Y0=Y(N) 100 CONTINUEREAD(*,*)dummy=setvideomode($DEFAULTMODE)STOP ENDINTEGER*2 FUNCTION NEWX(XCOORD)INTEGER*2 MAXX,MAXYREAL*4 TEMPXCOMMON MAXX,MAXY,TEMPX=200.0/MAXXTEMPX=XCOORD*TEMPXNEWX=TEMPXENDINTEGER*2 FUNCTION NEWY(YCOORD)INTEGER*2 MAXX,MAXYREAL*4 TEMPYCOMMON MAXX,MAXYTEMPY=200.0/MAXX TEMPY=YCOORD*TEMPYNEWY=TEMPYEND,程序运行结果,程序运行结果,
