PROGRAM EULER7 * * Double precision coding is used IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA E0, E1, E2, E3, E4, E5/1.000000000D0, 0.864652904D0, + 0.486589996D0,-0.038007436D0, + -0.550798926D0,-0.899193465D0/ * * Open a file for output OPEN(1,FILE='pend7.out') * * Problem parameters * * Pendulum length (m) XL=30.0D0 * * Acceleration of gravity (m/s**2) G=9.8D0 * * Step through a series of Euler integrations * (with varying step size, h) DO IH=1,4 IF(IH.EQ.1)H=1.0D0 IF(IH.EQ.2)H=0.1D0 IF(IH.EQ.3)H=0.01D0 IF(IH.EQ.4)H=0.001D0 * * Begin integration (set initial conditions) * * Initial time (sec) T=0.0D0 * * Initial angle (radians) * Nonlinear problem Y1=1.0D0 * * Initial velocity (radians/sec) * Nonlinear problem Y2=0.0D0 * * Final t (sec) TF=5.0D0 * * Number of steps for 0 le t le TF NSTEPS=INT(TF/H) * * Perform integration DO I=0,NSTEPS * * Write numerical solutions * * t=0 IF(T.LT.0.00001D0)THEN WRITE(*,2) WRITE(1,2) 2 FORMAT(//,9X,'h',7X,'t',13X,'y1',10X,'exact',10X,'error') EXACT=E0 * * t = 1, 2, 3, 4, 5 ELSE IF ( + ((T.GT.0.99999D0).AND.(T.LT.1.00001D0)))THEN EXACT=E1 ELSE IF ( + ((T.GT.1.99999D0).AND.(T.LT.2.00001D0)))THEN EXACT=E2 ELSE IF ( + ((T.GT.2.99999D0).AND.(T.LT.3.00001D0)))THEN EXACT=E3 ELSE IF ( + ((T.GT.3.99999D0).AND.(T.LT.4.00001D0)))THEN EXACT=E4 ELSE IF ( + ((T.GT.4.99999D0).AND.(T.LT.5.00001D0)))THEN EXACT=E5 END IF ERROR=Y1-EXACT WRITE(*,1)H,T,Y1,EXACT,ERROR WRITE(1,1)H,T,Y1,EXACT,ERROR 1 FORMAT(F10.4,F8.2,3F15.9) * * Derivative at base point * * Nonlinear problem Y1B=Y1 Y2B=Y2 DY1DTB=Y2B DY2DTB=-(G/XL)*DSIN(Y1B) * * Euler step T=DFLOAT(I+1)*H * * Nonlinear problem Y1=Y1B+DY1DTB*H Y2=Y2B+DY2DTB*H * * Derivative at advanced points * * Nonlinear problem DY1DT=Y2 DY2DT=-(G/XL)*DSIN(Y1) * * Modified Euler step * * Nonlinear problem Y1=Y1B+(DY1DTB+DY1DT)/2.0D0*H Y2=Y2B+(DY2DTB+DY2DT)/2.0D0*H * * Continue Euler integration END DO * * Integration complete; next integration step size WRITE(*,3)NSTEPS WRITE(1,3)NSTEPS 3 FORMAT(/,3X,'NSTEPS = ',I6) PAUSE END DO * * Calculation for all integration steps completed WRITE(*,4) WRITE(1,4) 4 FORMAT(///,'Conclusion: ',/,3X, + 'The modified Euler method is second-order ',/,3X, + 'correct for the nonlinear pendulum problem.') * * End of calculations STOP END DOUBLE PRECISION FUNCTION DFLOAT(I) * * Function DFLOAT converts a single precision integer to a * double precision float * DFLOAT=DBLE(REAL(I)) RETURN END