PROGRAM PRO2P13 C... C... PROGRAM PRO2P13 CALLS: (1) SUBROUTINE INITAL TO DEFINE THE ODE C... INITAL CONDITIONS, (2) SUBROUTINE RKF45 TO INTEGRATE THE ODES, C... AND (3) SUBROUTINE PRINT TO PRINT THE SOLUTION. C... C... THE FOLLOWING CODING IS FOR 250 ODES. IF MORE ODES ARE TO BE INTE- C... GRATED, ALL OF THE 250'S SHOULD BE CHANGED TO THE REQUIRED NUMBER IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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(1600), 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 C... C... DEFINE THE INPUT/OUTPUT UNIT NUMBERS NI=5 NO=6 C... C... OPEN INPUT AND OUTPUT FILES OPEN(NI,FILE='DATA' ,STATUS='OLD') OPEN(NO,FILE='OUTPUT',STATUS='NEW') C... C... TOTAL NUMBER OF RUNS NORUNS=22 C... C... STEP THROUGH NORUNS RUNS DO 1 NORUN=1,NORUNS C... C... INITIALIZE THE RUN TERMINATION VARIABLE NSTOP=0 C... C... READ THE FIRST LINE OF DATA IF(NORUN.EQ.1)READ(NI,1000,END=999)(TITLE(I),I=1,20) C... C... READ THE SECOND LINE OF DATA IF(NORUN.EQ.1)READ(NI,1001,END=999)T0,TF,TP C... C... READ THE THIRD LINE OF DATA IF(NORUN.EQ.1)READ(NI,1002,END=999)NEQN,ERROR C... C... PRINT A DATA SUMMARY WRITE(NO,1003)NORUN,(TITLE(I),I=1,20), 1 T0,TF,TP, 2 NEQN,ERROR C... C... INITIALIZE TIME T=T0 C... C... SET THE INITIAL CONDITIONS CALL INITAL C... C... SET THE INITIAL DERIVATIVES (FOR POSSIBLE PRINTING) CALL DERV C... C... PRINT THE INITIAL CONDITIONS CALL PRINT(NI,NO) C... C... SET THE INITIAL CONDITIONS FOR SUBROUTINE RKF45 TV=T0 DO 5 I=1,NEQN YV(I)=Y(I) 5 CONTINUE C... C... SET THE PARAMETERS FOR SUBROUTINE RKF45 RELERR=ERROR ABSERR=ERROR IFLAG=1 TOUT=T0+TP C... C... CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL C... CONDITION (IFLAG = 1) OR COMPUTE THE SOLUTION TO THE NEXT PRINT C... POINT (IFLAG = 2) 4 CALL RKF45(FCN,NEQN,YV,TV,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK) C... C... PRINT THE SOLUTION AT THE NEXT PRINT POINT T=TV DO 6 I=1,NEQN Y(I)=YV(I) 6 CONTINUE CALL DERV CALL PRINT(NI,NO) C... C... TEST FOR AN ERROR CONDITION IF(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 TOUT=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 1 CONTINUE C... C... ALL RUNS ARE COMPLETE 999 STOP C... C... ***************************************************************** C... C... FORMATS C... 1000 FORMAT(20A4) 1001 FORMAT(3E10.0) 1002 FORMAT(I5,20X,E10.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 ' MAXIMUM INTEGRATION ERROR - ',D10.3,//, 7 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... C... NOTE THAT THE SIZE OF ARRAYS Y AND F IN THE FOLLOWING COMMON AREA C... IS ACTUALLY SET BY THE CORRESPONDING COMMON STATEMENT IN MAIN C... PROGRAM STEEP IMPLICIT DOUBLE PRECISION (A-H,O-Z) 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 VARIABLE, 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