PROGRAM ORBIT2 C... C... FEYNMAN ORBIT PROBLEM - SOLUTION BY SUBROUTINE RKF45 C... C... SEE SUBROUTINE INITAL FOR A DISCUSSION OF THE PROBLEM. C... C... THIS MAIN PROGRAM ESSENTIALLY CALLS SUBROUTINE RKF45 FOR THE C... SOLUTION OF A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS (ODES). C... C... THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND C... THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV. THE C... NUMERICAL SOLUTION IS PRINTED AND PLOTTED IN SUBROUTINE PRINT. C... C... DOUBLE PRECISION USED TO PERMIT SMALLER ERROR TOLERANCES WHEN C... CALLING RKF45 IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... SUBROUTINES INITAL, DERV AND PRINT. THE FOLLOWING CODING IS FOR C... 250 ORDINARY DIFFERENTIAL EQUATIONS (ODES). IF MORE ODES ARE TO C... BE INTEGRATED, ALL OF THE 250*S SHOULD BE CHANGED TO THE REQUIRED C... NUMBER COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(250) 2 /F/ F(250) C... C... THE NUMBER OF DIFFERENTIAL EQUATIONS IS IN COMMON/N/ FOR USE IN C... SUBROUTINE FCN COMMON/N/ NEQN C... C... COMMON AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER C... SUBROUTINES COMMON/IO/ NI, NO C... C... ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY RKF45 DIMENSION YV(250), WORK(1503), IWORK(5) C... C... EXTERNAL THE DERIVATIVE ROUTINE CALLED BY RKF45 EXTERNAL FCN C... C... ARRAY FOR THE TITLE (FIRST LINE OF DATA), CHARACTERS END OF RUNS CHARACTER TITLE(20)*4, ENDRUN(3)*4 C... C... DEFINE THE CHARACTERS END OF RUNS DATA ENDRUN/'END ','OF R','UNS '/ C... C... DEFINE THE INPUT/OUTPUT UNIT NUMBERS NI=5 NO=6 C... C... DEFINE INPUT/OUTPUT FILES OPEN(NI,FILE='DATA' ,STATUS='OLD') OPEN(NO,FILE='OUTPUT',STATUS='NEW') C... C... INITIALIZE THE RUN COUNTER NORUN=0 C... C... BEGIN A RUN 1 NORUN=NORUN+1 C... C... INITIALIZE THE RUN TERMINATION VARIABLE NSTOP=0 C... C... READ THE FIRST LINE OF DATA READ(NI,1000,END=999)(TITLE(I),I=1,20) C... C... TEST FOR END OF RUNS IN THE DATA DO 2 I=1,3 IF(TITLE(I).NE.ENDRUN(I))GO TO 3 2 CONTINUE C... C... AN END OF RUNS HAS BEEN READ, SO TERMINATE EXECUTION 999 STOP C... C... READ THE SECOND LINE OF DATA 3 READ(NI,1001,END=999)T0,TF,TP C... C... READ THE THIRD LINE OF DATA READ(NI,1002,END=999)NEQN,TOL C... C... PRINT A DATA SUMMARY WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),T0,TF,TP,NEQN,TOL C... C... INITIALIZE TIME T=T0 C... C... SET THE INITIAL CONDITIONS CALL INITAL C... C... SET THE INITIAL CONDITIONS FOR SUBROUTINE RKF45 DO 5 I=1,NEQN YV(I)=Y(I) 5 CONTINUE C... C... SET THE PARAMETERS FOR SUBROUTINE RKF45 C... C... INITIAL TIME TV=T0 C... C... ERROR TOLERANCE ABSERR=0.D0 RELERR=TOL C... C... INITIALIZATION FOR MULTI-STEP MODE (RKF45 NOT CONSTRAINED C... BY OUTPUT INTERVAL TP) IFLAG=1 C... C... INITIALIZATION FOR ONE-STEP MODE (RKF45 CONSTRAINED BY C... OUTPUT INTERVAL TP) C... IFLAG=-1 C... C... FINAL TIME TEND=T0 C... C... CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL C... CONDITION (IFLAG = 1 OR -1) OR COMPUTE THE SOLUTION TO THE NEXT C... PRINT POINT (IFLAG = 2 OR -2) 4 CALL RKF45(FCN,NEQN,YV,TV,TEND,RELERR,ABSERR,IFLAG,WORK,IWORK) C... C... PRINT THE INITIAL CONDITION (IFLAG = 1 OR -1) OR THE SOLUTION AT C... THE NEXT PRINT POINT (IFLAG = 2 OR -2) T=TV DO 6 I=1,NEQN Y(I)=YV(I) 6 CONTINUE CALL PRINT(NI,NO) C... C... TEST FOR AN ERROR CONDITION IF((IFLAG.NE.2).AND.(IFLAG.NE.-2))THEN C... C... PRINT A MESSAGE INDICATING AN ERROR CONDITION WRITE(NO,1004)IFLAG C... C... GO ON TO THE NEXT RUN GO TO 1 END IF C... C... CHECK FOR A RUN TERMINATION IF(NSTOP.NE.0)GO TO 1 C... C... CHECK FOR THE END OF THE RUN TEND=TV+TP IF(TV.LT.(TF-0.5D0*TP))GO TO 4 C... C... THE CURRENT RUN IS COMPLETE, SO GO ON TO THE NEXT RUN GO TO 1 C... C... ***************************************************************** C... C... FORMATS C... 1000 FORMAT(20A4) 1001 FORMAT(3D10.0) 1002 FORMAT(I5,20X,D10.0) 1003 FORMAT(1H1, 1 ' RUN NO. - ',I3,2X,20A4,//, 2 ' INITIAL T - ',D10.3,//, 3 ' FINAL T - ',D10.3,//, 4 ' PRINT T - ',D10.3,//, 5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//, 6 ' INTEGRATION ALGORITHM - RKF45 ',//, 7 ' MAXIMUM INTEGRATION ERROR - ',D10.3,//, 8 1H1) 1004 FORMAT(1H ,//,' IFLAG = ',I3,//, 1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN' ,/, 2 ' IS TERMINATED. PLEASE REFER TO THE DOCUMENTATION FOR' ,/, 3 ' SUBROUTINE',//,25X,'RKF45',//, 4 ' FOR AN EXPLANATION OF THESE ERROR INDICATORS' ) END SUBROUTINE FCN(TV,YV,YDOT) C... C... SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES RKF45 C... AND DERV C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... ODE COMMON AREA COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(1) 2 /F/ F(1) C... C... THE NUMBER OF DIFFERENTIAL EQUATIONS IS AVAILABLE THROUGH COMMON C... /N/ COMMON/N/ NEQN C... C... ABSOLUTE DIMENSION THE DEPENDENT, DERIVATIVE VECTORS DIMENSION YV(250), YDOT(250) C... C... TRANSFER THE INDEPENDENT VARIABLE, DEPENDENT VARIABLE VECTOR C... FOR USE IN SUBROUTINE DERV T=TV DO 1 I=1,NEQN Y(I)=YV(I) 1 CONTINUE C... C... EVALUATE THE DERIVATIVE VECTOR CALL DERV C... C... TRANSFER THE DERIVATIVE VECTOR FOR USE BY SUBROUTINE RKF45 DO 2 I=1,NEQN YDOT(I)=F(I) 2 CONTINUE RETURN END SUBROUTINE INITAL C... C... THE FOLLOWING SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS C... (ODES) FOR PLANETARY MOTION IS DISCUSSED BY FEYNMAN ET AL C... (1963) C... C... 2 2 C... D X/DT = - X/R**3 (1) C... C... 2 2 C... D Y/DT = - Y/R**3 (2) C... C... R = (X**2 + Y**2)**0.5 (3) C... C... X(0) = 0.5, DX(0)/DT = 0 (4)(5) C... C... Y(0) = 0, DY(0)/DT = 1.63 (6)(7) C... C... THE TOTAL ENERGY FOR THE SYSTEM, WHICH IS CONSTANT, IS GIVEN C... BY C... C... E = (1/2)*(V **2 + V **2) - 1/R (8) C... X Y C... WHERE C... C... V = DX/DT, V = DY/DT C... X Y C... C... E FROM EQUATION (8) IS COMPUTED AND PRINTED IN SUBROUTINE C... PRINT TO TEST THE OPERATION OF THE INTEGRATOR. C... C... EQUATIONS (1) AND (2) CAN BE STATED IN TERMS OF NEW VARIABLES C... DEFINED AS C... C... Y1 = X, Y2 = DX/DT C... C... Y3 = Y, Y4 = DY/DT C... C... THE FOLLOWING PROGRAMMING IS IN TERMS OF X1, X2, Y1 AND Y2. C... C... EQUATIONS (1) TO (7) HAVE AN EXACT SOLUTION (GREENSPAN, 1990), C... WHICH, IF WRITTEN IN PARAMETRIC FORM (STENGLE, 1991), IS C... C... X(T) = -C + A*COS(THETA(T)) (9) C... C... Y(T) = B*SIN(THETA(T)) (10) C... C... WHERE THETA(T) IS GIVEN BY THE ODE (2) C... C... DTHETA/DT = SQRT((1/TERM1)*TERM2) (11) C... C... TERM1 = (1/2)*((A*SIN(THETA))**2 + (B*COS(THETA))**2) (12) C... C... TERM2 = -1/(2*A) + 1/SQRT(TERM3) (13) C... C... TERM3 = SQRT((A*COS(THETA) - C)**2 + (B*SIN(THETA))**2) (14) C... C... THE INITIAL CONDITION FOR EQUATION (11) IS C... C... THETA(0) = 0 (15) C... C... THE SOLUTION TO EQUATIONS (1), (2) AND (11) IS BY ODE SOLVER C... RKF45 DISCUSSED IN FORSYTHE (1976). C... C... REFERENCES: C... C... FEYNMAN, RICHARD P., ROBERT B. LEIGHTON AND MATTHEW SANDS, LEC- C... TURES ON PHYSICS, ADDISON-WESLEY PUBLISHING COMPANY, READING, C... MA, PP 9-6 TO 9-9 (1963) C... C... FORSYTHE, GEORGE E., MICHAEL A. MALCOLM, AND CLEVE B. MOLER, C... COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-HALL, C... TATIONS, PRENTICE-HALL, ENGLEWOOD CLIFF, NJ, PP 129-147 (1976) C... C... GREENSPAN, DONALD, A COUNTEREXAMPLE OF THE USE OF ENERGY AS A C... MEASURE OF COMPUTATIONAL EFFICIENCY, J. COMPUTATIONAL PHYSICS, C... V. 91, 490-494 (1990) C... C... STENGLE, GILBERT, PRIVATE COMMUNICATION, OCTOBER 8, 1991 C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N=5) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(N) 2 /F/ DYDT(N) 3 /C/ R, A, B, C 4 /I/ IP, NP C... C... MODEL PARAMETERS A=10000.0D0/13431.0D0 B=DSQRT(664225.0D0/1343100.0D0) C=6569.0D0/26862.0D0 C... C... INITIAL CONDITIONS, EQUATIONS (4) TO (7), (11) Y(1)=0.5D0 Y(2)=0.D0 Y(3)=0.D0 Y(4)=DSQRT(2.0D0*(-1.0D0/(2.0D0*A)+1.0D0/(A-C))) Y(5)=0.D0 C... C... COMPUTE THE INITIAL DERIVATIVES CALL DERV C... C... INITIALIZE INTEGER VARIABLES TO CONTROL THE PRINTING IP=0 NP=21 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N=5) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(N) 2 /F/ DYDT(N) 3 /C/ R, A, B, C 4 /I/ IP, NP C... C... RADIUS (WITH CUBING TO MINIMIZE CALCULATIONS) R=DSQRT(Y(1)**2+Y(3)**2)**3 C... C... EQUATION (1) DYDT(1)=Y(2) DYDT(2)=-Y(1)/R C... C... EQUATION (2) DYDT(3)=Y(4) DYDT(4)=-Y(3)/R C... C... EQUATION (11) SINE=DSIN(Y(5)) COSE=DCOS(Y(5)) TERM1=0.5D0*((A*SINE)**2+(B*COSE)**2) TERM2=-1.0D0/(2.0D0*A)+1.0D0/DSQRT((A*COSE-C)**2+(B*SINE)**2) DYDT(5)=DSQRT((1.0D0/TERM1)*TERM2) RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N=5) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(N) 2 /F/ DYDT(N) 3 /C/ R, A, B, C 4 /I/ IP, NP DIMENSION YE(N), ER(N) C... C... PRINT THE NUMERICAL SOLUTION C... C... ENERGY, EQUATION (8) EN=ENERGY(N,Y) C... C... EXACT SOLUTION AND ERROR IN NUMERICAL SOLUTION C... (THETA = Y(5), DTHETA/DT = DYDT(5) IN EQUATION (11)) THETA= Y(5) DTHDT=DYDT(5) YE(1)=-C+A*DCOS(THETA) YE(2)=-A*DSIN(THETA)*DTHDT YE(3)= B*DSIN(THETA) YE(4)= B*DCOS(THETA)*DTHDT DO 2 I=1,N-1 ER(I)=Y(I)-YE(I) 2 CONTINUE C... C... PRINT NUMERICAL AND EXACT SOLUTIONS, NUMERICAL ERROR AND C... ENERGY WRITE(NO,1)T,THETA,(I,Y(I),YE(I),ER(I),I=1,N-1),EN 1 FORMAT(' T = ',F6.3,' THETA = ',F6.3,//, + 4(2X, ' I = ',I2, ' Y = ',F9.6, + ' YE = ',F9.6, ' ER = ',E9.3,/),/,2X, + ' EN = ',F9.6,///) C... C... STORE THE NUMERICAL SOLUTION FOR PLOTTING IP=IP+1 CALL PLOT RETURN END DOUBLE PRECISION FUNCTION ENERGY(N,Y) C... C... FUNCTION ENERGY COMPUTES THE TOTAL ENERGY OF THE SYSTEM C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N) C... C... ENERGY, EQUATION (8) R=DSQRT(Y(1)**2+Y(3)**2) ENERGY=0.5D0*(Y(2)**2+Y(4)**2)-1.0D0/R END SUBROUTINE PLOT IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (N=5) COMMON/T/ T, NSTOP, NORUN 1 /Y/ Y(N) 2 /F/ DYDT(N) 3 /C/ R, A, B, C 4 /I/ IP, NP C... C... INITIATE PLOTTING IF(IP.EQ.1)THEN C... C... OPEN FILE FOR TOP DRAWER PLOTTING IT=4 OPEN(IT,FILE= 'ORBIT.TOP',STATUS='UNKNOWN') C... C... SCALE THE AXES WRITE(IT,7) 7 FORMAT(' SET LIMITS X FROM -1 TO 1 Y FROM -1 TO 1',/, 1 ' SET FONT DUPLEX') WRITE(IT,4) 4 FORMAT(' SET WINDOW X 2 TO 8 Y 2 TO 8') END IF C... C... WRITE Y(3) VS Y(1) FOR PLOTTING WRITE(IT,5)Y(1),Y(3) 5 FORMAT(2F10.4) C... C... COMPLETE PLOTTING IF(IP.EQ.NP)THEN C... C... CONNECT POINTS WRITE(IT,6) 6 FORMAT(' JOIN 1') C... C... LABEL THE AXES WRITE(IT,8)NP 8 FORMAT( 1 ' TITLE 2.75 8.5 "',I3,'-point orbit from RKF45 - DP"' 2 ,/,' TITLE LEFT " y3(t)"' 3 ,/,' TITLE BOTTOM "y1(t)"' 4 ,/,' TITLE 2.5 0.7 "Greenspan/Feynman Orbit Problem"') C... C... END PLOT WRITE(IT,10) 10 FORMAT(' NEW FRAME') END IF RETURN END