C This file contains the internal documentation only for DVERK. The C documentation also is relevant to the complete suite of codes developed C with the same design and calling sequence as DVERK. In this latter case C the calling sequence is intepreted in the same way, but the underlying C strategies may have been modified. C C SUBROUTINE DVERK(N, FCN, X, Y, XEND, TOL, IND, C, NW, W) C C*********************************************************************** C * C PURPOSE - THIS IS A RUNGE-KUTTA SUBROUTINE BASED ON VERNER'S * C FIFTH AND SIXTH ORDER PAIR OF FORMULAS FOR FINDING APPROXIMATIONS TO * C THE SOLUTION OF A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL * C EQUATIONS WITH INITIAL CONDITIONS. IT ATTEMPTS TO KEEP THE GLOBAL * C ERROR PROPORTIONAL TO A TOLERANCE SPECIFIED BY THE USER. (THE * C PROPORTIONALITY DEPENDS ON THE KIND OF ERROR CONTROL THAT IS USED, * C AS WELL AS THE DIFFERENTIAL EQUATION AND THE RANGE OF INTEGRATION.) * C * C VARIOUS OPTIONS ARE AVAILABLE TO THE USER, INCLUDING DIFFERENT * C KINDS OF ERROR CONTROL, RESTRICTIONS ON STEP SIZES, AND INTERRUPTS * C WHICH PERMIT THE USER TO EXAMINE THE STATE OF THE CALCULATION (AND * C PERHAPS MAKE MODIFICATIONS) DURING INTERMEDIATE STAGES. * C * C THE PROGRAM IS EFFICIENT FOR NON-STIFF SYSTEMS. HOWEVER, A GOOD * C VARIABLE-ORDER-ADAMS METHOD WILL PROBABLY BE MORE EFFICIENT IF THE * C FUNCTION EVALUATIONS ARE VERY COSTLY. SUCH A METHOD WOULD ALSO BE * C MORE SUITABLE IF ONE WANTED TO OBTAIN A LARGE NUMBER OF INTERMEDIATE * C SOLUTION VALUES BY INTERPOLATION, AS MIGHT BE THE CASE FOR EXAMPLE * C WITH GRAPHICAL OUTPUT. * C * C HULL-ENRIGHT-JACKSON 1/10/76 * C * C*********************************************************************** C * C USE - THE USER MUST SPECIFY EACH OF THE FOLLOWING * C * C N NUMBER OF EQUATIONS * C * C FCN NAME OF SUBROUTINE FOR EVALUATING FUNCTIONS - THE SUBROUTINE * C ITSELF MUST ALSO BE PROVIDED BY THE USER - IT SHOULD BE OF * C THE FOLLOWING FORM * C SUBROUTINE FCN(N, X, Y, YPRIME) * C INTEGER N * C DOUBLE PRECISION X, Y(N), YPRIME(N) * C *** ETC *** * C AND IT SHOULD EVALUATE YPRIME, GIVEN N, X AND Y * C * C X INDEPENDENT VARIABLE - INITIAL VALUE SUPPLIED BY USER * C * C Y DEPENDENT VARIABLE - INITIAL VALUES OF COMPONENTS Y(1), Y(2), * C ..., Y(N) SUPPLIED BY USER * C * C XEND VALUE OF X TO WHICH INTEGRATION IS TO BE CARRIED OUT - IT MAY * C BE LESS THAN THE INITIAL VALUE OF X * C * C TOL TOLERANCE - THE SUBROUTINE ATTEMPTS TO CONTROL A NORM OF THE * C LOCAL ERROR IN SUCH A WAY THAT THE GLOBAL ERROR IS * C PROPORTIONAL TO TOL. IN SOME PROBLEMS THERE WILL BE ENOUGH * C DAMPING OF ERRORS, AS WELL AS SOME CANCELLATION, SO THAT * C THE GLOBAL ERROR WILL BE LESS THAN TOL. ALTERNATIVELY, THE * C CONTROL CAN BE VIEWED AS ATTEMPTING TO PROVIDE A * C CALCULATED VALUE OF Y AT XEND WHICH IS THE EXACT SOLUTION * C TO THE PROBLEM Y' = F(X,Y) + E(X) WHERE THE NORM OF E(X) * C IS PROPORTIONAL TO TOL. (THE NORM IS A MAX NORM WITH * C WEIGHTS THAT DEPEND ON THE ERROR CONTROL STRATEGY CHOSEN * C BY THE USER. THE DEFAULT WEIGHT FOR THE K-TH COMPONENT IS * C 1/MAX(1,ABS(Y(K))), WHICH THEREFORE PROVIDES A MIXTURE OF * C ABSOLUTE AND RELATIVE ERROR CONTROL.) * C * C IND INDICATOR - ON INITIAL ENTRY IND MUST BE SET EQUAL TO EITHER * C 1 OR 2. IF THE USER DOES NOT WISH TO USE ANY OPTIONS, HE * C SHOULD SET IND TO 1 - ALL THAT REMAINS FOR THE USER TO DO * C THEN IS TO DECLARE C AND W, AND TO SPECIFY NW. THE USER * C MAY ALSO SELECT VARIOUS OPTIONS ON INITIAL ENTRY BY * C SETTING IND = 2 AND INITIALIZING THE FIRST 9 COMPONENTS OF * C C AS DESCRIBED IN THE NEXT SECTION. HE MAY ALSO RE-ENTER * C THE SUBROUTINE WITH IND = 3 AS MENTIONED AGAIN BELOW. IN * C ANY EVENT, THE SUBROUTINE RETURNS WITH IND EQUAL TO * C 3 AFTER A NORMAL RETURN * C 4, 5, OR 6 AFTER AN INTERRUPT (SEE OPTIONS C(8), C(9)) * C -1, -2, OR -3 AFTER AN ERROR CONDITION (SEE BELOW) * C * C C COMMUNICATIONS VECTOR - THE DIMENSION MUST BE GREATER THAN OR * C EQUAL TO 24, UNLESS OPTION C(1) = 4 OR 5 IS USED, IN WHICH * C CASE THE DIMENSION MUST BE GREATER THAN OR EQUAL TO N+30 * C * C NW FIRST DIMENSION OF WORKSPACE W - MUST BE GREATER THAN OR * C EQUAL TO N * C * C W WORKSPACE MATRIX - FIRST DIMENSION MUST BE NW AND SECOND MUST * C BE GREATER THAN OR EQUAL TO 9 * C * C THE SUBROUTINE WILL NORMALLY RETURN WITH IND = 3, HAVING * C REPLACED THE INITIAL VALUES OF X AND Y WITH, RESPECTIVELY, THE VALUE * C OF XEND AND AN APPROXIMATION TO Y AT XEND. THE SUBROUTINE CAN BE * C CALLED REPEATEDLY WITH NEW VALUES OF XEND WITHOUT HAVING TO CHANGE * C ANY OTHER ARGUMENT. HOWEVER, CHANGES IN TOL, OR ANY OF THE OPTIONS * C DESCRIBED BELOW, MAY ALSO BE MADE ON SUCH A RE-ENTRY IF DESIRED. * C * C THREE ERROR RETURNS ARE ALSO POSSIBLE, IN WHICH CASE X AND Y * C WILL BE THE MOST RECENTLY ACCEPTED VALUES - * C WITH IND = -3 THE SUBROUTINE WAS UNABLE TO SATISFY THE ERROR * C REQUIREMENT WITH A PARTICULAR STEP-SIZE THAT IS LESS THAN OR * C EQUAL TO HMIN, WHICH MAY MEAN THAT TOL IS TOO SMALL * C WITH IND = -2 THE VALUE OF HMIN IS GREATER THAN HMAX, WHICH * C PROBABLY MEANS THAT THE REQUESTED TOL (WHICH IS USED IN THE * C CALCULATION OF HMIN) IS TOO SMALL * C WITH IND = -1 THE ALLOWED MAXIMUM NUMBER OF FCN EVALUATIONS HAS * C BEEN EXCEEDED, BUT THIS CAN ONLY OCCUR IF OPTION C(7), AS * C DESCRIBED IN THE NEXT SECTION, HAS BEEN USED * C * C THERE ARE SEVERAL CIRCUMSTANCES THAT WILL CAUSE THE CALCULATIONS * C TO BE TERMINATED, ALONG WITH OUTPUT OF INFORMATION THAT WILL HELP * C THE USER DETERMINE THE CAUSE OF THE TROUBLE. THESE CIRCUMSTANCES * C INVOLVE ENTRY WITH ILLEGAL OR INCONSISTENT VALUES OF THE ARGUMENTS, * C SUCH AS ATTEMPTING A NORMAL RE-ENTRY WITHOUT FIRST CHANGING THE * C VALUE OF XEND, OR ATTEMPTING TO RE-ENTER WITH IND LESS THAN ZERO. * C * C*********************************************************************** C * C OPTIONS - IF THE SUBROUTINE IS ENTERED WITH IND = 1, THE FIRST 9 * C COMPONENTS OF THE COMMUNICATIONS VECTOR ARE INITIALIZED TO ZERO, AND * C THE SUBROUTINE USES ONLY DEFAULT VALUES FOR EACH OPTION. IF THE * C SUBROUTINE IS ENTERED WITH IND = 2, THE USER MUST SPECIFY EACH OF * C THESE 9 COMPONENTS - NORMALLY HE WOULD FIRST SET THEM ALL TO ZERO, * C AND THEN MAKE NON-ZERO THOSE THAT CORRESPOND TO THE PARTICULAR * C OPTIONS HE WISHES TO SELECT. IN ANY EVENT, OPTIONS MAY BE CHANGED ON * C RE-ENTRY TO THE SUBROUTINE - BUT IF THE USER CHANGES ANY OF THE * C OPTIONS, OR TOL, IN THE COURSE OF A CALCULATION HE SHOULD BE CAREFUL * C ABOUT HOW SUCH CHANGES AFFECT THE SUBROUTINE - IT MAY BE BETTER TO * C RESTART WITH IND = 1 OR 2. (COMPONENTS 10 TO 24 OF C ARE USED BY THE * C PROGRAM - THE INFORMATION IS AVAILABLE TO THE USER, BUT SHOULD NOT * C NORMALLY BE CHANGED BY HIM.) * C * C C(1) ERROR CONTROL INDICATOR - THE NORM OF THE LOCAL ERROR IS THE * C MAX NORM OF THE WEIGHTED ERROR ESTIMATE VECTOR, THE * C WEIGHTS BEING DETERMINED ACCORDING TO THE VALUE OF C(1) - * C IF C(1)=1 THE WEIGHTS ARE 1 (ABSOLUTE ERROR CONTROL) * C IF C(1)=2 THE WEIGHTS ARE 1/ABS(Y(K)) (RELATIVE ERROR * C CONTROL) * C IF C(1)=3 THE WEIGHTS ARE 1/MAX(ABS(C(2)),ABS(Y(K))) * C (RELATIVE ERROR CONTROL, UNLESS ABS(Y(K)) IS LESS * C THAN THE FLOOR VALUE, ABS(C(2)) ) * C IF C(1)=4 THE WEIGHTS ARE 1/MAX(ABS(C(K+30)),ABS(Y(K))) * C (HERE INDIVIDUAL FLOOR VALUES ARE USED) * C IF C(1)=5 THE WEIGHTS ARE 1/ABS(C(K+30)) * C FOR ALL OTHER VALUES OF C(1), INCLUDING C(1) = 0, THE * C DEFAULT VALUES OF THE WEIGHTS ARE TAKEN TO BE * C 1/MAX(1,ABS(Y(K))), AS MENTIONED EARLIER * C (IN THE TWO CASES C(1) = 4 OR 5 THE USER MUST DECLARE THE * C DIMENSION OF C TO BE AT LEAST N+30 AND MUST INITIALIZE THE * C COMPONENTS C(31), C(32), ..., C(N+30).) * C * C C(2) FLOOR VALUE - USED WHEN THE INDICATOR C(1) HAS THE VALUE 3 * C * C C(3) HMIN SPECIFICATION - IF NOT ZERO, THE SUBROUTINE CHOOSES HMIN * C TO BE ABS(C(3)) - OTHERWISE IT USES THE DEFAULT VALUE * C 10*MAX(DWARF,RREB*MAX(WEIGHTED NORM Y/TOL,ABS(X))), * C WHERE DWARF IS A VERY SMALL POSITIVE MACHINE NUMBER AND * C RREB IS THE RELATIVE ROUNDOFF ERROR BOUND * C * C C(4) HSTART SPECIFICATION - IF NOT ZERO, THE SUBROUTINE WILL USE * C AN INITIAL HMAG EQUAL TO ABS(C(4)), EXCEPT OF COURSE FOR * C THE RESTRICTIONS IMPOSED BY HMIN AND HMAX - OTHERWISE IT * C USES THE DEFAULT VALUE OF HMAX*(TOL)**(1/6) * C * C C(5) SCALE SPECIFICATION - THIS IS INTENDED TO BE A MEASURE OF THE * C SCALE OF THE PROBLEM - LARGER VALUES OF SCALE TEND TO MAKE * C THE METHOD MORE RELIABLE, FIRST BY POSSIBLY RESTRICTING * C HMAX (AS DESCRIBED BELOW) AND SECOND, BY TIGHTENING THE * C ACCEPTANCE REQUIREMENT - IF C(5) IS ZERO, A DEFAULT VALUE * C OF 1 IS USED. FOR LINEAR HOMOGENEOUS PROBLEMS WITH * C CONSTANT COEFFICIENTS, AN APPROPRIATE VALUE FOR SCALE IS A * C NORM OF THE ASSOCIATED MATRIX. FOR OTHER PROBLEMS, AN * C APPROXIMATION TO AN AVERAGE VALUE OF A NORM OF THE * C JACOBIAN ALONG THE TRAJECTORY MAY BE APPROPRIATE * C * C C(6) HMAX SPECIFICATION - FOUR CASES ARE POSSIBLE * C IF C(6).NE.0 AND C(5).NE.0, HMAX IS TAKEN TO BE * C MIN(ABS(C(6)),2/ABS(C(5))) * C IF C(6).NE.0 AND C(5).EQ.0, HMAX IS TAKEN TO BE ABS(C(6)) * C IF C(6).EQ.0 AND C(5).NE.0, HMAX IS TAKEN TO BE * C 2/ABS(C(5)) * C IF C(6).EQ.0 AND C(5).EQ.0, HMAX IS GIVEN A DEFAULT VALUE * C OF 2 * C * C C(7) MAXIMUM NUMBER OF FUNCTION EVALUATIONS - IF NOT ZERO, AN * C ERROR RETURN WITH IND = -1 WILL BE CAUSED WHEN THE NUMBER * C OF FUNCTION EVALUATIONS EXCEEDS ABS(C(7)) * C * C C(8) INTERRUPT NUMBER 1 - IF NOT ZERO, THE SUBROUTINE WILL * C INTERRUPT THE CALCULATIONS AFTER IT HAS CHOSEN ITS * C PRELIMINARY VALUE OF HMAG, AND JUST BEFORE CHOOSING HTRIAL * C AND XTRIAL IN PREPARATION FOR TAKING A STEP (HTRIAL MAY * C DIFFER FROM HMAG IN SIGN, AND MAY REQUIRE ADJUSTMENT IF * C XEND IS NEAR) - THE SUBROUTINE RETURNS WITH IND = 4, AND * C WILL RESUME CALCULATION AT THE POINT OF INTERRUPTION IF * C RE-ENTERED WITH IND = 4 * C * C C(9) INTERRUPT NUMBER 2 - IF NOT ZERO, THE SUBROUTINE WILL * C INTERRUPT THE CALCULATIONS IMMEDIATELY AFTER IT HAS * C DECIDED WHETHER OR NOT TO ACCEPT THE RESULT OF THE MOST * C RECENT TRIAL STEP, WITH IND = 5 IF IT PLANS TO ACCEPT, OR * C IND = 6 IF IT PLANS TO REJECT - Y(*) IS THE PREVIOUSLY * C ACCEPTED RESULT, WHILE W(*,9) IS THE NEWLY COMPUTED TRIAL * C VALUE, AND W(*,2) IS THE UNWEIGHTED ERROR ESTIMATE VECTOR. * C THE SUBROUTINE WILL RESUME CALCULATIONS AT THE POINT OF * C INTERRUPTION ON RE-ENTRY WITH IND = 5 OR 6. (THE USER MAY * C CHANGE IND IN THIS CASE IF HE WISHES, FOR EXAMPLE TO FORCE * C ACCEPTANCE OF A STEP THAT WOULD OTHERWISE BE REJECTED, OR * C VICE VERSA. HE CAN ALSO RESTART WITH IND = 1 OR 2.) * C * C*********************************************************************** C * C SUMMARY OF THE COMPONENTS OF THE COMMUNICATIONS VECTOR * C * C PRESCRIBED AT THE OPTION DETERMINED BY THE PROGRAM * C OF THE USER * C * C C(10) RREB(REL ROUNDOFF ERR BND) * C C(1) ERROR CONTROL INDICATOR C(11) DWARF (VERY SMALL MACH NO) * C C(2) FLOOR VALUE C(12) WEIGHTED NORM Y * C C(3) HMIN SPECIFICATION C(13) HMIN * C C(4) HSTART SPECIFICATION C(14) HMAG * C C(5) SCALE SPECIFICATION C(15) SCALE * C C(6) HMAX SPECIFICATION C(16) HMAX * C C(7) MAX NO OF FCN EVALS C(17) XTRIAL * C C(8) INTERRUPT NO 1 C(18) HTRIAL * C C(9) INTERRUPT NO 2 C(19) EST * C C(20) PREVIOUS XEND * C C(21) FLAG FOR XEND * C C(22) NO OF SUCCESSFUL STEPS * C C(23) NO OF SUCCESSIVE FAILURES * C C(24) NO OF FCN EVALS * C * C IF C(1) = 4 OR 5, C(31), C(32), ... C(N+30) ARE FLOOR VALUES * C * C RREB and DWARF are machine dependent constants currently set so * C that they should be appropriate for most machines. However, it may * C be appropriate to change them when this program is installed on a * C new machine. K.R.J. 3 Oct 1991. * C * C*********************************************************************** C * C AN OVERVIEW OF THE PROGRAM * C * C BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES * C ......ABORT IF IND OUT OF RANGE 1 TO 6 * C . CASES - INITIAL ENTRY, NORMAL RE-ENTRY, INTERRUPT RE-ENTRIES * C . CASE 1 - INITIAL ENTRY (IND .EQ. 1 OR 2) * C V........ABORT IF N.GT.NW OR TOL.LE.0 * C . IF INITIAL ENTRY WITHOUT OPTIONS (IND .EQ. 1) * C . SET C(1) TO C(9) EQUAL TO ZERO * C . ELSE INITIAL ENTRY WITH OPTIONS (IND .EQ. 2) * C . MAKE C(1) TO C(9) NON-NEGATIVE * C . MAKE FLOOR VALUES NON-NEGATIVE IF THEY ARE TO BE USED * C . END IF * C . INITIALIZE RREB, DWARF, PREV XEND, FLAG, COUNTS * C . CASE 2 - NORMAL RE-ENTRY (IND .EQ. 3) * C .........ABORT IF XEND REACHED, AND EITHER X CHANGED OR XEND NOT * C . RE-INITIALIZE FLAG * C . CASE 3 - RE-ENTRY FOLLOWING AN INTERRUPT (IND .EQ. 4 TO 6) * C V TRANSFER CONTROL TO THE APPROPRIATE RE-ENTRY POINT....... * C . END CASES . * C . END INITIALIZATION, ETC. . * C . V * C . LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL STEP . * C . STAGE 1 - PREPARE . * C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT . * C . CALC SLOPE (ADDING 1 TO NO OF FCN EVALS) IF IND .NE. 6 . * C . CALC HMIN, SCALE, HMAX . * C***********ERROR RETURN (WITH IND=-2) IF HMIN .GT. HMAX . * C . CALC PRELIMINARY HMAG . * C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED.......RE-ENTRY.V * C . CALC HMAG, XTRIAL AND HTRIAL . * C . END STAGE 1 . * C V STAGE 2 - CALC YTRIAL (ADDING 7 TO NO OF FCN EVALS) . * C . STAGE 3 - CALC THE ERROR ESTIMATE . * C . STAGE 4 - MAKE DECISIONS . * C . SET IND=5 IF STEP ACCEPTABLE, ELSE SET IND=6 . * C***********INTERRUPT NO 2 IF REQUESTED....................RE-ENTRY.V * C . IF STEP ACCEPTED (IND .EQ. 5) * C . UPDATE X, Y FROM XTRIAL, YTRIAL * C . ADD 1 TO NO OF SUCCESSFUL STEPS * C . SET NO OF SUCCESSIVE FAILURES TO ZERO * C**************RETURN(WITH IND=3, XEND SAVED, FLAG SET) IF X .EQ. XEND * C . ELSE STEP NOT ACCEPTED (IND .EQ. 6) * C . ADD 1 TO NO OF SUCCESSIVE FAILURES * C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN * C . END IF * C . END STAGE 4 * C . END LOOP * C . * C BEGIN ABORT ACTION * C OUTPUT APPROPRIATE MESSAGE ABOUT STOPPING THE CALCULATIONS, * C ALONG WITH VALUES OF IND, N, NW, TOL, HMIN, HMAX, X, XEND, * C PREVIOUS XEND, NO OF SUCCESSFUL STEPS, NO OF SUCCESSIVE * C FAILURES, NO OF FCN EVALS, AND THE COMPONENTS OF Y * C STOP * C END ABORT ACTION * C * C***********************************************************************