PROGRAM EULER1 * * Double precision coding is used IMPLICIT DOUBLE PRECISION (A-H,O-Z) * * Open a file for output OPEN(1,FILE='pend1.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) * Linear problem Y3=1.0D0 Y30=Y3 * * Initial velocity (radians/sec) * Linear problem Y4=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(//,7X,'h',7X,'t',8X,'y3',7X,'y3e',6X,'error') ERROR=0.0D0 WRITE(*,1)H,T,Y3,Y30,ERROR WRITE(1,1)H,T,Y3,Y30,ERROR 1 FORMAT(F8.4,F8.2,2F10.4,F11.6) * * t = 1, 2, 3, 4, 5 ELSE IF ( + ((T.GT.0.99999D0).AND.(T.LT.1.00001D0)).OR. + ((T.GT.1.99999D0).AND.(T.LT.2.00001D0)).OR. + ((T.GT.2.99999D0).AND.(T.LT.3.00001D0)).OR. + ((T.GT.3.99999D0).AND.(T.LT.4.00001D0)).OR. + ((T.GT.4.99999D0).AND.(T.LT.5.00001D0)))THEN * * Exact solution and error Y3E=Y30*DCOS(DSQRT(G/XL)*T) ERROR=Y3-Y3E * * Write numerical and exact solutions, error WRITE(*,1)H,T,Y3,Y3E,ERROR WRITE(1,1)H,T,Y3,Y3E,ERROR END IF * * Derivative at base point * Linear problem DY3DT=Y4 DY4DT=-(G/XL)*Y3 * * Euler step T=DFLOAT(I+1)*H * * Linear problem Y3=Y3+DY3DT*H Y4=Y4+DY4DT*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, + 'Euler''s method is first order correct,',/,3X, + 'i.e., O(h)') * * 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