SUBROUTINE DDEABM(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,RWORK,LRW, * IWORK,LIW,RPAR,IPAR) C***BEGIN PROLOGUE DDEABM C***DATE WRITTEN 820301 (YYMMDD) C***REVISION DATE 851111 (YYMMDD) C***CATEGORY NO. I1A1B C***KEYWORDS ADAMS METHOD,DEPAC,DOUBLE PRECISION, C INITIAL VALUE PROBLEMS,ODE, C ORDINARY DIFFERENTIAL EQUATIONS,PREDICTOR-CORRECTOR C***AUTHOR SHAMPINE, L. F., (SNLA) C WATTS, H. A., (SNLA) C***PURPOSE SOLVE INITIAL VALUE PROBLEMS IN ORDINARY DIFFERENTIAL C EQUATIONS USING AN ADAMS-BASHFORTH METHOD. C***DESCRIPTION C C ****** DOUBLE PRECISION VERSION OF DEABM ****** C C THIS IS THE ADAMS CODE IN THE PACKAGE OF DIFFERENTIAL EQUATION C SOLVERS DEPAC, CONSISTING OF THE CODES DDERKF, DDEABM, AND DDEBDF. C DESIGN OF THE PACKAGE WAS BY L. F. SHAMPINE AND H. A. WATTS. C IT IS DOCUMENTED IN C SAND79-2374 , DEPAC - DESIGN OF A USER ORIENTED PACKAGE OF ODE C SOLVERS. C DDEABM IS A DRIVER FOR A MODIFICATION OF THE CODE ODE WRITTEN BY C L. F. SHAMPINE AND M. K. GORDON C SANDIA LABORATORIES C ALBUQUERQUE, NEW MEXICO 87185 C C ********************************************************************** C * ABSTRACT * C ************ C C SUBROUTINE DDEABM USES THE ADAMS-BASHFORTH-MOULTON C PREDICTOR-CORRECTOR FORMULAS OF ORDERS ONE THROUGH TWELVE TO C INTEGRATE A SYSTEM OF NEQ FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS OF THE FORM C DU/DX = DF(X,U) C WHEN THE VECTOR Y(*) OF INITIAL VALUES FOR U(*) AT X=T IS GIVEN. C THE SUBROUTINE INTEGRATES FROM T TO TOUT. IT IS EASY TO CONTINUE THE C INTEGRATION TO GET RESULTS AT ADDITIONAL TOUT. THIS IS THE INTERVAL C MODE OF OPERATION. IT IS ALSO EASY FOR THE ROUTINE TO RETURN WITH C THE SOLUTION AT EACH INTERMEDIATE STEP ON THE WAY TO TOUT. THIS IS C THE INTERMEDIATE-OUTPUT MODE OF OPERATION. C C DDEABM USES SUBPROGRAMS DDES, DSTEPS, DINTP, DHSTRT, DVNORM, C D1MACH, AND THE ERROR HANDLING ROUTINE XERRWV. THE ONLY MACHINE C DEPENDENT PARAMETERS TO BE ASSIGNED APPEAR IN D1MACH. C C ********************************************************************** C * DESCRIPTION OF THE ARGUMENTS TO DDEABM (AN OVERVIEW) * C ********************************************************************** C C THE PARAMETERS ARE C C DF -- THIS IS THE NAME OF A SUBROUTINE WHICH YOU PROVIDE TO C DEFINE THE DIFFERENTIAL EQUATIONS. C C NEQ -- THIS IS THE NUMBER OF (FIRST ORDER) DIFFERENTIAL C EQUATIONS TO BE INTEGRATED. C C T -- THIS IS A VALUE OF THE INDEPENDENT VARIABLE. C C Y(*) -- THIS ARRAY CONTAINS THE SOLUTION COMPONENTS AT T. C C TOUT -- THIS IS A POINT AT WHICH A SOLUTION IS DESIRED. C C INFO(*) -- THE BASIC TASK OF THE CODE IS TO INTEGRATE THE C DIFFERENTIAL EQUATIONS FROM T TO TOUT AND RETURN AN C ANSWER AT TOUT. INFO(*) IS AN INTEGER ARRAY WHICH IS USED C TO COMMUNICATE EXACTLY HOW YOU WANT THIS TASK TO BE C CARRIED OUT. C C RTOL, ATOL -- THESE QUANTITIES REPRESENT RELATIVE AND ABSOLUTE C ERROR TOLERANCES WHICH YOU PROVIDE TO INDICATE HOW C ACCURATELY YOU WISH THE SOLUTION TO BE COMPUTED. YOU MAY C CHOOSE THEM TO BE BOTH SCALARS OR ELSE BOTH VECTORS. C C IDID -- THIS SCALAR QUANTITY IS AN INDICATOR REPORTING WHAT C THE CODE DID. YOU MUST MONITOR THIS INTEGER VARIABLE TO C DECIDE WHAT ACTION TO TAKE NEXT. C C RWORK(*), LRW -- RWORK(*) IS A DOUBLE PRECISION WORK ARRAY OF C LENGTH LRW WHICH PROVIDES THE CODE WITH NEEDED STORAGE C SPACE. C C IWORK(*), LIW -- IWORK(*) IS AN INTEGER WORK ARRAY OF LENGTH LIW C WHICH PROVIDES THE CODE WITH NEEDED STORAGE SPACE. C C RPAR, IPAR -- THESE ARE DOUBLE PRECISION AND INTEGER PARAMETER C ARRAYS WHICH YOU CAN USE FOR COMMUNICATION BETWEEN YOUR C CALLING PROGRAM AND THE DF SUBROUTINE. C C QUANTITIES WHICH ARE USED AS INPUT ITEMS ARE C NEQ, T, Y(*), TOUT, INFO(*), C RTOL, ATOL, RWORK(1), LRW AND LIW. C C QUANTITIES WHICH MAY BE ALTERED BY THE CODE ARE C T, Y(*), INFO(1), RTOL, ATOL, C IDID, RWORK(*) AND IWORK(*). C C ********************************************************************** C * INPUT -- WHAT TO DO ON THE FIRST CALL TO DDEABM * C ********************************************************************** C C THE FIRST CALL OF THE CODE IS DEFINED TO BE THE START OF EACH NEW C PROBLEM. READ THROUGH THE DESCRIPTIONS OF ALL THE FOLLOWING ITEMS, C PROVIDE SUFFICIENT STORAGE SPACE FOR DESIGNATED ARRAYS, SET C APPROPRIATE VARIABLES FOR THE INITIALIZATION OF THE PROBLEM, AND C GIVE INFORMATION ABOUT HOW YOU WANT THE PROBLEM TO BE SOLVED. C C C DF -- PROVIDE A SUBROUTINE OF THE FORM C DF(X,U,UPRIME,RPAR,IPAR) C TO DEFINE THE SYSTEM OF FIRST ORDER DIFFERENTIAL EQUATIONS C WHICH IS TO BE SOLVED. FOR THE GIVEN VALUES OF X AND THE C VECTOR U(*)=(U(1),U(2),...,U(NEQ)) , THE SUBROUTINE MUST C EVALUATE THE NEQ COMPONENTS OF THE SYSTEM OF DIFFERENTIAL C EQUATIONS DU/DX=DF(X,U) AND STORE THE DERIVATIVES IN THE C ARRAY UPRIME(*), THAT IS, UPRIME(I) = * DU(I)/DX * FOR C EQUATIONS I=1,...,NEQ. C C SUBROUTINE DF MUST NOT ALTER X OR U(*). YOU MUST DECLARE C THE NAME DF IN AN EXTERNAL STATEMENT IN YOUR PROGRAM THAT C CALLS DDEABM. YOU MUST DIMENSION U AND UPRIME IN DF. C C RPAR AND IPAR ARE DOUBLE PRECISION AND INTEGER PARAMETER C ARRAYS WHICH YOU CAN USE FOR COMMUNICATION BETWEEN YOUR C CALLING PROGRAM AND SUBROUTINE DF. THEY ARE NOT USED OR C ALTERED BY DDEABM. IF YOU DO NOT NEED RPAR OR IPAR, C IGNORE THESE PARAMETERS BY TREATING THEM AS DUMMY C ARGUMENTS. IF YOU DO CHOOSE TO USE THEM, DIMENSION THEM IN C YOUR CALLING PROGRAM AND IN DF AS ARRAYS OF APPROPRIATE C LENGTH. C C NEQ -- SET IT TO THE NUMBER OF DIFFERENTIAL EQUATIONS. C (NEQ .GE. 1) C C T -- SET IT TO THE INITIAL POINT OF THE INTEGRATION. C YOU MUST USE A PROGRAM VARIABLE FOR T BECAUSE THE CODE C CHANGES ITS VALUE. C C Y(*) -- SET THIS VECTOR TO THE INITIAL VALUES OF THE NEQ SOLUTION C COMPONENTS AT THE INITIAL POINT. YOU MUST DIMENSION Y AT C LEAST NEQ IN YOUR CALLING PROGRAM. C C TOUT -- SET IT TO THE FIRST POINT AT WHICH A SOLUTION C IS DESIRED. YOU CAN TAKE TOUT = T, IN WHICH CASE THE CODE C WILL EVALUATE THE DERIVATIVE OF THE SOLUTION AT T AND C RETURN. INTEGRATION EITHER FORWARD IN T (TOUT .GT. T) OR C BACKWARD IN T (TOUT .LT. T) IS PERMITTED. C C THE CODE ADVANCES THE SOLUTION FROM T TO TOUT USING C STEP SIZES WHICH ARE AUTOMATICALLY SELECTED SO AS TO C ACHIEVE THE DESIRED ACCURACY. IF YOU WISH, THE CODE WILL C RETURN WITH THE SOLUTION AND ITS DERIVATIVE FOLLOWING C EACH INTERMEDIATE STEP (INTERMEDIATE-OUTPUT MODE) SO THAT C YOU CAN MONITOR THEM, BUT YOU STILL MUST PROVIDE TOUT IN C ACCORD WITH THE BASIC AIM OF THE CODE. C C THE FIRST STEP TAKEN BY THE CODE IS A CRITICAL ONE C BECAUSE IT MUST REFLECT HOW FAST THE SOLUTION CHANGES NEAR C THE INITIAL POINT. THE CODE AUTOMATICALLY SELECTS AN C INITIAL STEP SIZE WHICH IS PRACTICALLY ALWAYS SUITABLE FOR C THE PROBLEM. BY USING THE FACT THAT THE CODE WILL NOT STEP C PAST TOUT IN THE FIRST STEP, YOU COULD, IF NECESSARY, C RESTRICT THE LENGTH OF THE INITIAL STEP SIZE. C C FOR SOME PROBLEMS IT MAY NOT BE PERMISSIBLE TO INTEGRATE C PAST A POINT TSTOP BECAUSE A DISCONTINUITY OCCURS THERE C OR THE SOLUTION OR ITS DERIVATIVE IS NOT DEFINED BEYOND C TSTOP. WHEN YOU HAVE DECLARED A TSTOP POINT (SEE INFO(4) C AND RWORK(1)), YOU HAVE TOLD THE CODE NOT TO INTEGRATE C PAST TSTOP. IN THIS CASE ANY TOUT BEYOND TSTOP IS INVALID C INPUT. C C INFO(*) -- USE THE INFO ARRAY TO GIVE THE CODE MORE DETAILS ABOUT C HOW YOU WANT YOUR PROBLEM SOLVED. THIS ARRAY SHOULD BE C DIMENSIONED OF LENGTH 15 TO ACCOMODATE OTHER MEMBERS OF C DEPAC OR POSSIBLE FUTURE EXTENSIONS, THOUGH DDEABM USES C ONLY THE FIRST FOUR ENTRIES. YOU MUST RESPOND TO ALL OF C THE FOLLOWING ITEMS WHICH ARE ARRANGED AS QUESTIONS. THE C SIMPLEST USE OF THE CODE CORRESPONDS TO ANSWERING ALL C QUESTIONS AS YES ,I.E. SETTING ALL ENTRIES OF INFO TO 0. C C INFO(1) -- THIS PARAMETER ENABLES THE CODE TO INITIALIZE C ITSELF. YOU MUST SET IT TO INDICATE THE START OF EVERY C NEW PROBLEM. C C **** IS THIS THE FIRST CALL FOR THIS PROBLEM ... C YES -- SET INFO(1) = 0 C NO -- NOT APPLICABLE HERE. C SEE BELOW FOR CONTINUATION CALLS. **** C C INFO(2) -- HOW MUCH ACCURACY YOU WANT OF YOUR SOLUTION C IS SPECIFIED BY THE ERROR TOLERANCES RTOL AND ATOL. C THE SIMPLEST USE IS TO TAKE THEM BOTH TO BE SCALARS. C TO OBTAIN MORE FLEXIBILITY, THEY CAN BOTH BE VECTORS. C THE CODE MUST BE TOLD YOUR CHOICE. C C **** ARE BOTH ERROR TOLERANCES RTOL, ATOL SCALARS ... C YES -- SET INFO(2) = 0 C AND INPUT SCALARS FOR BOTH RTOL AND ATOL C NO -- SET INFO(2) = 1 C AND INPUT ARRAYS FOR BOTH RTOL AND ATOL **** C C INFO(3) -- THE CODE INTEGRATES FROM T IN THE DIRECTION C OF TOUT BY STEPS. IF YOU WISH, IT WILL RETURN THE C COMPUTED SOLUTION AND DERIVATIVE AT THE NEXT C INTERMEDIATE STEP (THE INTERMEDIATE-OUTPUT MODE) OR C TOUT, WHICHEVER COMES FIRST. THIS IS A GOOD WAY TO C PROCEED IF YOU WANT TO SEE THE BEHAVIOR OF THE SOLUTION. C IF YOU MUST HAVE SOLUTIONS AT A GREAT MANY SPECIFIC C TOUT POINTS, THIS CODE WILL COMPUTE THEM EFFICIENTLY. C C **** DO YOU WANT THE SOLUTION ONLY AT C TOUT (AND NOT AT THE NEXT INTERMEDIATE STEP) ... C YES -- SET INFO(3) = 0 C NO -- SET INFO(3) = 1 **** C C INFO(4) -- TO HANDLE SOLUTIONS AT A GREAT MANY SPECIFIC C VALUES TOUT EFFICIENTLY, THIS CODE MAY INTEGRATE PAST C TOUT AND INTERPOLATE TO OBTAIN THE RESULT AT TOUT. C SOMETIMES IT IS NOT POSSIBLE TO INTEGRATE BEYOND SOME C POINT TSTOP BECAUSE THE EQUATION CHANGES THERE OR IT IS C NOT DEFINED PAST TSTOP. THEN YOU MUST TELL THE CODE C NOT TO GO PAST. C C **** CAN THE INTEGRATION BE CARRIED OUT WITHOUT ANY C RESTRICTIONS ON THE INDEPENDENT VARIABLE T ... C YES -- SET INFO(4)=0 C NO -- SET INFO(4)=1 C AND DEFINE THE STOPPING POINT TSTOP BY C SETTING RWORK(1)=TSTOP **** C C RTOL, ATOL -- YOU MUST ASSIGN RELATIVE (RTOL) AND ABSOLUTE (ATOL) C ERROR TOLERANCES TO TELL THE CODE HOW ACCURATELY YOU WANT C THE SOLUTION TO BE COMPUTED. THEY MUST BE DEFINED AS C PROGRAM VARIABLES BECAUSE THE CODE MAY CHANGE THEM. YOU C HAVE TWO CHOICES -- C BOTH RTOL AND ATOL ARE SCALARS. (INFO(2)=0) C BOTH RTOL AND ATOL ARE VECTORS. (INFO(2)=1) C IN EITHER CASE ALL COMPONENTS MUST BE NON-NEGATIVE. C C THE TOLERANCES ARE USED BY THE CODE IN A LOCAL ERROR TEST C AT EACH STEP WHICH REQUIRES ROUGHLY THAT C DABS(LOCAL ERROR) .LE. RTOL*DABS(Y)+ATOL C FOR EACH VECTOR COMPONENT. C (MORE SPECIFICALLY, A EUCLIDEAN NORM IS USED TO MEASURE C THE SIZE OF VECTORS, AND THE ERROR TEST USES THE MAGNITUDE C OF THE SOLUTION AT THE BEGINNING OF THE STEP.) C C THE TRUE (GLOBAL) ERROR IS THE DIFFERENCE BETWEEN THE TRUE C SOLUTION OF THE INITIAL VALUE PROBLEM AND THE COMPUTED C APPROXIMATION. PRACTICALLY ALL PRESENT DAY CODES, C INCLUDING THIS ONE, CONTROL THE LOCAL ERROR AT EACH STEP C AND DO NOT EVEN ATTEMPT TO CONTROL THE GLOBAL ERROR C DIRECTLY. ROUGHLY SPEAKING, THEY PRODUCE A SOLUTION Y(T) C WHICH SATISFIES THE DIFFERENTIAL EQUATIONS WITH A C RESIDUAL R(T), DY(T)/DT = DF(T,Y(T)) + R(T) , C AND, ALMOST ALWAYS, R(T) IS BOUNDED BY THE ERROR C TOLERANCES. USUALLY, BUT NOT ALWAYS, THE TRUE ACCURACY OF C THE COMPUTED Y IS COMPARABLE TO THE ERROR TOLERANCES. THIS C CODE WILL USUALLY, BUT NOT ALWAYS, DELIVER A MORE ACCURATE C SOLUTION IF YOU REDUCE THE TOLERANCES AND INTEGRATE AGAIN. C BY COMPARING TWO SUCH SOLUTIONS YOU CAN GET A FAIRLY C RELIABLE IDEA OF THE TRUE ERROR IN THE SOLUTION AT THE C BIGGER TOLERANCES. C C SETTING ATOL=0.D0 RESULTS IN A PURE RELATIVE ERROR TEST ON C THAT COMPONENT. SETTING RTOL=0. RESULTS IN A PURE ABSOLUTE C ERROR TEST ON THAT COMPONENT. A MIXED TEST WITH NON-ZERO C RTOL AND ATOL CORRESPONDS ROUGHLY TO A RELATIVE ERROR C TEST WHEN THE SOLUTION COMPONENT IS MUCH BIGGER THAN ATOL C AND TO AN ABSOLUTE ERROR TEST WHEN THE SOLUTION COMPONENT C IS SMALLER THAN THE THRESHOLD ATOL. C C PROPER SELECTION OF THE ABSOLUTE ERROR CONTROL PARAMETERS C ATOL REQUIRES YOU TO HAVE SOME IDEA OF THE SCALE OF THE C SOLUTION COMPONENTS. TO ACQUIRE THIS INFORMATION MAY MEAN C THAT YOU WILL HAVE TO SOLVE THE PROBLEM MORE THAN ONCE. IN C THE ABSENCE OF SCALE INFORMATION, YOU SHOULD ASK FOR SOME C RELATIVE ACCURACY IN ALL THE COMPONENTS (BY SETTING RTOL C VALUES NON-ZERO) AND PERHAPS IMPOSE EXTREMELY SMALL C ABSOLUTE ERROR TOLERANCES TO PROTECT AGAINST THE DANGER OF C A SOLUTION COMPONENT BECOMING ZERO. C C THE CODE WILL NOT ATTEMPT TO COMPUTE A SOLUTION AT AN C ACCURACY UNREASONABLE FOR THE MACHINE BEING USED. IT WILL C ADVISE YOU IF YOU ASK FOR TOO MUCH ACCURACY AND INFORM C YOU AS TO THE MAXIMUM ACCURACY IT BELIEVES POSSIBLE. C C RWORK(*) -- DIMENSION THIS DOUBLE PRECISION WORK ARRAY OF LENGTH C LRW IN YOUR CALLING PROGRAM. C C RWORK(1) -- IF YOU HAVE SET INFO(4)=0, YOU CAN IGNORE THIS C OPTIONAL INPUT PARAMETER. OTHERWISE YOU MUST DEFINE A C STOPPING POINT TSTOP BY SETTING RWORK(1) = TSTOP. C (FOR SOME PROBLEMS IT MAY NOT BE PERMISSIBLE TO INTEGRATE C PAST A POINT TSTOP BECAUSE A DISCONTINUITY OCCURS THERE C OR THE SOLUTION OR ITS DERIVATIVE IS NOT DEFINED BEYOND C TSTOP.) C C LRW -- SET IT TO THE DECLARED LENGTH OF THE RWORK ARRAY. C YOU MUST HAVE LRW .GE. 130+21*NEQ C C IWORK(*) -- DIMENSION THIS INTEGER WORK ARRAY OF LENGTH LIW IN C YOUR CALLING PROGRAM. C C LIW -- SET IT TO THE DECLARED LENGTH OF THE IWORK ARRAY. C YOU MUST HAVE LIW .GE. 50 C C RPAR, IPAR -- THESE ARE PARAMETER ARRAYS, OF DOUBLE PRECISION AND C INTEGER TYPE, RESPECTIVELY. YOU CAN USE THEM FOR C COMMUNICATION BETWEEN YOUR PROGRAM THAT CALLS DDEABM AND C THE DF SUBROUTINE. THEY ARE NOT USED OR ALTERED BY C DDEABM. IF YOU DO NOT NEED RPAR OR IPAR, IGNORE THESE C PARAMETERS BY TREATING THEM AS DUMMY ARGUMENTS. IF YOU DO C CHOOSE TO USE THEM, DIMENSION THEM IN YOUR CALLING PROGRAM C AND IN DF AS ARRAYS OF APPROPRIATE LENGTH. C C ********************************************************************** C * OUTPUT -- AFTER ANY RETURN FROM DDEABM * C ********************************************************************** C C THE PRINCIPAL AIM OF THE CODE IS TO RETURN A COMPUTED SOLUTION AT C TOUT, ALTHOUGH IT IS ALSO POSSIBLE TO OBTAIN INTERMEDIATE RESULTS C ALONG THE WAY. TO FIND OUT WHETHER THE CODE ACHIEVED ITS GOAL C OR IF THE INTEGRATION PROCESS WAS INTERRUPTED BEFORE THE TASK WAS C COMPLETED, YOU MUST CHECK THE IDID PARAMETER. C C C T -- THE SOLUTION WAS SUCCESSFULLY ADVANCED TO THE C OUTPUT VALUE OF T. C C Y(*) -- CONTAINS THE COMPUTED SOLUTION APPROXIMATION AT T. C YOU MAY ALSO BE INTERESTED IN THE APPROXIMATE DERIVATIVE C OF THE SOLUTION AT T. IT IS CONTAINED IN C RWORK(21),...,RWORK(20+NEQ). C C IDID -- REPORTS WHAT THE CODE DID C C *** TASK COMPLETED *** C REPORTED BY POSITIVE VALUES OF IDID C C IDID = 1 -- A STEP WAS SUCCESSFULLY TAKEN IN THE C INTERMEDIATE-OUTPUT MODE. THE CODE HAS NOT C YET REACHED TOUT. C C IDID = 2 -- THE INTEGRATION TO TOUT WAS SUCCESSFULLY C COMPLETED (T=TOUT) BY STEPPING EXACTLY TO TOUT. C C IDID = 3 -- THE INTEGRATION TO TOUT WAS SUCCESSFULLY C COMPLETED (T=TOUT) BY STEPPING PAST TOUT. C Y(*) IS OBTAINED BY INTERPOLATION. C C *** TASK INTERRUPTED *** C REPORTED BY NEGATIVE VALUES OF IDID C C IDID = -1 -- A LARGE AMOUNT OF WORK HAS BEEN EXPENDED. C (500 STEPS ATTEMPTED) C C IDID = -2 -- THE ERROR TOLERANCES ARE TOO STRINGENT. C C IDID = -3 -- THE LOCAL ERROR TEST CANNOT BE SATISFIED C BECAUSE YOU SPECIFIED A ZERO COMPONENT IN ATOL C AND THE CORRESPONDING COMPUTED SOLUTION C COMPONENT IS ZERO. THUS, A PURE RELATIVE ERROR C TEST IS IMPOSSIBLE FOR THIS COMPONENT. C C IDID = -4 -- THE PROBLEM APPEARS TO BE STIFF. C C IDID = -5,-6,-7,..,-32 -- NOT APPLICABLE FOR THIS CODE C BUT USED BY OTHER MEMBERS OF DEPAC OR POSSIBLE C FUTURE EXTENSIONS. C C *** TASK TERMINATED *** C REPORTED BY THE VALUE OF IDID=-33 C C IDID = -33 -- THE CODE HAS ENCOUNTERED TROUBLE FROM WHICH C IT CANNOT RECOVER. A MESSAGE IS PRINTED C EXPLAINING THE TROUBLE AND CONTROL IS RETURNED C TO THE CALLING PROGRAM. FOR EXAMPLE, THIS OCCURS C WHEN INVALID INPUT IS DETECTED. C C RTOL, ATOL -- THESE QUANTITIES REMAIN UNCHANGED EXCEPT WHEN C IDID = -2. IN THIS CASE, THE ERROR TOLERANCES HAVE BEEN C INCREASED BY THE CODE TO VALUES WHICH ARE ESTIMATED TO BE C APPROPRIATE FOR CONTINUING THE INTEGRATION. HOWEVER, THE C REPORTED SOLUTION AT T WAS OBTAINED USING THE INPUT VALUES C OF RTOL AND ATOL. C C RWORK, IWORK -- CONTAIN INFORMATION WHICH IS USUALLY OF NO C INTEREST TO THE USER BUT NECESSARY FOR SUBSEQUENT CALLS. C HOWEVER, YOU MAY FIND USE FOR C C RWORK(11)--WHICH CONTAINS THE STEP SIZE H TO BE C ATTEMPTED ON THE NEXT STEP. C C RWORK(12)--IF THE TOLERANCES HAVE BEEN INCREASED BY THE C CODE (IDID = -2) , THEY WERE MULTIPLIED BY THE C VALUE IN RWORK(12). C C RWORK(13)--WHICH CONTAINS THE CURRENT VALUE OF THE C INDEPENDENT VARIABLE, I.E. THE FARTHEST POINT C INTEGRATION HAS REACHED. THIS WILL BE DIFFERENT C FROM T ONLY WHEN INTERPOLATION HAS BEEN C PERFORMED (IDID=3). C C RWORK(20+I)--WHICH CONTAINS THE APPROXIMATE DERIVATIVE C OF THE SOLUTION COMPONENT Y(I). IN DDEABM, IT C IS OBTAINED BY CALLING SUBROUTINE DF TO C EVALUATE THE DIFFERENTIAL EQUATION USING T AND C Y(*) WHEN IDID=1 OR 2, AND BY INTERPOLATION C WHEN IDID=3. C C ********************************************************************** C * INPUT -- WHAT TO DO TO CONTINUE THE INTEGRATION * C * (CALLS AFTER THE FIRST) * C ********************************************************************** C C THIS CODE IS ORGANIZED SO THAT SUBSEQUENT CALLS TO CONTINUE THE C INTEGRATION INVOLVE LITTLE (IF ANY) ADDITIONAL EFFORT ON YOUR C PART. YOU MUST MONITOR THE IDID PARAMETER IN ORDER TO DETERMINE C WHAT TO DO NEXT. C C RECALLING THAT THE PRINCIPAL TASK OF THE CODE IS TO INTEGRATE C FROM T TO TOUT (THE INTERVAL MODE), USUALLY ALL YOU WILL NEED C TO DO IS SPECIFY A NEW TOUT UPON REACHING THE CURRENT TOUT. C C DO NOT ALTER ANY QUANTITY NOT SPECIFICALLY PERMITTED BELOW, C IN PARTICULAR DO NOT ALTER NEQ, T, Y(*), RWORK(*), IWORK(*) OR C THE DIFFERENTIAL EQUATION IN SUBROUTINE DF. ANY SUCH ALTERATION C CONSTITUTES A NEW PROBLEM AND MUST BE TREATED AS SUCH, I.E. C YOU MUST START AFRESH. C C YOU CANNOT CHANGE FROM VECTOR TO SCALAR ERROR CONTROL OR VICE C VERSA (INFO(2)) BUT YOU CAN CHANGE THE SIZE OF THE ENTRIES OF C RTOL, ATOL. INCREASING A TOLERANCE MAKES THE EQUATION EASIER C TO INTEGRATE. DECREASING A TOLERANCE WILL MAKE THE EQUATION C HARDER TO INTEGRATE AND SHOULD GENERALLY BE AVOIDED. C C YOU CAN SWITCH FROM THE INTERMEDIATE-OUTPUT MODE TO THE C INTERVAL MODE (INFO(3)) OR VICE VERSA AT ANY TIME. C C IF IT HAS BEEN NECESSARY TO PREVENT THE INTEGRATION FROM GOING C PAST A POINT TSTOP (INFO(4), RWORK(1)), KEEP IN MIND THAT THE C CODE WILL NOT INTEGRATE TO ANY TOUT BEYOND THE CURRENTLY C SPECIFIED TSTOP. ONCE TSTOP HAS BEEN REACHED YOU MUST CHANGE C THE VALUE OF TSTOP OR SET INFO(4)=0. YOU MAY CHANGE INFO(4) C OR TSTOP AT ANY TIME BUT YOU MUST SUPPLY THE VALUE OF TSTOP IN C RWORK(1) WHENEVER YOU SET INFO(4)=1. C C THE PARAMETER INFO(1) IS USED BY THE CODE TO INDICATE THE C BEGINNING OF A NEW PROBLEM AND TO INDICATE WHETHER INTEGRATION C IS TO BE CONTINUED. YOU MUST INPUT THE VALUE INFO(1) = 0 C WHEN STARTING A NEW PROBLEM. YOU MUST INPUT THE VALUE C INFO(1) = 1 IF YOU WISH TO CONTINUE AFTER AN INTERRUPTED TASK. C DO NOT SET INFO(1) = 0 ON A CONTINUATION CALL UNLESS YOU C WANT THE CODE TO RESTART AT THE CURRENT T. C C *** FOLLOWING A COMPLETED TASK *** C IF C IDID = 1, CALL THE CODE AGAIN TO CONTINUE THE INTEGRATION C ANOTHER STEP IN THE DIRECTION OF TOUT. C C IDID = 2 OR 3, DEFINE A NEW TOUT AND CALL THE CODE AGAIN. C TOUT MUST BE DIFFERENT FROM T. YOU CANNOT CHANGE C THE DIRECTION OF INTEGRATION WITHOUT RESTARTING. C C *** FOLLOWING AN INTERRUPTED TASK *** C TO SHOW THE CODE THAT YOU REALIZE THE TASK WAS C INTERRUPTED AND THAT YOU WANT TO CONTINUE, YOU C MUST TAKE APPROPRIATE ACTION AND RESET INFO(1) = 1 C IF C IDID = -1, THE CODE HAS ATTEMPTED 500 STEPS. C IF YOU WANT TO CONTINUE, SET INFO(1) = 1 AND C CALL THE CODE AGAIN. AN ADDITIONAL 500 STEPS C WILL BE ALLOWED. C C IDID = -2, THE ERROR TOLERANCES RTOL, ATOL HAVE BEEN C INCREASED TO VALUES THE CODE ESTIMATES APPROPRIATE C FOR CONTINUING. YOU MAY WANT TO CHANGE THEM C YOURSELF. IF YOU ARE SURE YOU WANT TO CONTINUE C WITH RELAXED ERROR TOLERANCES, SET INFO(1)=1 AND C CALL THE CODE AGAIN. C C IDID = -3, A SOLUTION COMPONENT IS ZERO AND YOU SET THE C CORRESPONDING COMPONENT OF ATOL TO ZERO. IF YOU C ARE SURE YOU WANT TO CONTINUE, YOU MUST FIRST C ALTER THE ERROR CRITERION TO USE POSITIVE VALUES C FOR THOSE COMPONENTS OF ATOL CORRESPONDING TO ZERO C SOLUTION COMPONENTS, THEN SET INFO(1)=1 AND CALL C THE CODE AGAIN. C C IDID = -4, THE PROBLEM APPEARS TO BE STIFF. IT IS VERY C INEFFICIENT TO SOLVE SUCH PROBLEMS WITH DDEABM. C THE CODE DDEBDF IN DEPAC HANDLES THIS TASK C EFFICIENTLY. IF YOU ARE ABSOLUTELY SURE YOU WANT C TO CONTINUE WITH DDEABM, SET INFO(1)=1 AND CALL C THE CODE AGAIN. C C IDID = -5,-6,-7,..,-32 --- CANNOT OCCUR WITH THIS CODE C BUT USED BY OTHER MEMBERS OF DEPAC OR POSSIBLE C FUTURE EXTENSIONS. C C *** FOLLOWING A TERMINATED TASK *** C IF C IDID = -33, YOU CANNOT CONTINUE THE SOLUTION OF THIS C PROBLEM. AN ATTEMPT TO DO SO WILL RESULT IN YOUR C RUN BEING TERMINATED. C C ********************************************************************** C***LONG DESCRIPTION C C ********************************************************************** C * DEPAC PACKAGE OVERVIEW * C ********************************************************************** C C .... YOU HAVE A CHOICE OF THREE DIFFERENTIAL EQUATION SOLVERS FROM C .... DEPAC. THE FOLLOWING BRIEF DESCRIPTIONS ARE MEANT TO AID YOU IN C .... CHOOSING THE MOST APPROPRIATE CODE FOR YOUR PROBLEM. C C .... DDERKF IS A FIFTH ORDER RUNGE-KUTTA CODE. IT IS THE SIMPLEST OF C .... THE THREE CHOICES, BOTH ALGORITHMICALLY AND IN THE USE OF THE C .... CODE. DDERKF IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND C .... MILDLY STIFF DIFFERENTIAL EQUATIONS WHEN DERIVATIVE EVALUATIONS C .... ARE NOT EXPENSIVE. IT SHOULD GENERALLY NOT BE USED TO GET HIGH C .... ACCURACY RESULTS NOR ANSWERS AT A GREAT MANY SPECIFIC POINTS. C .... BECAUSE DDERKF HAS VERY LOW OVERHEAD COSTS, IT WILL USUALLY C .... RESULT IN THE LEAST EXPENSIVE INTEGRATION WHEN SOLVING C .... PROBLEMS REQUIRING A MODEST AMOUNT OF ACCURACY AND HAVING C .... EQUATIONS THAT ARE NOT COSTLY TO EVALUATE. DDERKF ATTEMPTS TO C .... DISCOVER WHEN IT IS NOT SUITABLE FOR THE TASK POSED. C C .... DDEABM IS A VARIABLE ORDER (ONE THROUGH TWELVE) ADAMS CODE. C .... ITS COMPLEXITY LIES SOMEWHERE BETWEEN THAT OF DDERKF AND C .... DDEBDF. DDEABM IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND C .... MILDLY STIFF DIFFERENTIAL EQUATIONS WHEN DERIVATIVE EVALUATIONS C .... ARE EXPENSIVE, HIGH ACCURACY RESULTS ARE NEEDED OR ANSWERS AT C .... MANY SPECIFIC POINTS ARE REQUIRED. DDEABM ATTEMPTS TO DISCOVER C .... WHEN IT IS NOT SUITABLE FOR THE TASK POSED. C C .... DDEBDF IS A VARIABLE ORDER (ONE THROUGH FIVE) BACKWARD C .... DIFFERENTIATION FORMULA CODE. IT IS THE MOST COMPLICATED OF C .... THE THREE CHOICES. DDEBDF IS PRIMARILY DESIGNED TO SOLVE STIFF C .... DIFFERENTIAL EQUATIONS AT CRUDE TO MODERATE TOLERANCES. C .... IF THE PROBLEM IS VERY STIFF AT ALL, DDERKF AND DDEABM WILL BE C .... QUITE INEFFICIENT COMPARED TO DDEBDF. HOWEVER, DDEBDF WILL BE C .... INEFFICIENT COMPARED TO DDERKF AND DDEABM ON NON-STIFF PROBLEMS C .... BECAUSE IT USES MUCH MORE STORAGE, HAS A MUCH LARGER OVERHEAD, C .... AND THE LOW ORDER FORMULAS WILL NOT GIVE HIGH ACCURACIES C .... EFFICIENTLY. C C .... THE CONCEPT OF STIFFNESS CANNOT BE DESCRIBED IN A FEW WORDS. C .... IF YOU DO NOT KNOW THE PROBLEM TO BE STIFF, TRY EITHER DDERKF C .... OR DDEABM. BOTH OF THESE CODES WILL INFORM YOU OF STIFFNESS C .... WHEN THE COST OF SOLVING SUCH PROBLEMS BECOMES IMPORTANT. C C ********************************************************************* C***REFERENCES SHAMPINE, L. F., WATTS, H. A., *DEPAC - DESIGN OF A USER C ORIENTED PACKAGE OF ODE SOLVERS*, SAND79-2374, C SANDIA LABORATORIES, 1979. C***ROUTINES CALLED DDES,XERRWV C***END PROLOGUE DDEABM C