C* *MAIN 10
C*---------- THIS IS A MAIN PROGRAM FOR STINT-TYPE SUBROUTINES -------*MAIN 20
C*--------------------- SET HERE FOR N = 4 PROBLEM --------------------*MAIN 30
C* *MAIN 40
DOUBLE PRECISION D, EPS, ERO, ERORI, HI, HMAX, HMIN, HNEXT, MAIN 50
1 H0, RJ, RW, S, SAVE, T, TF, TI, TOUT, MAIN 60
2 TOUTP, TS, YDOT, YI, YMAX, YOUT, Y0 MAIN 70
C DOUBLE PRECISION DABS, DBLE, DMAX1 MAIN 80
DIMENSION YI(4), ERORI(4), Y0(4,8,4), YDOT(4,4), SAVE(52), RJ(16),MAIN 90
1 YMAX(4), RW(16), ERO(4), YOUT(4), IPIV(4) MAIN 100
COMMON /STAT/ NSTEP,NFE,NJE,NINVS MAIN 110
DATA LIN/5/,LOUT/6/ MAIN 120
C* *MAIN 130
C* SET THE INPUT PARAMETERS *MAIN 140
C* *MAIN 150
C* NP PROBLEM IDENTIFIER *MAIN 160
C* NC NUMBER OF TIMES THE SAME PROBLEM IS SOLVED WITH DIFFERENT*MAIN 170
C* VALUE OF EPS *MAIN 180
C* IP NUMBER OF TEST POINTS AT WHICH THE SOLUTION IS SOUGHT *MAIN 190
C* (DOES NOT INCLUDE THE FINAL POINT) *MAIN 200
C* N NUMBER OF DIFFERENTIAL EQUATIONS TO BE SOLVED *MAIN 210
C* YI THE INITIAL VALUES OF THE DEPENDENT VARIABLE *MAIN 220
C* TI THE INITIAL VALUE OF THE INDEPENDENT VARIABLE *MAIN 230
C* TF THE FINAL VALUE OF THE INDEPENDENT VARIABLE *MAIN 240
C* HI THE INITIAL STEP-SIZE *MAIN 250
C* MF METHOD FLAG DETERMINING THE MODE OF THE JACOBIAN *MAIN 260
C* MATRIX EVALUATION *MAIN 270
C* EPS USERS SPECIFIED ERROR PER STEP *MAIN 280
C* *MAIN 290
READ(LIN,10) NP MAIN 300
READ(LIN,10) NC MAIN 310
READ(LIN,10) IP MAIN 320
READ(LIN,10) N MAIN 330
10 FORMAT(I2) MAIN 340
READ(LIN,20) (YI(I),I=1,N) MAIN 350
READ(LIN,20) TI, TF, HI MAIN 360
DO 1000 K=1,NC MAIN 370
READ(LIN,10) MF MAIN 380
READ(LIN,20) EPS MAIN 390
20 FORMAT(3D22.15) MAIN 400
WRITE(LOUT,30) NP, MF, EPS MAIN 410
30 FORMAT(1H1,10X,13H PROBLEM NO. ,I2// MAIN 420
1 11X,16H METHOD FLAG = ,I2// MAIN 430
2 11X,18H ERROR PER STEP = ,D10.3// MAIN 440
3 5X,4HTIME,15X,29HS O L U T I O N / E R R O R,26X, MAIN 450
4 7H H-USED,3X,7HNQ-USED,3X,6H STEPS,3X,5H FNS ,3X, MAIN 460
5 5H JACS//) MAIN 470
DO 40 I=1,N MAIN 480
C* *MAIN 490
C* SET THE NORMALIZING VECTOR YMAX *MAIN 500
C* *MAIN 510
YMAX(I)=DMAX1(DABS(YI(I)),1.D0) MAIN 520
C* *MAIN 530
ERORI(I)=0.D0 MAIN 540
40 Y0(I,1,1)=YI(I) MAIN 550
C* *MAIN 560
C* THE FOLLOWING PARAMETERS ARE USED FOR STATISTICAL PURPOSES *MAIN 570
C* *MAIN 580
C* NSTEP NUMBER OF ACCEPTED STEPS *MAIN 590
C* NFE NUMBER OF FUNCTION EVALUATIONS *MAIN 600
C* NJE NUMBER OF JACOBIAN EVALUATIONS *MAIN 610
C* NINVS NUMBER OF CONVERGENCE FACTOR INVERSIONS *MAIN 620
C* *MAIN 630
NSTEP=0 MAIN 640
NFE=0 MAIN 650
NJE=0 MAIN 660
NINVS=0
C* *
C* THE FIRST POINT AT WHICH THE SOLUTION IS SOUGHT IS SET BELOW. *
C* FOR TESTING PURPOSES THE SET OF POINTS TOUT AT WHICH THE *
C* SOLUTION IS SOUGHT IS DESIGNED SO THAT ITS DENSITY IS LARGE *
C* AT THE BEGINNING OF THE INTEGRATION INTERVAL AND THEN *
C* GEOMETRICALLY DECREASES. TOUTP IS THE LAST INTERPOLATION POINT*
C* HMAX, THE MAXIMUM STEP-SIZE, IS EQUAL TO THE DISTANCE *
C* BETWEEN THE NEXT AND LAST INTERPOLATION POINTS, MULTIPLIED *
C* BY 10. *
C* *
TOUT=TI+(TF-TI)/2.D0**IP
TOUTP=TI
HMAX=(TOUT-TOUTP)*1.D1
C* *
T =TI
HNEXT=HI
HMIN=HI*1.D-2
MAXDER=7
JSTART=0
50 CALL STINT (N, T, Y0, YDOT, SAVE, H0, HNEXT, HMIN, HMAX, EPS,
1 YMAX, KFLAG, KNEXT, JSTART, MAXDER, RW, RJ, MF, IPIV)
IF(KFLAG.GE.0) GO TO 200
WRITE(LOUT,60) KFLAG
60 FORMAT(1H0/20H COMPLETION CODE IS ,I4)
GO TO 800
200 KFLGP1=KFLAG+1
C* *
C* CHECK WHETHER THE COMPUTED SOLUTION REACHED BEYOND THE *
C* INTERPOLATION POINT TOUT *
C* *
II=0
DO 500 I=1,KFLGP1
II=II+1
TS=TOUT-T+H0*DBLE(FLOAT(I-1))
IF(TS.LT.0.D0) GO TO 500
IF(I-1) 50, 50, 510
500 CONTINUE
C* *
C* THE SOLUTION REACHED BEYOND TOUT. *
C* PERFORM INTERPOLATION AT TOUT. *
C* *
510 IND=KFLAG+3-II
IF(II.EQ.2) IND=1
S=(TS-H0)/H0
DO 600 I=1,N
D=1.D0
YOUT(I)=Y0(I,1,IND)
DO 600 J=1,JSTART
D=D*(DBLE(FLOAT(J-1))+S)/DBLE(FLOAT(J))
600 YOUT(I)=YOUT(I)+D*Y0(I,J+1,IND)
CALL ERROR(TOUT,YOUT,N,ERO)
C* *
C* DETERMINE THE MAXIMUM ERROR OF THE INTERPOLATED SOLUTION *
C* COMMITTED SO FAR AND PRINT IT OUT WITH THE SOLUTION AT TOUT *
C* *
DO 610 I=1,N
610 ERORI(I)=DMAX1(ERORI(I),DABS(ERO(I)))
WRITE(LOUT,620) TOUT, (YOUT(I),I=1,N), H0, JSTART, NSTEP, NFE,
1 NJE, (ERO(I),I=1,N)
620 FORMAT( 4H T =,D11.4,4D15.6,D12.3,3X,I4,6X,I4,4X,I4,4X,I4/15X,
1 4D15.6/)
C* *
C* SET THE NEW INTERPOLATION POINT AND CHECK WHETHER AN *
C* INTERPOLATION CAN BE PERFORMED. SET NEW HMAX *
C* *
TOUTP=TOUT
TOUT=TI+(TOUT-TI)*2.D0
IF(TOUT.GT.TF) GO TO 800
IF(TOUT.LE.T) GO TO 200
HMAX=(TOUT-TOUTP)*1.D1
GO TO 50
C* *
C* REGARDLESS OF WHETHER OR NOT THE PROBLEM HAS BEEN SOLVED, *
C* PRINT OUT THE FINAL STATISTICS *
C* *
800 WRITE(LOUT,810)
810 FORMAT(1H0/20X,20HINTERPOLATION ERRORS)
WRITE(LOUT,820) (ERORI(I),I=1,N)
820 FORMAT(1H0/10X,23H MAXIMUM ABSOLUTE ERROR//11X,4D22.8)
DO 830 I=1,N
830 ERORI(I)=ERORI(I)/YMAX(I)
WRITE(LOUT,840) (ERORI(I),I=1,N)
840 FORMAT(1H0/10X,23H MAXIMUM RELATIVE ERROR//11X,4D22.8)
WRITE(LOUT,850) NSTEP, NFE, NJE, NINVS
850 FORMAT(1H0/22H PROBLEM COMPLETED IN ,I5,6H STEPS/
1 22X,I5,21H FUNCTION EVALUATIONS/
2 22X,I5,21H JACOBIAN EVALUATIONS/
3 22X,I5,18H LU DECOMPOSITIONS)
1000 CONTINUE
STOP
END
SUBROUTINE STINT (N, T, Y, DY, SAVE, H, HNEXT, HMIN, HMAX, EPS, STNT 10
+ YMAX, KFLAG, KNEXT, JSTART, MAXORD, RW, RJ, MF,
+ IPIV)
C
C***********************************************************************
C* *
C* THIS PROGRAM INTEGRATES A SET OF N FIRST ORDER ORDINARY *
C* DIFFERENTIAL EQUATIONS. A BLOCK OF THREE OR FOUR SOLUTION POINTS,*
C* EACH SEPARATED BY A STEP-SIZE H, IS COMPUTED AT EACH CALL. THE *
C* STEP-SIZE MAGNITUDE MAY BE SPECIFIED BY THE USER AT EACH CALL. *
C* ALTERNATIVELY, IT MAY BE INCREASED OR DECREASED BY STINT WITHIN *
C* THE RANGE DABS(HMIN) TO DABS(HMAX), *
C* IN ORDER TO ACHIEVE AS LARGE A STEP AS POSSIBLE, WHILE NOT *
C* COMMITING A SINGLE STEP ERROR WHICH IS LARGER THAN EPS IN THE *
C* RMS NORM, WHERE EACH COMPONENT OF THE ERROR VECTOR IS DIVIDED BY *
C* THE CORRESPONDING COMPONENT OF YMAX. *
C* *
C* THE PROGRAM REQUIRES 4 SUBROUTINES NAMED *
C* *
C* DIFFUN(N, T, Y, DY) *
C* PEDERV(N, T, Y, RJ, N0) *
C* DEC(N, N0, A, IPIV, IER) *
C* SOL(N, N0, A, B, IPIV) *
C* *
C* DIFFUN EVALUATES THE DERIVATIVES OF THE DEPENDENT VARIABLE Y(I) *
C* FOR I=1,...,N AND T AND STORES THE RESULT IN THE ARRAY DY. *
C* *
C* PEDERV COMPUTES THE PARTIAL DERIVATIVES OF THE DIFFERENTIAL *
C* EQUATIONS AT THE VALUES Y(I) FOR I=1,...,N AND T AND STORES THE *
C* RESULT IN ARRAY RJ. THUS, RJ(J+(K-1)*N0) IS THE PARTIAL OF DY(J) *
C* WITH RESPECT TO Y(K) FOR J,K=1,...,N EVALUATED AT Y(I) FOR I=1,.. *
C* ..,N AND T. NOTE--N0 IS THE VALUE OF N ON THE FIRST CALL. IF THE*
C* ANALYTIC EXPRESSIONS FOR THE PARTIAL DERIVATIVES ARE NOT *
C* AVAILABLE, THEIR APPROXIMATE VALUES CAN BE OBTAINED BY NUMERICAL *
C* DIFFERENCING. (SEE PARAMETER MF.) IN THIS CASE SUBROUTINE PEDERV*
C* MAY SIMPLY BE *
C* *
C* SUBROUTINE PEDERV(N, T, Y, RJ, N0) *
C* RETURN *
C* END *
C* *
C* DEC PERFORMS AN LU DECOMPOSITION OF A MATRIX A. IF THE *
C* DECOMPOSITION IS SUCCESSFUL IER SHOULD BE SET TO 0, OTHERWISE IT *
C* SHOULD BE SET TO +1. *
C* *
C* SOL SOLVES THE LINEAR ALGEBRAIC SYSTEM A*X=B, FOR WHICH THE *
C* MATRIX A WAS PROCESSED BY DEC. *
C* *
C* THIS PROGRAM USES DOUBLE PRECISION FOR ALL FLOATING POINT *
C* VARIABLES, EXCEPT FOR THOSE STARTING WITH P, WHICH ARE IN SINGLE *
C* PRECISION. *
C* *
C* TEMPORARY STORAGE SPACE IS PROVIDED (BY THE USER) IN THE *
C* ARRAYS IPIV, RJ, RW, AND SAVE. THE ARRAY IPIV HOLDS A VECTOR OF *
C* THE SAME NAME. THE ARRAYS RJ AND RW ARE USED TO HOLD MATRICES OF *
C* THE SAME NAMES. THE ARRAY SAVE IS PARTITIONED AS FOLLOWS *
C* *
C* SAVE(J,I) 1.LE.J.LE.8 AND 1.LE.I.LE.N IS USED TO SAVE *
C* Y(I,J) IN CASE A STEP (AND HENCE THE WHOLE CYCLE) *
C* HAS TO BE REPEATED. *
C* SAVE(9,I) 1.LE.I.LE.N IS USED TO STORE THE DERIVATIVE OF THE *
C* I-TH DEPENDENT VARIABLE SCALED BY H. *
C* SAVE(N9+I,1) IS USED TO STORE THE DERIVATIVES AS THEY ARE *
C* COMPUTED BY DIFFUN FOR THE CORRECTOR. IT IS ALSO *
C* ACCESSED AS A COMPLETE ARRAY SAVE(N9P1,1). *
C* IN ADDITION IT IS USED IN THE ERROR CONTROL TEST. *
C* (N9 = N0 * 9 AND N9P1 = N9 + 1.) *
C* SAVE(N10+I,1) IS USED TO HOLD THE CORRECTION TERMS FOR THE ENTIRE*
C* CORRECTOR ITERATION IN THE CASE, THE CORRECTOR HAS *
C* TO BE REPEATED. (N10 = N0 * 10.) *
C* SAVE(N11+I,1) IS USED TO HOLD THE DERIVATIVES EVALUATED AT *
C* Y(I)+D AND T, WHERE D IS THE INCREMENT USED IN THE *
C* NUMERICAL DIFFERENCING SCHEME INVOKED, IN ORDER TO *
C* OBTAIN APPROXIMATE VALUES OF THE PARTIAL *
C* DERIVATIVES. (N11 = N0 * 11.) *
C* SAVE(N12+I,1) IS USED TO HOLD THE DERIVATIVES EVALUATED AT Y(I) *
C* AND T IN ORDER TO OBTAIN APPROXIMATE VALUES OF THE *
C* PARTIAL DERIVATIVES. (N12 = N0 * 12.) *
C* *
C* THE PARAMETERS IN THE CALLING SEQUENCE HAVE THE FOLLOWING MEANING *
C* *
C* N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS TO BE *
C* INTEGRATED. N MAY BE DECREASED ON SUBSEQUENT CALLS IF *
C* THE NUMBER OF ACTIVE EQUATIONS DECREASES, WITH THE *
C* FIRST ONES BEING THOSE RETAINED. BUT, IT MUST *
C* NOT BE INCREASED WITHOUT CALLING WITH JSTART = 0. *
C* T THE INDEPENDENT VARIABLE. ON ENTRY T IS THE CURRENT *
C* SETTING OF THE INDEPENDENT VARIABLE. ON RETURN TO THE *
C* CALLING PROGRAM, T CORRESPONDS TO THE SETTING OF THE *
C* INDEPENDENT VARIABLE FOR THE MOST FORWARD POINT OBTAINED*
C* THUS FAR. *
C* Y AN N BY 8 BY 4 ARRAY CONTAINING THE DEPENDENT VARIABLES *
C* AND THEIR BACKWARD DIFFERENCES. ON EACH CALL UP TO FOUR*
C* SOLUTION POINTS ARE OBTAINED. THE MOST FORWARD POINT IS *
C* ALWAYS AT Y(I,1,1). THE POINT NEXT TO THE MOST FORWARD *
C* POINT IS RETURNED IN Y(I,1,KFLAG). THE MOST BACKWARD *
C* POINT IN THE NEW BLOCK IS RETURNED IN Y(I,1,2). ONLY THE*
C* INITIAL SOLUTION VALUES, ENTERED IN Y(I,1,1) FOR *
C* I=1,...,N, NEED TO BE SUPPLIED ON THE FIRST *
C* CALL (JSTART = 0). *
C* Y(I,J+1,K) CONTAINS THE J-TH BACKWARD DIFFERENCE OF THE *
C* I-TH DEPENDENT VARIABLE (FOR K=1,...,KFLAG). *
C* IF IT IS DESIRED TO INTERPOLATE TO NON-MESH POINTS, *
C* THESE VALUES CAN BE USED. IF THE CURRENT STEP-SIZE IS *
C* H AND THE VALUE AT T-E (0.LT.DABS(E).LT.DABS(H)) IS *
C* NEEDED, FORM S=E/H AND THEN COMPUTE *
C* NQ *
C* Y(I)(T-E) = SUM Y(I,J+1,K)*B(J) *
C* J=0 *
C* WHERE K, WHICH CORRESPONDS TO T, IS THE FIRST MESH POINT*
C* BEYOND THE POINT T-E, AND B(J+1) = B(J)*(J-S)/(J+1) *
C* WITH B(0) = 1. *
C* DY AN N BY 4 ARRAY. DY(I,K), 1.LE.I.LE.N, 1.LE.K.LE.KFLAG, *
C* CONTAINS THE DERIVATIVE OF THE I-TH DEPENDENT VARIABLE *
C* SCALED BY H. THE DY(I,1) ARRAY NEED NOT BE SUPPLIED AT *
C* THE FIRST CALL (JSTART = 0). *
C* SAVE A BLOCK OF AT LEAST 13*N0 DOUBLE PRECISION FLOATING POINT*
C* LOCATIONS. *
C* H THE STEP-SIZE USED FOR THE JUST COMPLETED BLOCK. *
C* HNEXT THE STEP-SIZE FOR THE NEXT BLOCK. ON THE *
C* FIRST CALL (JSTART=0) THE USER MUST SPECIFY AN INITIAL *
C* STEP-SIZE. A GOOD ESTIMATE OF ITS MAGNITUDE IS GIVEN *
C* BY 0.2*(EPS/DABS(E))**0.5, WHERE E IS THE LARGEST *
C* EIGENVALUE OF THE JACOBIAN EVALUATED AT THE INITIAL *
C* VALUES OF T AND Y. THE SIGN OF THE INITIAL HNEXT *
C* SHOULD BE POSITIVE (NEGATIVE) IF THE FINAL TIME IS *
C* GREATER (LESS) THAN THE INITIAL TIME. IF THE INITIAL *
C* STEP-SIZE CHOICE DOES NOT CAUSE AN ERROR GREATER THAN *
C* EPS, IN THE RMS NORM, IT WILL BE ACCEPTED. OTHERWISE, *
C* ITS MAGNITUDE WILL BE DECREASED UNTIL AN ERROR LESS *
C* THAN EPS IS ACHIEVED. THE SUBROUTINE AUTOMATICALLY *
C* ADJUSTS THE STEP-SIZE AFTER THE INITIAL AND SUBSEQUENT *
C* CALLS FOR THE STEP-SIZE OF LARGEST POSSIBLE MAGNITUDE. *
C* THE MAGNITUDE MAY BE ADJUSTED DOWN ON ANY SUBSEQUENT *
C* CALL. NOTE--A MAGNITUDE ADJUSTMENT UP OR A SIGN CHANGE *
C* WILL BE IGNORED. *
C* HMIN DABS(HMIN) IS THE MINIMUM STEP-SIZE MAGNITUDE TO BE *
C* ALLOWED FOR THE NEXT INTEGRATION CYCLE. (ON THE FIRST *
C* CALL (JSTART=0), DABS(HMIN) SHOULD BE CHOSEN SIGNIF- *
C* ICANTLY SMALLER THAN DABS(HNEXT) SO AS TO AVOID START-UP*
C* PROBLEMS IF THE ERROR CRITERION IS NOT MET WITH THE USER*
C* SPECIFIED HNEXT.) NOTE--THE SIGN OF HMIN IS IGNORED. *
C* HMIN MAY BE CHANGED ON SUBSEQUENT CALLS. *
C* HMAX DABS(HMAX) IS THE MAXIMUM STEP-SIZE MAGNITUDE TO BE *
C* ALLOWED FOR THE NEXT INTEGRATION CYCLE. NOTE--IF *
C* DABS(HMAX) IS LESS THAN DABS(HMIN), THEN THE SUBROUTINE *
C* FUNCTIONS AS THOUGH DABS(HMAX) EQUALS DABS(HMIN). HMAX *
C* MAY BE CHANGED ON SUBSEQUENT CALLS. *
C* EPS THE ERROR TEST CONSTANT. THE SINGLE STEP ERROR ESTIMATE *
C* FOR Y, COMPUTED AS A WEIGHTED RMS NORM, MUST NOT EXCEED *
C* EPS. THE WEIGHT FOR THE I-TH ELEMENT IN THE ERROR *
C* ESTIMATE IS 1/YMAX(I). (SEE PARAMETER YMAX.) THE *
C* STEP-SIZE AND/OR ORDER ARE ADJUSTED TO ACHIEVE THIS. *
C* YMAX AN ARRAY OF N0 LOCATIONS, WITH--FOR I=1,...,N--YMAX(I) *
C* BEING THE MAXIMUM OF UNITY AND THE MAXIMUM VALUE OF *
C* DABS(Y(I)) SEEN THUS FAR. ON THE FIRST CALL IT SHOULD *
C* BE SET TO THE MAXIMUM OF UNITY AND THE INITIAL VALUE OF *
C* DABS(Y(I)). *
C* KFLAG A COMPLETION CODE. IF KFLAG IS GREATER THAN 0,THEN *
C* KFLAG POINTS HAVE BEEN COMPUTED. IF KFLAG IS *
C* -2, -3, OR -4 THEN 2, 3, OR 4 MESH POINTS, RESPECTIVELY,*
C* HAVE BEEN COMPUTED WITH DABS(H) EQUAL TO DABS(HMIN), *
C* BUT THE REQUESTED ERROR WAS NOT ACHIEVED. OTHER VALUES *
C* KFLAG CAN ASSUME ARE AS FOLLOWS *
C* -5 THE REQUESTED ERROR WAS SMALLER THAN CAN BE HANDLED*
C* FOR THIS PROBLEM. *
C* -6 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED FOR *
C* DABS(H).GT.DABS(HMIN). *
C* -7 THE MAXIMUM ORDER SPECIFIED WAS TOO LARGE. *
C* KNEXT AFTER THE INITIAL CALL (JSTART=0), THE VALUE OF KNEXT IS *
C* THE NUMBER OF POINTS TO BE COMPUTED DURING THE NEXT *
C* CYCLE. THE VALUE IS SUPPLIED FOR INFORMATION ONLY. THE*
C* USER CANNOT CONTROL THE NUMBER OF POINTS BY SETTING *
C* KNEXT. *
C* JSTART AN INPUT INDICATOR WITH THE FOLLOWING MEANINGS *
C* .EQ.0 INITIALIZATION CALL. (JSTART MUST BE SET TO 0 *
C* ON THE FIRST CALL.) *
C* .GT.0 CONTINUE FROM THE LAST STEP. *
C* ON RETURN JSTART IS SET TO NQ, THE MAXIMUM BACKWARD *
C* DIFFERENCE AVAILABLE IN THE Y ARRAY. (THIS ALSO CORRES-*
C* PONDS TO THE ORDER OF THE METHOD USED TO COMPUTE THE *
C* JUST COMPLETED BLOCK OF POINTS.) *
C* MAXORD THE MAXIMUM ORDER (1.LE.MAXORD.LE.7) THAT MAY BE USED. *
C* NOTE--IF MAXORD IS RESET BETWEEN CYCLES TO A VALUE LESS *
C* THAN THE ORDER DETERMINED FOR THE NEXT CYCLE, THEN THE *
C* ORDER MAY FOR SEVERAL CYCLES EXCEED MAXORD. HOWEVER, *
C* IT CANNOT EXCEED THE ABOVE DETERMINED VALUE AND, ONCE *
C* THE ORDER IS LESS THAN OR EQUAL TO MAXORD, IT CANNOT *
C* THEN EXCEED MAXORD. *
C* RJ A BLOCK OF AT LEAST N0**2 DOUBLE PRECISION FLOATING POINT*
C* LOCATIONS, WHICH CONTAINS AN ESTIMATE OF THE JACOBIAN *
C* OF THE DIFFERENTIAL EQUATION. *
C* RW A BLOCK OF AT LEAST N0**2 DOUBLE PRECISION FLOATING POINT*
C* LOCATIONS. *
C* IPIV A BLOCK OF AT LEAST N0 INTEGERS USED TO HOLD PIVOT DATA *
C* GENERATED DURING AN LU DECOMPOSITION. *
C* MF METHOD FLAG. IT DETERMINES THE MODE BY WHICH THE PARTIAL*
C* DERIVATIVES ARE OBTAINED WITH THE FOLLOWING MEANINGS *
C* 1 ANALYTIC EXPRESSIONS FOR THE PARTIAL DERIVATIVES ARE *
C* SUPPLIED BY THE USER IN THE SUBROUTINE PEDERV. *
C* 2 THE ANALYTIC EXPRESSIONS ARE NOT AVAILABLE. *
C* APPROXIMATE VALUES OF THE PARTIAL DERIVATIVES ARE *
C* OBTAINED BY NUMERICAL DIFFERENCING. *
C* *
C***********************************************************************
C
DOUBLE PRECISION B, BND, B1, B2, C, CRATE, D, DF, DI, DY,
+ D1, D2, D3, E, EDOWN, ENQDWN, ENQSAM,
+ ENQUP, EPS, ES, EUP, FN, H, HMAX, HMIN,
+ HNEW, HNEXT, HOLD, Q, RATIO, RC, RJ, RMAX,
+ RRDOWN, RRSAME, RRUP, RW, SAVE, SQRTUR, T,
+ TDL, TOLD, UROUND, Y, YJ1, YMAX
C DOUBLE PRECISION DABS, DBLE, DMAX1, DMIN1, DSQRT
DIMENSION Y(N,8,4), DY(N,4), SAVE(9,1), YMAX(1), RW(1), IPIV(1),
+ INDEX(7,2), B(82), C(16), PERROR(9), PC(16), PD(7),
+ B1(48), B2(34), RJ(1)
COMMON /STAT/ NSTEP,NFE,NJE,NINVS
EQUIVALENCE(B(1),B1(1)),(B(49),B2(1))
C
C ****************************************************
C ***** SEVEN OF EIGHT DATA STATEMENTS CONTAIN *****
C ***** ARRAY NAMES AND VARIOUS LINES OF CODE *****
C ***** INCLUDE SUBSCRIPT EXPRESSIONS. *****
C ****************************************************
C
C***********************************************************************
C* THE CONSTANT UROUND SHOULD BE SET EQUAL TO THE UNIT ROUND-OFF *
C* FOR THE MACHINE. THE CONSTANT SQRTUR SHOULD BE SET EQUAL TO THE *
C* SQUARE ROOT OF UROUND. *
C***********************************************************************
C
DATA UROUND/4.D-16/, SQRTUR/2.D-8/
C
C***********************************************************************
C* THE ARRAY INDEX HOLDS POINTERS AND CONSTANTS FOR THE VARIOUS *
C* ORDER METHODS. FOR NQ=1,...,7 THE ENTRIES ARE AS FOLLOWS *
C* INDEX(NQ,1) BASE INDEX FOR B ARRAY (H*DY PREDICTOR). *
C* INDEX(NQ,2) BASE INDEX FOR C ARRAY (CORRECTOR). *
C***********************************************************************
C
DATA INDEX/ 1, 2, 4, 11, 20, 38, 59,
+ 1, 2, 3, 5, 7, 10, 14/
C
C***********************************************************************
C* THE COEFFICIENTS APPEARING IN THE NEXT DATA STATEMENTS FOR THE *
C* B ARRAY SHOULD BE DEFINED TO THE MAXIMUM ACCURACY PERMITTED BY *
C* THE MACHINE. THEY ARE, IN THE ORDER SPECIFIED,... *
C* 1 *
C* -2, 3 *
C* -9/2, -5/4, 11/2 *
C* -15/2, -3/4, 13/2, 2 *
C* -22/3, -8/3, -17/18, 25/3 *
C* -9, -9/4, -7/8, 35/4, 5/4 *
C* -125/12, -101/24, -71/36, -37/48, 137/12 *
C* -253/24, -1001/240, -707/360, -123/160, 1373/120, 1/10 *
C* -57/4, -25/8, -17/10, -169/240, 61/5, -1/20, 31/10 *
C* -137/10, -117/20, -46/15, -191/120, -197/300, 147/10 *
C* -353/25, -571/100, -1819/600, -79/50, -3919/6000, 1477/100, 7/20 *
C* -3529/200, -1889/400, -1609/600, -3527/2400, -931/1500, 3079/200,*
C* -3/20, 17/5 *
C* -343/20, -303/40, -253/60, -589/240, -101/75, -23/40, 363/20 *
C* -266/15, -221/30, -749/180, -73/30, -803/600, -103/180, 547/30, *
C* 1/2 *
C* -1316/75, -1151/150, -3689/900, -121/50, -4003/3000, -257/450, *
C* 2737/150, -1/5, 1/2 *
C***********************************************************************
C
DATA B1/
+ 1.0D0, -2.0D0, 3.0D0, -4.5D0, -1.25D0, 5.5D0, -7.5D0,
+ -.75D0, 6.5D0, 2.0D0, -7.3333333333333333D0,
+ -2.6666666666666667D0, -.94444444444444444D0,
+ 8.3333333333333333D0, -9.0D0, -2.25D0, -.875D0, 8.75D0,
+ 1.25D0, -10.416666666666667D0, -4.2083333333333333D0,
+ -1.9722222222222222D0, -.77083333333333333D0,
+ 11.416666666666667D0, -10.541666666666667D0,
+ -4.1708333333333333D0, -1.9638888888888889D0, -.76875D0,
+ 11.441666666666667D0, 0.1D0, -14.25D0, -3.125D0, -1.7D0,
+ -.70416666666666667D0, 12.2D0, -.05D0, 3.1D0, -13.7D0,
+ -5.85D0, -3.0666666666666667D0, -1.5916666666666667D0,
+ -.65666666666666667D0, 14.7D0, -14.12D0, -5.71D0,
+ -3.0316666666666667D0, -1.58D0, -.65316666666666667D0/
DATA B2/
+ 14.77D0, 0.35D0, -17.645D0, -4.7225D0,
+ -2.6816666666666667D0, -1.4695833333333333D0,
+ -.62066666666666667D0, 15.395D0, -.15D0, 3.4D0, -17.15D0,
+ -7.575D0, -4.2166666666666667D0, -2.4541666666666667D0,
+ -1.3466666666666667D0, -.575D0, 18.15D0,
+ -17.733333333333333D0, -7.3666666666666667D0,
+ -4.1611111111111111D0, -2.4333333333333333D0,
+ -1.3383333333333333D0, -.57222222222222222D0,
+ 18.233333333333333D0, 0.5D0, -17.546666666666667D0,
+ -7.6733333333333333D0, -4.0988888888888889D0, -2.42D0,
+ -1.3343333333333333D0, -.57111111111111111D0,
+ 18.246666666666667D0, -0.2D0, 0.5D0/
C
C***********************************************************************
C* THE COEFFICIENTS APPEARING IN THE NEXT DATA STATEMENTS FOR THE *
C* C ARRAY SHOULD BE DEFINED TO THE MAXIMUM ACCURACY PERMITTED BY *
C* THE MACHINE. THEY ARE, IN THE ORDER SPECIFIED,... *
C* -1 *
C* -2/3 *
C* -6/11, -2/3 *
C* -12/25, -16/31 *
C* -60/137, -600/1373, -100/193 *
C* -20/49, -120/293, -75/184, -1200/2299 *
C* -140/363, -60/143, -1050/2437 *
C***********************************************************************
C
DATA C/ -1.0D0, -.66666666666666667D0, -.54545454545454545D0,
+ -.66666666666666667D0, -.48D0, -.51612903225806452D0,
+ -.43795620437956204D0, -.43699927166788056D0,
+ -.51813471502590674D0, -.40816326530612245D0,
+ -.40955631399317406D0, -.40760869565217391D0,
+ -.52196607220530666D0, -.38567493112947659D0,
+ -.41958041958041958D0, -.43085761181780879D0/
C
C***********************************************************************
C* THE COEFFICIENTS IN THE PERROR ARRAY ARE USED IN THE ERROR TEST, *
C* THE FIRST TIME IT IS PERFORMED, AS WELL AS IN THE STEP-SIZE/ORDER *
C* SELECTION SEGMENT. PERROR(I+1) = 1/D(I+1), I=1,...,7, WHERE *
C* D(I+1) IS THE DISCRETIZATION ERROR CONSTANT CORRESPONDING TO THE *
C* SECOND PASS OF THE INTEGRATION CYCLE OF ORDER I. PERROR(1) AND *
C* PERROR(9) ARE DEFINED SOLELY FOR PROGRAMMING EASE. THEY ARE NOT *
C* USED. *
C***********************************************************************
C
DATA PERROR/ 1.0, 1.0, 1.92857, 2.78161, 3.56735,
+ 4.29497, 4.9065, 5.6066, 1.0/
C
C***********************************************************************
C* THE COEFFICIENTS IN THE ARRAY PC ARE USED BOTH IN THE CONVERGENCE *
C* AND ERROR TESTS. THEY ARE THE RECIPROCAL VALUES OF THE *
C* DISCRETIZATION ERROR CONSTANTS FOR EQUATIONS CONSTITUTING THE *
C* METHODS OF ORDER 1 THRU 7. *
C***********************************************************************
C
DATA PC/ 2.0, 4.5, 7.3333, 6.0, 10.4167, 9.3, 13.7, 13.8687,
+ 9.6904, 17.15, 16.9504, 17.4349, 9.472, 20.7429, 15.921,
+ 14.7809/
C
C***********************************************************************
C* THE COEFFICIENTS APPEARING IN THE ARRAY PD ARE USED IN THE TESTING*
C* OF THE *OUTDATEDNESS* OF THE ARRAY RW. THE PD(I) ELEMENT CONTAINS*
C* THE AVERAGE VALUE OF THE COEFFICIENTS IN ARRAY C CORRESPONDING TO *
C* ORDER I. *
C***********************************************************************
C
DATA PD/ 1.0, 0.6667, 0.6061, 0.4981, 0.4644, 0.4368, 0.412/
C
C***********************************************************************
C***********************************************************************
C
KFLAG = 0
IFAIL = 0
IF(JSTART.NE.0) GO TO 80
C* *
C* INITIALIZATION --- FIRST CALL *
C* *
20 N0 = N
FN=DBLE(FLOAT(N))
N9 = N0 * 9
N9P1 = N9 + 1
N10 = N9 + N0
N11 = N10 + N0
N11P1 = N11 + 1
N12 = N11 + N0
N12P1 = N12 + 1
NSQ = N0*N0
IWEVAL=+1
TDL=T
H=DMAX1(DABS(HMIN),DMIN1(DABS(HNEXT),DABS(HMAX)))
IF(HNEXT.LT.0.D0) H=-H
30 NQOLD= 0
ISW1=0
ISW2=0
CRATE=1.D0
CALL DIFFUN (N, T, Y, DY)
NFE=NFE+1
DO 40 I=1,N
40 DY(I,1) = H*DY(I,1)
NQ = 1
IST = 1
IDEL = 0
GO TO 180
C* *
C* CONTINUE WITH THE STEP-SIZE H *
C* *
80 HNEW=DMAX1(DABS(HMIN),DMIN1(DABS(HNEW),DABS(HNEXT),DABS(HMAX)))
IF(H.LT.0.D0) HNEW=-HNEW
IF(H.NE.HNEW) GO TO 100
IST = NQP1
IDEL = NQP1
ISW2=0
GO TO 180
C* *
C* NEW STEP-SIZE --- INTERPOLATE FOR NEW POINTS *
C* *
100 RATIO=HNEW/H
H = RATIO*H
RC=RC*RATIO*DBLE(PD(NQ)/PDOLD)
PDOLD = PD(NQ)
IF (NQ.EQ.1) GO TO 150
IEQ = NEQ
D = 0.0D0
DO 140 J=2,NQ
D = D + RATIO
IF (D.GT.DBLE(FLOAT(NEQ+NQP1-IEQ))) IEQ=2
D1 = DBLE(FLOAT(NEQ-IEQ-1)) - D
DO 140 I=1,N
D2 = 1.0D0
D3 = 0.0D0
DO 120 J1=2,NQP1
D2 = D2*(DBLE(FLOAT(J1)) + D1)/DBLE(FLOAT(J1-1))
120 D3 = D3 + D2*Y(I,J1,IEQ)
140 Y(I,J,1) = D3 + Y(I,1,IEQ)
150 IST = NQ
IDEL = 0
DO 160 I=1,N
160 DY(I,1) = DY(I,1)*RATIO
IRET = 3
GO TO 4000
C* *
C* INITIALIZE SAVE ARRAY *
C* *
180 DO 200 I=1,N
SAVE(9,I) = DY(I,1)
DO 200 J=1,IST
200 SAVE(J,I) = Y(I,J,1)
NQST = NQ
RATIO = 1.0D0
TOLD = T
HOLD = H
220 IF ((NQ.EQ.NQOLD).AND.(FN.EQ.DBLE(FLOAT(N)))) GO TO 300
IF ((NQ.EQ.NQOLD).AND.(FN.NE.DBLE(FLOAT(N)))) GO TO 280
IF (MAXORD.LT.8) GO TO 260
KFLAG = -7
GO TO 4185
C* *
C* SET APPROPRIATE PARAMETERS AND CONSTANTS *
C* FOR NEW CYCLE OF ORDER NQ *
C* *
260 INDB = INDEX(NQ,1)
INDC = INDEX(NQ,2)
NEQ = 3+NQ/5
JSTART = NQ
NQOLD = NQ
NQM1=NQ-1
NQP1 = NQ+1
Q = DBLE(FLOAT(NQ))
ENQDWN = 0.5D0/Q
ENQSAM = 0.5D0/(Q+1.0D0)
ENQUP = 0.5D0/(Q+2.0D0)
280 FN = DBLE(FLOAT(N))
EDOWN = FN*(DBLE(PERROR(NQ))*EPS)**2
EUP = FN*(DBLE(PERROR(NQ+2))*EPS)**2
ES = FN*(DBLE(PERROR(NQP1))*EPS)**2
IF(EDOWN.GT.0.D0) GO TO 300
KFLAG = -5
GO TO 4185
C* *
C* CHECK FOR REEVALUATION OF JACOBIAN *
C* *
300 IF(IWEVAL.GT.0) GO TO 320
IF(DABS(RC-1.D0).GE.4.D-1) IWEVAL = +2
IF((TDL-TOLD)*H.GT.0.D0) GO TO 320
IF(DABS(RC-1.D0).GE.8.D-1) IWEVAL = +1
320 INDBB = INDB
INDCC = INDC
DO 1000 IEQ=1,NEQ
IST = MOD (IEQ, NEQ) + 1
T = T + H
BND=FN*(DBLE(PC(INDCC))*ENQUP*EPS)**2
E=ES
IF(IEQ.GT.2) E=FN*(DBLE(PC(INDCC))*EPS)**2
C* *
C* PREDICT Y AND DY FOR THE NEXT MESH POINT *
C* *
DO 460 I=1,N
D = DY(I,IEQ)
D1 = Y(I,1,IEQ) + Q*D
D2 = B(INDBB+NQM1)*D
IF (NQ-2) 440, 400, 340
340 IF (IEQ-3) 400, 380, 360
360 D2 = D2 + B(INDBB+NQP1)*DY(I,3)
380 D2 = D2 + B(INDBB+NQ)*DY(I,2)
400 DO 420 J=2,NQ
D = Y(I,J,IEQ)
D3 = DBLE(FLOAT(J-1))
D1 = D1 + D*(D3-Q)/D3
420 D2 = D2 + D*B(INDBB+J-2)
440 Y(I,1,IST) = D1
460 DY(I,IST) = D2
IF(NQ.LE.2) GO TO 470
IF(IEQ.GT.1) INDBB = INDBB + NQ + IEQ - 2
C* *
C* ITERATE THE CORRECTOR UP TO THREE TIMES. ACCUMULATE THE *
C* CORRECTION TERMS IN SAVE(N10+I,1) FOR REDOING THE ENTIRE *
C* CORRECTOR IF CONVERGENCE IS NOT ACHIEVED *
C* *
470 D1 = C(INDCC)
500 DO 510 I=1,N
510 SAVE(N10+I,1) = 0.D0
DO 780 J=1,3
CALL DIFFUN (N, T, Y(1,1,IST), SAVE(N9P1,1))
NFE=NFE+1
IF(IWEVAL-1) 680, 520, 620
520 IND = 1
IF(IEQ.EQ.2) IND = 1
GO TO (610, 530), MF
530 CALL DIFFUN (N, T-H*DBLE(FLOAT(1+IEQ-IND)), Y(1,1,IND),
+ SAVE(N12P1,1))
NFE=NFE+1
D = 0.D0
DO 540 I=1,N
540 D = D + SAVE(N12+I,1)**2
D = DABS(H)*1.D3*UROUND*DSQRT(D)
NJE = NJE + 1
J0 = 0
DO 560 J1=1,N
DI = SQRTUR*YMAX(J1)
DI = DMAX1(DI,D)
YJ1 = Y(J1,1,IND)
Y(J1,1,IND) = Y(J1,1,IND) + DI
CALL DIFFUN (N, T-H*DBLE(FLOAT(1+IEQ-IND)), Y(1,1,IND),
+ SAVE(N11P1,1))
NFE = NFE + 1
DO 550 I=1,N
550 RJ(I+J0) = (SAVE(N11+I,1) - SAVE(N12+I,1))/DI
J0 = J0 + N0
Y(J1,1,IND) = YJ1
560 CONTINUE
GO TO 620
610 CALL PEDERV(N, T-H*DBLE(FLOAT(1+IEQ-IND)), Y(1,1,IND), RJ, N0)
NJE=NJE+1
620 D = D1*H
DO 640 I=1,NSQ
640 RW(I) = RJ(I)*D
J1 = 1
DO 660 I=1,N
RW(J1) = 1.D0 + RW(J1)
660 J1 = J1 + N0 + 1
CALL DEC (N, N0, RW, IPIV, IER)
NINVS=NINVS+1
IWEVAL = -IEQ
RC=1.D0
PDOLD=PD(NQ)
IF(IER.NE.0) GO TO 800
680 DO 700 I=1,N
700 SAVE(N9+I,1) = DY(I,IST) - H*SAVE(N9+I,1)
CALL SOL (N, N0, RW, SAVE(N9P1,1), IPIV)
D2 = 0.0D0
J3 = N9
DO 760 I=1,N
J3 = J3 + 1
J2 = J3 + N0
SAVE(J2,1) = SAVE(J2,1) + SAVE(J3,1)
Y(I,1,IST) = Y(I,1,IST) + D1*SAVE(J3,1)
DY(I,IST) = DY(I,IST) - SAVE(J3,1)
760 D2 = D2 + (SAVE(J3,1)/YMAX(I))**2
IF(J.NE.1) CRATE = DMAX1(CRATE*9.D-1,D2/D3)
IF (D2*DMIN1(1.D0,2.D0*CRATE).GT.BND) GO TO 770
GO TO 940
770 D3 = D2
780 CONTINUE
C* *
C* CORRECTOR FAILED TO CONVERGE IN THREE ITERATIONS. *
C* IF JACOBIAN WAS REEVALUATED DURING THIS CYCLE, STEP-SIZE IS *
C* REDUCED TO 3/10 OF H. OTHERWISE, JACOBIAN IS REEVALUATED *
C* *
TDL = TOLD
IF(IWEVAL.EQ.0) GO TO 900
800 IF(DABS(H).LE.(1.00001D0*DABS(HMIN))) IF(NQ-2) 4180, 1080, 4140
RATIO = RATIO*0.3D0
IRET = 1
ISW1=0
ISW2=1
IWEVAL = +2
GO TO 3000
900 DO 920 I=1,N
D = SAVE(N10+I,1)
Y(I,1,IST) = Y(I,1,IST) - D1*D
920 DY(I,IST) = DY(I,IST) + D
IWEVAL = +1
GO TO 500
C* *
C* CORRECTOR CONVERGED. THE BACKWARD DIFFERENCES OF ORDER *
C* ONE THROUGH NQ AT THE NEW MESH POINT ARE COMPUTED *
C* *
940 DO 980 I=1,N
DO 960 J=1,NQ
960 Y(I,J+1,IST) = Y(I,J,IST) - Y(I,J,IEQ)
IF(IEQ.EQ.1) GO TO 980
C* *
C* THE (NQ+1)-ST BACKWARD DIFFERENCE FOR ALL BUT THE FIRST *
C* MESH POINT IS ESTABLISHED *
C* *
SAVE(N9+I,1) = Y(I,NQP1,IST) - Y(I,NQP1,IEQ)
980 CONTINUE
IF (IEQ.EQ.NEQ) GO TO 1015
IF (NQ.LE.2) GO TO 1010
IF ((IEQ.GT.1).OR.(NQ.EQ.6)) INDCC = INDCC + 1
1010 IF(IEQ.EQ.1) GO TO 1000
C* *
C* ERROR TEST FOR ALL BUT THE FIRST MESH POINT IS PERFORMED *
C* *
1015 D = 0.D0
DO 1020 I=1,N
1020 D = D + (SAVE(N9+I,1)/YMAX(I))**2
IF(D.GT.E) GO TO 1030
1000 CONTINUE
GO TO 1160
C* *
C* THE ERROR CRITERION WAS NOT MET *
C* *
1030 IFAIL=IFAIL+1
IF(IFAIL.GT.2) GO TO 1080
IF(DABS(H).LE.(DABS(HMIN)*1.00001D0)) IF(NQ-2) 1140, 4120, 4140
IWEVAL = +2
IF(IFAIL.EQ.1) GO TO 1200
TDL=T
RATIO = RATIO*5.D-1
IRET = 1
ISW1=0
ISW2=1
GO TO 3000
C* *
C* START OVER WITH ORDER 1 METHOD *
C* *
1080 IFAIL = 0
DO 1090 I=1,N
1090 Y(I,1,1) = SAVE(1,I)
T = TOLD
H = H*DMAX1(1.D-1,DABS(HMIN/H))
IWEVAL = +2
GO TO 30
C* *
C* A NEW BLOCK OF SOLUTION POINTS HAS BEEN ACCEPTED *
C* *
1140 IWEVAL = 0
NSTEP = NSTEP + IEQ
KFLAG = -IEQ
RETURN
1160 IWEVAL = 0
E = ES
KFLAG = NEQ
NSTEP = NSTEP + KFLAG
HNEW=H
C* *
C* CHECK FOR CONTINUATION WITH THE SAME H AND NQ *
C* *
IF(ISW2.EQ.1) GO TO 1420
IF(NQ.GT.3) ISW1=1-ISW1
IF(ISW1.EQ.1) GO TO 1420
C* *
C* NEW STEP-SIZE AND/OR ORDER SELECTION *
C* *
1200 RRSAME=1.2D0*(D/E)**ENQSAM
IF(IFAIL.NE.0) GO TO 1380
RMAX = 1.D-4
DF=DBLE(FLOAT(NEQ+NQM1))
IF(NQ.NE.1) RMAX=(Q-1.D0)/DF
RRSAME=DMAX1(RRSAME,RMAX)
RRUP = 1.D20
RRDOWN = 1.D20
IF(NQ.GE.MAXORD) GO TO 1240
D = 0.0D0
DO 1220 I=1,N
D1 = Y(I,NQP1,NEQ) - Y(I,NQP1,NEQ-1)
1220 D = D + ((SAVE(N9+I,1) - D1)/YMAX(I))**2
RRUP = 1.2D0*(D/EUP)**ENQUP
RMAX=Q/DF
RRUP=DMAX1(RRUP,RMAX)
1240 IF(NQ.EQ.1) GO TO 1280
D = 0.0D0
DO 1260 I=1,N
1260 D = D + (Y(I,NQP1,1)/YMAX(I))**2
RRDOWN = 1.2D0*(D/EDOWN)**ENQDWN
RMAX = 1.D-4
IF(NQ.NE.2) RMAX=(Q-2.D0)/DF
RRDOWN=DMAX1(RRDOWN,RMAX)
1280 IF (RRSAME.GT.RRUP) IF (RRUP-RRDOWN) 1340, 1300, 1300
IF (RRSAME.LE.RRDOWN) GO TO 1320
1300 NEWQ = NQM1
D=1.D0/RRDOWN
GO TO 1360
1320 NEWQ = NQ
D=1.D0/RRSAME
GO TO 1360
1340 NEWQ = NQP1
D=1.D0/RRUP
1360 IF (D-1.1D0) 1420, 1400, 1400
1380 RATIO=RATIO/RRSAME
IRET = 1
ISW1=0
ISW2=1
GO TO 3000
1400 HNEW = H*D
NQ = NEWQ
1420 DO 1460 I=1,N
D=YMAX(I)
DO 1440 J=1,NEQ
1440 D=DMAX1(D,DABS(Y(I,1,J)))
1460 YMAX(I)=D
HNEXT = HNEW
KNEXT = 3+NQ/5
RETURN
C* *
C* THIS SECTION IS USED WHEN STEP-SIZE OR ORDER IS CHANGED *
C* DURING THE CYCLE. STARTING VALUES ARE RETRIEVED FROM THE *
C* SAVE ARRAY. *
C* *
3000 RATIO = DMAX1 (DABS(HMIN/HOLD), DMIN1(RATIO, 1.0D0))
T = TOLD
IF (RATIO.LT.1.0D0) IF (IDEL) 3030, 3030, 3080
DO 3020 I=1,N
DY(I,1) = SAVE(9,I)
DO 3020 J=1,NQ
3020 Y(I,J,1) = SAVE(J,I)
GO TO (320,260), IRET
C* *
C* THE (NQST+1)-ST ORDER BACKWARD DIFFERENCE IS ESTABLISHED *
C* *
3030 IDEL = NQST+1
IF ((NQST.LT.2) .OR. (NQ.LT.2)) GO TO 3080
D = DBLE(FLOAT(NQST))
DO 3060 I=1,N
D1 = SAVE(9,I)
DO 3040 J=2,NQST
3040 D1 = D1 - SAVE(J,I)/DBLE(FLOAT(J-1))
3060 SAVE(IDEL,I) = D*D1
C* *
C* INTERPOLATE FOR NEW POINTS *
C* *
3080 DO 3140 I=1,N
DY(I,1) = RATIO*SAVE(9,I)
IF (NQ.LT.2) GO TO 3140
D1 = 0.0D0
DO 3120 J=2,NQ
D1 = D1 + RATIO
D2 = 1.0D0
D = 0.0D0
DO 3100 J1=2,IDEL
D2 = D2*(DBLE(FLOAT(J1-2)) - D1)/DBLE(FLOAT(J1-1))
3100 D = D + D2*SAVE(J1,I)
3120 Y(I,J,1) = D + SAVE(1,I)
3140 Y(I,1,1) = SAVE(1,I)
H = HOLD*RATIO
C* *
C* FORM THE BACKWARD DIFFERENCES *
C* *
4000 IF (NQ.LT.2) GO TO 4040
NQM1 = NQ-1
DO 4020 I=1,N
DO 4020 J=1,NQM1
J0=J + 1
DO 4020 J1=J0,NQ
J2 = NQ-J1+J+1
4020 Y(I,J2,1) = Y(I,J2-1,1) - Y(I,J2,1)
4040 GO TO (320,260,180), IRET
C* *
C* CONVERGENCE OR ERROR PROBLEMS *
C* *
4120 NQ = 1
GO TO 4160
4140 NQ = 2
4160 IFAIL = 0
IRET = 2
IWEVAL = +2
GO TO 3000
4180 KFLAG = -6
4185 J1 = NQST + 1
DO 4190 I=1,N
DY(I,1) = SAVE(9,I)
DO 4190 J=1,J1
4190 Y(I,J,1) = SAVE(J,I)
H = HOLD
T = TOLD
JSTART = NQST
RETURN
C
C***********************************************************************
C* *
C* --- END OF SUBROUTINE STINT --- *
C* *
C***********************************************************************
C
END
SUBROUTINE DEC (N, NDIM, A, IP, IER) DEC0 10
C
INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
DOUBLE PRECISION A, T
C DOUBLE PRECISION DABS
DIMENSION A(NDIM,N), IP(N)
C-----------------------------------------------------------------------
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C INPUT..
C N = ORDER OF MATRIX.
C NDIM = DECLARED DIMENSION OF ARRAY A .
C A = MATRIX TO BE TRIANGULARIZED.
C OUTPUT..
C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
C SINGULAR AT STAGE K.
C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
C
C REFERENCE..
C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
C C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
IER = 0
IP(N) = 1
IF (N .EQ. 1) GO TO 70
NM1 = N - 1
DO 60 K = 1,NM1
KP1 = K + 1
M = K
DO 10 I = KP1,N
10 IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
IP(K) = M
T = A(M,K)
IF (M .EQ. K) GO TO 20
IP(N) = -IP(N)
A(M,K) = A(K,K)
A(K,K) = T
20 IF (T .EQ. 0.D0) GO TO 80
T = 1.D0/T
DO 30 I = KP1,N
30 A(I,K) = -A(I,K)*T
DO 50 J = KP1,N
T = A(M,J)
A(M,J) = A(K,J)
A(K,J) = T
IF (T .EQ. 0.D0) GO TO 50
DO 40 I = KP1,N
40 A(I,J) = A(I,J) + A(I,K)*T
50 CONTINUE
60 CONTINUE
70 K = N
IF (A(N,N) .EQ. 0.D0) GO TO 80
RETURN
80 IER = K
IP(N) = 0
RETURN
C----------------------- END OF SUBROUTINE DEC -------------------------
END
SUBROUTINE SOL (N, NDIM, A, B, IP) SOL0 10
C
INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
DOUBLE PRECISION A,B,T
DIMENSION A(NDIM,N), B(N), IP(N)
C-----------------------------------------------------------------------
C SOLUTION OF LINEAR SYSTEM, A*X = B .
C INPUT..
C N = ORDER OF MATRIX.
C NDIM = DECLARED DIMENSION OF ARRAY A .
C A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
C B = RIGHT HAND SIDE VECTOR.
C IP = PIVOT VECTOR OBTAINED FROM DEC.
C DO NOT USE IF DEC HAS SET IER .NE. 0.
C OUTPUT..
C B = SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
IF (N .EQ. 1) GO TO 50
NM1 = N - 1
DO 20 K = 1,NM1
KP1 = K + 1
M = IP(K)
T = B(M)
B(M) = B(K)
B(K) = T
DO 10 I = KP1,N
10 B(I) = B(I) + A(I,K)*T
20 CONTINUE
DO 40 KB = 1,NM1
KM1 = N - KB
K = KM1 + 1
B(K) = B(K)/A(K,K)
T = -B(K)
DO 30 I = 1,KM1
30 B(I) = B(I) + A(I,K)*T
40 CONTINUE
50 B(1) = B(1)/A(1,1)
RETURN
C----------------------- END OF SUBROUTINE SOL -------------------------
END
SUBROUTINE DIFFUN(N, T, Y, YDOT) TEST 10
DOUBLE PRECISION T, T1, T2, W1, W2, W3, W4, S, Y, YDOT
DIMENSION Y(N), YDOT(N)
T1=-Y(1)+Y(2)+Y(3)+Y(4)
T2=Y(1)-Y(2)+Y(3)+Y(4)
W1=5.D-1*(T1*T1-T2*T2)
W2=T1*T2
W3=(Y(1)+Y(2)-Y(3)+Y(4))**2
W4=(Y(1)+Y(2)+Y(3)-Y(4))**2
S=0.D0
DO 10 I=1,3
10 S=S-1.D-3*Y(I)
S=S+1.D-3*Y(4)
YDOT(1)=-2.5D-1*(9.8D2*Y(1)+1.02D3*Y(2)-9.8D2*Y(3)+1.02D3*Y(4)-S)
+ +1.25D-1*(-W1+W2+W3+W4)
YDOT(2)=-2.5D-1*(1.02D3*Y(1)+9.8D2*Y(2)-1.02D3*Y(3)+9.8D2*Y(4)-S)
+ +1.25D-1*( W1-W2+W3+W4)
YDOT(3)=-2.5D-1*(-1.02D3*Y(1)-9.8D2*Y(2)+9.8D2*Y(3)-1.02D3*Y(4)-S)
+ +1.25D-1*( W1+W2-W3+W4)
YDOT(4)=-2.5D-1*(9.8D2*Y(1)+1.02D3*Y(2)-1.02D3*Y(3)+9.8D2*Y(4)+S)
+ +1.25D-1*( W1+W2+W3-W4)
RETURN
END
SUBROUTINE PEDERV(N, T, Y, P0, N0) TEST 240
DOUBLE PRECISION P0, Y, T
DIMENSION P0(16), Y(N)
P0(1)=(-9.80001D2+Y(1)+Y(3)+Y(4)+3.D0*Y(2))*2.5D-1
P0(5)=(-1.020001D3+Y(2)-Y(3)-Y(4)+3.D0*Y(1))*2.5D-1
P0(9)=(9.79999D2+Y(1)-Y(2)-Y(4)+3.D0*Y(3))*2.5D-1
P0(13)=(-1.019999D3+Y(1)-Y(2)-Y(3)+3.D0*Y(4))*2.5D-1
P0(2)=P0(5)
P0(6)=P0(1)
P0(10)=-P0(13)
P0(14)=-P0(9)
P0(3)=P0(10)
P0(7)=P0(9)
P0(11)=P0(1)
P0(15)=-P0(5)
P0(4)=P0(14)
P0(8)=P0(13)
P0(12)=P0(15)
P0(16)=P0(1)
RETURN
END
SUBROUTINE ERROR(T, Y, N, ER0) TEST 450
DOUBLE PRECISION T, T1, T2, T3, W1, W2, W3, W4, X, Y, ER0
C DOUBLE PRECISION DEXP, DSIN, DCOS, DABS
DIMENSION ER0(N), Y(N)
X=1.D1*T
IF(X.GT.1.65D2) GO TO 20
T1=1.D0+DEXP(-X)*(9.D0*DCOS(X)+1.D1*DSIN(X))
T2=DEXP(-X)*(1.D1*DCOS(X)-9.0D0*DSIN(X))
IF(DABS(T2).LE.1.D-15) GO TO 10
T3=1.D1/(T1*T1+T2*T2)
GO TO 11
10 T3=1.D1/T1**2
11 W1=-(T1+T2)*T3
W2=(T1-T2)*T3
GO TO 21
20 W1=-1.D1
W2=1.D1
21 X=1.D3*T
IF(X.GT.1.65D2) GO TO 22
W3=5.D2/(1.D0-1.001D3*DEXP(X))
GO TO 23
22 W3=0.D0
23 W4=5.D-4/(1.D0-1.001D0*DEXP(1.D-3*T))
ER0(1)=-W1+W2+W3+W4
ER0(2)= W1-W2+W3+W4
ER0(3)= W1+W2-W3+W4
ER0(4)= W1+W2+W3-W4
DO 30 I=1,N
30 ER0(I)=ER0(I)-Y(I)
RETURN
END
20 TEST 760
6 TEST 770
15 TEST 780
04 TEST 790
0.000000000000000D+00-2.000000000000000D+00-1.000000000000000D+00 TEST 800
-1.000000000000000D+00 TEST 810
0.000000000000000D+00 1.000000000000000D+02 5.000000000000000D-09 TEST 820
1 TEST 830
1.000000000000000D-03 TEST 840
1 TEST 850
1.000000000000000D-05 TEST 860
1 TEST 870
1.000000000000000D-07 TEST 880
2 TEST 890
1.000000000000000D-03 TEST 900
2 TEST 910
1.000000000000000D-05 TEST 920
2 TEST 930
1.000000000000000D-07 TEST 940
.