C ALGORITHM 540R (REMARK ON ALG.540), COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 343-344. SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF, * INDEX,WORK,IWORK) C C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C THIS IS THE MARCH 24, 1978 VERSION OF PDECOL. C C THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY ANSI-FORTRAN C RULES AS WAS CONVENIENTLY POSSIBLE. THE FORTRAN USED VIOLATES ANSI C STANDARDS IN THE TWO WAYS LISTED BELOW.... C C 1. SUBSCRIPTS OF THE GENERAL FORM C*V1 + V2 + V3 ARE USED C (POSSIBLY IN A PERMUTED ORDER), WHERE C IS AN INTEGER CONSTANT C AND V1, V2, AND V3 ARE INTEGER VARIABLES. C C 2. ARRAY NAMES APPEAR SINGLY IN DATA STATEMENTS IN THE ROUTINES C BSPLVN AND COSET. C C MACHINE DEPENDENT FEATURES...... C C THIS VERSION OF PDECOL WAS DESIGNED FOR USE ON CDC MACHINES WITH C A WORD LENGTH OF 60 BITS. WE DO NOT RECOMMEND THE USE OF PDECOL WITH C WORD LENGTHS OF LESS THAN 48 BITS. THE MOST IMPORTANT MACHINE C AND WORD LENGTH DEPENDENT CONSTANTS ARE DEFINED IN THE BLOCK DATA C AND IN SUBROUTINES COLPNT AND COSET. THE USER SHOULD CHECK THESE C CAREFULLY FOR APPROPRIATENESS FOR HIS LOCAL SITUATION. THE FORTRAN C FUNCTIONS USED BY EACH ROUTINE ARE LISTED IN THE COMMENTS TO C FACILITATE CONVERSION TO DOUBLE PRECISION. C C----------------------------------------------------------------------- C----------------------------------------------------------------------- C C----------------------------------------------------------------------- C PDECOL IS THE DRIVER ROUTINE FOR A SOPHISTICATED PACKAGE OF C SUBROUTINES WHICH IS DESIGNED TO SOLVE THE GENERAL SYSTEM OF C NPDE NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS OF AT MOST SECOND C ORDER ON THE INTERVAL (XLEFT,XRIGHT) FOR T .GT. T0 WHICH IS OF THE C FORM.... C C DU/DT = F( T, X, U, UX, UXX ) C C WHERE C C U = ( U(1), U(2), ... , U(NPDE) ) C UX = ( UX(1), UX(2), ... , UX(NPDE) ) C UXX = (UXX(1),UXX(2), ... ,UXX(NPDE) ) . C C EACH U(K) IS A FUNCTION OF THE SCALAR QUANTITIES T AND X. C UX(K) REPRESENTS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT C TO THE VARIABLE X, UXX(K) REPRESENTS THE SECOND PARTIAL DERIVATIVE C OF U(K) WITH RESPECT TO THE VARIABLE X, AND DU/DT IS THE VECTOR OF C PARTIAL DERIVATIVES OF U WITH RESPECT TO THE TIME VARIABLE T. C F REPRESENTS AN ARBITRARY VECTOR VALUED FUNCTION WHOSE NPDE C COMPONENTS DEFINE THE RESPECTIVE PARTIAL DIFFERENTIAL EQUATIONS OF C THE PDE SYSTEM. SEE SUBROUTINE F DESCRIPTION BELOW. C C BOUNDARY CONDITIONS C C DEPENDING ON THE TYPE OF PDE(S), 0, 1, OR 2 BOUNDARY CONDITIONS C ARE REQUIRED FOR EACH PDE IN THE SYSTEM. THESE ARE IMPOSED AT XLEFT C AND/OR XRIGHT AND EACH MUST BE OF THE FORM.... C C B(U,UX) = Z(T) C C WHERE B AND Z ARE ARBITRARY VECTOR VALUED FUNCTIONS WITH C NPDE COMPONENTS AND U, UX, AND T ARE AS ABOVE. THESE BOUNDARY C CONDITIONS MUST BE CONSISTENT WITH THE INITIAL CONDITIONS WHICH ARE C DESCRIBED NEXT. C C INITIAL CONDITIONS C C EACH SOLUTION COMPONENT U(K) IS ASSUMED TO BE A KNOWN (USER C PROVIDED) FUNCTION OF X AT THE INITIAL TIME T = T0. THE C INITIAL CONDITION FUNCTIONS MUST BE CONSISTENT WITH THE BOUNDARY C CONDITIONS ABOVE, I.E. THE INITIAL CONDITION FUNCTIONS MUST C SATISFY THE BOUNDARY CONDITIONS FOR T = T0. SEE SUBROUTINE UINIT C DESCRIPTION BELOW. C----------------------------------------------------------------------- C C REQUIRED USER SUPPLIED SUBROUTINES C C THE USER IS REQUIRED TO CONSTRUCT THREE SUBPROGRAMS AND A MAIN C PROGRAM WHICH DEFINE THE PDE PROBLEM WHOSE SOLUTION IS TO BE C ATTEMPTED. THE THREE SUBPROGRAMS ARE... C C 1) SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE ) C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) C THIS ROUTINE DEFINES THE DESIRED PARTIAL DIFFERENTIAL C EQUATIONS TO BE SOLVED. THE PACKAGE PROVIDES VALUES OF THE C INPUT SCALARS T AND X AND INPUT ARRAYS (LENGTH NPDE) U, UX, C AND UXX, AND THE USER MUST CONSTRUCT THIS ROUTINE TO COMPUTE C THE OUTPUT ARRAY FVAL (LENGTH NPDE) WHICH CONTAINS THE C CORRESPONDING VALUES OF THE RIGHT HAND SIDES OF THE DESIRED C PARTIAL DIFFERENTIAL EQUATIONS, I.E. C C FVAL(K) = THE VALUE OF THE RIGHT HAND SIDE OF THE K-TH PDE IN C THE PDE SYSTEM ABOVE, FOR K = 1 TO NPDE. C C THE INCOMING VALUE OF THE SCALAR QUANTITY X WILL BE A C COLLOCATION POINT VALUE (SEE INITAL AND COLPNT) AND THE C INCOMING VALUES IN THE ARRAYS U, UX AND UXX CORRESPOND TO THIS C POINT X AND TIME T. C RETURN C END C C 2) SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE ) C DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) C DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH NEEDED C INFORMATION ABOUT THE BOUNDARY CONDITION FUNCTIONS B AND Z C ABOVE. THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES C T, X, U, AND UX, AND THE USER IS TO DEFINE THE CORRESPONDING C OUTPUT VALUES OF THE DERIVATIVES OF THE FUNCTIONS B AND Z C WHERE.... C DBDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO C THE J-TH VARIABLE U(J). C DBDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO C THE J-TH VARIABLE UX(J). C DZDT(K) = DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR C FUNCTION Z(T) ABOVE WITH RESPECT TO THE C VARIABLE T. C NOTE... THE INCOMING VALUE OF X WILL BE EITHER XLEFT OR XRIGHT. C IF NO BOUNDARY CONDITION IS DESIRED FOR SAY THE K-TH PDE AT C ONE OR BOTH OF THE ENDPOINTS XLEFT OR XRIGHT, THEN DBDU(K,K) C AND DBDUX(K,K) SHOULD BOTH BE SET TO ZERO WHEN BNDRY IS C CALLED FOR THAT POINT. WE REFER TO THIS AS A NULL BOUNDARY C CONDITION. THIS ROUTINE CAN BE STRUCTURED AS FOLLOWS... C THE COMMON BLOCK /ENDPT/ IS NOT A PART OF PDECOL AND C MUST BE SUPPLIED AND DEFINED BY THE USER. C COMMON /ENDPT/ XLEFT C IF( X .NE. XLEFT ) GO TO 10 C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE LEFT BOUNDARY POINT C X = XLEFT. C RETURN C 10 CONTINUE C HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), C AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE RIGHT BOUNDARY POINT C X = XRIGHT. C RETURN C END C C 3) SUBROUTINE UINIT( X, U, NPDE ) C DIMENSION U(NPDE) C THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH THE C NEEDED INITIAL CONDITION FUNCTION VALUES. THE PACKAGE C PROVIDES A VALUE OF THE INPUT VARIABLE X, AND THE USER IS TO C DEFINE THE PROPER INITIAL VALUES (AT T = T0) FOR ALL OF THE C PDE COMPONENTS, I.E. C U(K) = DESIRED INITIAL VALUE OF PDE COMPONENT U(K) AT C X AND T = T0 FOR K = 1 TO NPDE. C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT C VALUE. THE INITIAL CONDITIONS AND BOUNDARY CONDITIONS C MUST BE CONSISTENT (SEE ABOVE). C RETURN C END C----------------------------------------------------------------------- C C OPTIONAL USER SUPPLIED SUBROUTINE C C IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE C ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST C PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE C DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT C JACOBIAN MATRIX GENERATION. ON MOST COMPUTER SYSTEMS, THE USER WILL C BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE C OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW). C C 1) SUBROUTINE DERIVF( T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE ) C DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) C DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) C THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX, C AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE C THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS C DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE... C DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE U(J). C DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE UX(J). C DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE UXX(J). C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT C VALUE. C RETURN C END C----------------------------------------------------------------------- C C METHODS USED C C THE PACKAGE PDECOL IS BASED ON THE METHOD OF LINES AND USES A C FINITE ELEMENT COLLOCATION PROCEDURE (WITH PIECEWISE POLYNOMIALS C AS THE TRIAL SPACE) FOR THE DISCRETIZATION OF THE SPATIAL VARIABLE C X. THE COLLOCATION PROCEDURE REDUCES THE PDE SYSTEM TO A SEMI- C DISCRETE SYSTEM WHICH THEN DEPENDS ONLY ON THE TIME VARIABLE T. C THE TIME INTEGRATION IS THEN ACCOMPLISHED BY USE OF SLIGHTLY C MODIFIED STANDARD TECHNIQUES (SEE REFS. 1,2). C C PIECEWISE POLYNOMIALS C C THE USER IS REQUIRED TO SELECT THE PIECEWISE POLYNOMIAL SPACE C WHICH IS TO BE USED TO COMPUTE HIS APPROXIMATE SOLUTION. FIRST, THE C ORDER, KORD, OF THE POLYNOMIALS TO BE USED MUST BE SPECIFIED C (KORD = POLYNOMIAL DEGREE + 1). NEXT, THE NUMBER OF PIECES C (INTERVALS), NINT, INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS C TO BE DIVIDED, IS CHOSEN. THE NINT + 1 DISTINCT BREAKPOINTS OF C THE DOMAIN MUST BE DEFINED AND SET INTO THE ARRAY XBKPT IN C STRICTLY INCREASING ORDER, I.E. C XLEFT=XBKPT(1) .LT. XBKPT(2) .LT. ... .LT. XBKPT(NINT+1)=XRIGHT. C THE APPROXIMATE SOLUTION AT ANY TIME T WILL BE A POLYNOMIAL OF C ORDER KORD OVER EACH SUBINTERVAL (XBKPT(I),XBKPT(I+1)). THE C NUMBER OF CONTINUITY CONDITIONS, NCC, TO BE IMPOSED ACROSS ALL OF C THE BREAKPOINTS IS THE LAST PIECE OF USER SUPPLIED DATA WHICH IS C REQUIRED TO UNIQUELY DETERMINE THE DESIRED PIECEWISE POLYNOMIAL C SPACE. FOR EXAMPLE, NCC = 2 WOULD REQUIRE THAT THE APPROXIMATE C SOLUTION (MADE UP OF THE SEPARATE POLYNOMIAL PIECES) AND ITS FIRST C SPATIAL DERIVATIVE BE CONTINUOUS AT THE BREAKPOINTS AND HENCE ON C THE ENTIRE DOMAIN (XLEFT,XRIGHT). NCC = 3 WOULD REQUIRE THAT THE C APPROXIMATE SOLUTION AND ITS FIRST AND SECOND SPATIAL DERIVATIVES C BE CONTINUOUS AT THE BREAKPOINTS, ETC. THE DIMENSION OF THIS LINEAR C SPACE IS KNOWN AND FINITE AND IS NCPTS = KORD*NINT - NCC*(NINT-1). C THE WELL-KNOWN B-SPLINE BASIS (SEE REF. 3) FOR THIS SPACE IS USED C BY PDECOL AND IT CONSISTS OF NCPTS KNOWN PIECEWISE POLYNOMIAL C FUNCTIONS BF(I,X), FOR I=1 TO NCPTS, WHICH DO NOT DEPEND ON THE C TIME VARIABLE T. WE WISH TO EMPHASIZE THAT THE PIECEWISE POLYNOMIAL C SPACE USED IN PDECOL (WHICH IS SELECTED BY THE USER) WILL DETERMINE C THE MAGNITUDE OF THE SPATIAL DISCRETIZATION ERRORS IN THE COMPUTED C APPROXIMATE SOLUTION. THE PACKAGE HAS NO CONTROL OVER ERRORS C INTRODUCED BY THE USERS CHOICE OF THIS SPACE. SEE INPUT PARAMETERS C BELOW. C C COLLOCATION OVER PIECEWISE POLYNOMIALS C C THE BASIC ASSUMPTION MADE IS THAT THE APPROXIMATE SOLUTION C SATISFIES C NCPTS C U(T,X) = SUM C(I,T) * BF(I,X) C I=1 C C WHERE THE UNKNOWN COEFFICIENTS C DEPEND ONLY ON THE TIME T AND C THE KNOWN BASIS FUNCTIONS DEPEND ONLY ON X (WE HAVE ASSUMED THAT C NPDE = 1 FOR CONVENIENCE). SO, AT ANY GIVEN TIME T THE APPROX- C IMATE SOLUTION IS A PIECEWISE POLYNOMIAL IN THE USER CHOSEN SPACE. C THE SEMI-DISCRETE EQUATIONS (ACTUALLY ORDINARY DIFFERENTIAL C EQUATIONS) WHICH DETERMINE THE COEFFICIENTS C ARE OBTAINED BY C REQUIRING THAT THE ABOVE APPROXIMATE U(T,X) SATISFY THE PDE AND C BOUNDARY CONDITIONS EXACTLY AT A SET OF NCPTS COLLOCATION POINTS C (SEE COLPNT). THUS, PDECOL ACTUALLY COMPUTES THE BASIS FUNCTION C COEFFICIENTS RATHER THAN SPECIFIC APPROXIMATE SOLUTION VALUES. C C REFERENCES C C 1. MADSEN, N.K. AND R.F. SINCOVEC, PDECOL - COLLOCATION SOFTWARE C FOR PARTIAL DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. , NO. , C C 2. SINCOVEC, R.F. AND N.K. MADSEN, SOFTWARE FOR NONLINEAR PARTIAL C DIFFERENTIAL EQUATIONS, ACM-TOMS, VOL. 1, NO. 3, C SEPTEMBER 1975, PP. 232-260. C 3. HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY C 1976. C 4. DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J. C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472. C----------------------------------------------------------------------- C C USE OF PDECOL C C PDECOL IS CALLED ONCE FOR EACH DESIRED OUTPUT VALUE (TOUT) OF THE C TIME T, AND IT IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, C STIFIB, WHICH ADVANCES THE TIME BY TAKING SINGLE STEPS UNTIL C T .GE. TOUT. INTERPOLATION TO THE EXACT TIME TOUT IS THEN DONE. C SEE TOUT BELOW. C C C SUMMARY OF SUGGESTED INPUT VALUES C C IT IS OF COURSE IMPOSSIBLE TO SUGGEST INPUT PARAMETER VALUES WHICH C ARE APPROPRIATE FOR ALL PROBLEMS. THE FOLLOWING SUGGESTIONS ARE TO C BE USED ONLY IF YOU HAVE NO IDEA OF BETTER VALUES FOR YOUR PROBLEM. C C DT = 1.E-10 C XBKPT = CHOOSE NINT+1 EQUALLY SPACED VALUES SUCH THAT XBKPT(1) = C XLEFT AND XBKPT(NINT+1) = XRIGHT. C EPS = 1.E-4 C NINT = ENOUGH SO THAT ANY FINE STRUCTURE OF THE PROBLEM MAY BE C RESOLVED. C KORD = 4 C NCC = 2 C MF = 22 C INDEX = 1 (ON FIRST CALL ONLY, THEN 0 THEREAFTER). C C C THE INPUT PARAMETERS ARE.. C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE C (USED ONLY ON FIRST CALL). C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. SINCE C THE PACKAGE CHOOSES ITS OWN TIME STEP SIZES, THE C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT. C DT = THE INITIAL STEP SIZE IN T, IF INDEX = 1, OR, THE C MAXIMUM STEP SIZE ALLOWED (MUST BE .GT. 0), IF INDEX = 3. C USED FOR INPUT ONLY WHEN INDEX = 1 OR 3. SEE BELOW. C XBKPT = THE ARRAY OF PIECEWISE POLYNOMIAL BREAKPOINTS. C THE NINT+1 VALUES MUST BE STRICTLY INCREASING WITH C XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT (USED ONLY C ON FIRST CALL). C EPS = THE RELATIVE TIME ERROR BOUND (USED ONLY ON THE C FIRST CALL, UNLESS INDEX = 4). SINGLE STEP ERROR C ESTIMATES DIVIDED BY CMAX(I) WILL BE KEPT LESS THAN C EPS IN ROOT-MEAN-SQUARE NORM. THE VECTOR CMAX OF WEIGHTS C IS COMPUTED IN PDECOL. INITIALLY CMAX(I) IS SET TO C ABS(C(I)), WITH A DEFAULT VALUE OF 1 IF ABS(C(I)) .LT. 1. C THEREAFTER, CMAX(I) IS THE LARGEST VALUE C OF ABS(C(I)) SEEN SO FAR, OR THE INITIAL CMAX(I) IF C THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE C APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT C STATEMENTS 50 AND 130 BELOW. THE USER SHOULD EXERCISE C SOME DISCRETION IN CHOOSING EPS. IN GENERAL, THE C OVERALL RUNNING TIME FOR A PROBLEM WILL BE GREATER IF C EPS IS CHOSEN SMALLER. THERE IS USUALLY LITTLE REASON TO C CHOOSE EPS MUCH SMALLER THAN THE ERRORS WHICH ARE BEING C INTRODUCED BY THE USERS CHOICE OF THE POLYNOMIAL SPACE. C SEE RELATED COMMENTS CONCERNING CMAX BELOW STATEMENT 40. C NINT = THE NUMBER OF SUBINTERVALS INTO WHICH THE SPATIAL DOMAIN C (XLEFT,XRIGHT) IS TO BE DIVIDED (MUST BE .GE. 1) C (USED ONLY ON FIRST CALL). C KORD = THE ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED. C ITS VALUE MUST BE GREATER THAN 2 AND LESS THAN 21. FOR C FIRST ATTEMPTS WE RECOMMEND KORD = 4. IF THE SOLUTION C IS SMOOTH AND MUCH ACCURACY IS DESIRED, HIGHER VALUES C MAY PROVE TO BE MORE EFFICIENT. WE HAVE SELDOM USED C VALUES OF KORD IN EXCESS OF 8 OR 9, THOUGH THEY ARE C AVAILABLE FOR USE IN PDECOL (USED ONLY ON FIRST CALL). C NCC = THE NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED ON THE C APPROXIMATE SOLUTION AT THE BREAKPOINTS IN XBKPT. C NCC MUST BE GREATER THAN 1 AND LESS THAN KORD. WE C RECOMMEND THE USE OF NCC = 2 C SINCE THEORY PREDICTS THAT DRAMATICALLY MORE C ACCURATE RESULTS CAN OFTEN BE OBTAINED USING THIS CHOICE C (USED ONLY ON FIRST CALL). C NPDE = THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS IN THE SYSTEM C TO BE SOLVED (USED ONLY ON FIRST CALL). C MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS C INDEX = 4). ALLOWED VALUES ARE 11, 12, 21, 22. C FOR FIRST ATTEMPTS WE RECOMMEND THE USE OF MF = 22. C MF HAS TWO DECIMAL DIGITS, METH AND MITER C (MF = 10*METH + MITER). C METH IS THE BASIC METHOD INDICATOR.. C METH = 1 MEANS THE ADAMS METHODS (GENERALIZATIONS OF C CRANK-NICOLSON). C METH = 2 MEANS THE BACKWARD DIFFERENTIATION C FORMULAS (BDF), OR STIFF METHODS OF GEAR. C MITER IS THE ITERATION METHOD INDICATOR C AND DETERMINES HOW THE JACOBIAN MATRIX IS C TO BE COMPUTED.. C MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN. C FOR THIS USER SUPPLIES SUBROUTINE DERIVF. C SEE DESCRIPTION ABOVE. C MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED C INTERNALLY BY FINITE DIFFERENCES. SEE C SUBROUTINES PSETIB AND DIFFF. C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 1 THIS IS THE FIRST CALL FOR THIS PROBLEM. C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM, C AND INTEGRATION IS TO CONTINUE. C 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT C EXACTLY (NO INTERPOLATION IS DONE). SEE NOTE C BELOW. ASSUMES TOUT .GE. THE CURRENT T. C IF TOUT IS .LT. THE CURRENT TIME, THEN TOUT IS C RESET TO THE CURRENT TIME AND CONTROL IS C RETURNED TO THE USER. A CALL TO VALUES WILL C PRODUCE ANSWERS FOR THE NEW VALUE OF TOUT. C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING C PROGRAM AFTER ONE STEP. TOUT IS IGNORED AND C DT MUST BE SET .GT. 0 TO A MAXIMUM ALLOWED C DT VALUE. SEE ABOVE. C 4 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, C AND THE USER HAS RESET EPS AND/OR MF. C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, C IT NEED NOT BE RESET FOR NORMAL CONTINUATION. C C NOTE.. THE PACKAGE MUST HAVE TAKEN AT LEAST ONE SUCCESSFUL TIME C STEP BEFORE A CALL WITH INDEX = 2 OR 4 IS ALLOWED. C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. C A CHANGE OF PARAMETERS WITH INDEX = 4 CAN BE MADE AFTER C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN PROVIDED AT LEAST C ONE SUCCESSFUL TIME STEP HAS BEEN MADE. C C WORK = FLOATING POINT WORKING ARRAY FOR PDECOL. WE RECOMMEND C THAT IT BE INITIALIZED TO ZERO BEFORE THE FIRST CALL C TO PDECOL. ITS TOTAL LENGTH MUST BE AT LEAST C C KORD + 4*NPDE + 9*NPDE**2 + NCPTS*(3*KORD + 2) + C NPDE*NCPTS*(3*ML + MAXDER + 7) C C WHERE ML AND MAXDER ARE DEFINED BELOW (SEE STORAGE C ALLOCATION). C IWORK = INTEGER WORKING ARRAY FOR PDECOL. THE FIRST TWO C LOCATIONS MUST BE DEFINED AS FOLLOWS... C IWORK(1) = LENGTH OF USERS ARRAY WORK C IWORK(2) = LENGTH OF USERS ARRAY IWORK C THE TOTAL LENGTH OF IWORK MUST BE AT LEAST C NCPTS*(NPDE + 1). C OUTPUT C C THE SOLUTION VALUES ARE NOT RETURNED DIRECTLY TO THE USER BY PDECOL. C THE METHODS USED IN PDECOL COMPUTE BASIS FUNCTION COEFFICIENTS, SO C THE USER (AFTER A RETURN FROM PDECOL) MUST CALL THE PACKAGE ROUTINE C VALUES TO OBTAIN HIS APPROXIMATE SOLUTION VALUES AT ANY DESIRED SPACE C POINTS X AT THE TIME T = TOUT. SEE THE COMMENTS IN SUBROUTINE VALUES C FOR DETAILS ON HOW TO PROPERLY MAKE THE CALL. C C EXECUTION ERROR MESSAGES WILL BE PRINTED BY PDECOL ON LOGICAL UNIT C LOUT WHICH IS THE ONLY VARIABLE IN THE COMMON BLOCK /IOUNIT/. A C DEFAULT OF LOUT = 3 IS SET IN THE BLOCK DATA. C C THE COMMON BLOCK /GEAR0/ CONTAINS THE VARIABLES DTUSED, NQUSED, C NSTEP, NFE, AND NJE AND CAN BE ACCESSED EXTERNALLY BY THE USER IF C DESIRED. RESPECTIVELY, IT CONTAINS THE STEP SIZE LAST USED (SUCCESS- C FULLY), THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN C SO FAR, THE NUMBER OF RESIDUAL EVALUATIONS (RES CALLS) SO FAR, C AND THE NUMBER OF MATRIX EVALUATIONS (PSETIB CALLS) SO FAR. C DIFFUN CALLS ARE COUNTED IN WITH RESIDUAL EVALUATIONS. C C THE OUTPUT PARAMETERS ARE.. C DT = THE STEP SIZE USED LAST, WHETHER SUCCESSFULLY OR NOT. C TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL, C AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS C UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT C IS THE CURRENT VALUE OF T TO WHICH THE INTEGRATION C HAS BEEN COMPLETED. C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, C WITH THE FOLLOWING VALUES AND MEANINGS.. C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE C ERROR TEST EVEN AFTER REDUCING DT BY A FACTOR OF C 1.E10 FROM ITS INITIAL VALUE. C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE C CORRECTOR CONVERGENCE EVEN AFTER REDUCING DT BY A C FACTOR OF 1.E10 FROM ITS INITIAL VALUE. C -4 SINGULAR MATRIX ENCOUNTERED. PROBABLY DUE TO STORAGE C OVERWRITES. C -5 INDEX WAS 4 ON INPUT, BUT THE DESIRED CHANGES OF C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT C WAS NOT BEYOND T. INTERPOLATION TO T = TOUT WAS C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, C SIMPLY CALL AGAIN WITH INDEX = 4 AND A NEW TOUT. C -6 ILLEGAL INDEX VALUE. C -7 ILLEGAL EPS VALUE. C -8 AN ATTEMPT TO INTEGRATE IN THE WRONG DIRECTION. THE C SIGN OF DT IS WRONG RELATIVE TO T0 AND TOUT. C -9 DT .EQ. 0.0. C -10 ILLEGAL NINT VALUE. C -11 ILLEGAL KORD VALUE. C -12 ILLEGAL NCC VALUE. C -13 ILLEGAL NPDE VALUE. C -14 ILLEGAL MF VALUE. C -15 ILLEGAL BREAKPOINTS - NOT STRICTLY INCREASING. C -16 INSUFFICIENT STORAGE FOR WORK OR IWORK. C----------------------------------------------------------------------- C C C SUMMARY OF ALL PACKAGE ROUTINES C C PDECOL - STORAGE ALLOCATION, ERROR CHECKING, INITIALIZATION, REPEATED C CALLS TO STIFIB TO ADVANCE THE TIME. C C INTERP - INTERPOLATES COMPUTED BASIS FUNCTION COEFFICIENTS TO THE C DESIRED OUTPUT TIMES, TOUT, FOR USE BY VALUES. C C INITAL - INITIALIZATION, GENERATION AND STORAGE OF PIECEWISE POLY- C NOMIAL SPACE BASIS FUNCTION VALUES AND DERIVATIVES, DET- C ERMINES THE BASIS FUNCTION COEFFICINTS OF THE PIECEWISE C POLYNOMIALS WHICH INTERPOLATE THE USERS INITIAL CONDITIONS. C C COLPNT - GENERATION OF REQUIRED COLLOCATION POINTS. C C BSPLVD - B-SPLINE PACKAGE ROUTINES WHICH ALLOW FOR EVALUATION OF C BSPLVN ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE. C INTERV C C VALUES - GENERATION AT ANY POINT(S) OF VALUES OF THE COMPUTED C APPROXIMATE SOLUTION AND ITS DERIVATIVES WHICH ARE C PIECEWISE POLYNOMIALS. THE SUBROUTINE IS CALLED ONLY BY C THE USER. C C STIFIB - CORE INTEGRATOR, TAKES SINGLE TIME STEPS TO ADVANCE THE C TIME. ASSEMBLES AND SOLVES THE PROPER NONLINEAR EQUATIONS C WHICH ARE RELATED TO USE OF ADAMS OR GEAR TYPE INTEGRATION C FORMULAS. CHOOSES PROPER STEP SIZE AND INTEGRATION FORMULA C ORDER TO MAINTAIN A DESIRED ACCURACY. DESIGNED FOR ODE C PROBLEMS OF THE FORM A * (DY/DT) = G(T,Y). C C COSET - GENERATES INTEGRATION FORMULA AND ERROR CONTROL COEFFICIENTS. C C RES - COMPUTES RESIDUAL VECTORS USED IN SOLVING THE NONLINEAR C EQUATIONS BY A MODIFIED NEWTON METHOD. C C DIFFUN - COMPUTES A**-1 * G(T,Y) WHERE A AND G ARE AS ABOVE (STIFIB). C C ADDA - ADDS THE A MATRIX TO A GIVEN MATRIX IN BAND FORM. C C EVAL - EVALUATES THE COMPUTED PIECEWISE POLYNOMIAL SOLUTION AND C DERIVATIVES AT COLLOCATION POINTS. C C GFUN - EVALUATES THE FUNCTION G(T,Y) BY CALLING EVAL AND THE USER C SUBROUTINES F AND BNDRY. C C PSETIB - GENERATES PROPER JACOBIAN MATRICES REQUIRED BY THE MODIFIED C NEWTON METHOD. C C DIFFF - PERFORMS SAME ROLE AS THE USER ROUTINE DERIVF. COMPUTES C DERIVATIVE APPROXIMATIONS BY USE OF FINITE DIFFERENCES. C C DECB - PERFORM AN LU DECOMPOSTION AND FORWARD AND BACKWARD C SOLB SUBSTITUTION FOR SOLVING BANDED SYSTEMS OF LINEAR EQUATIONS. C C----------------------------------------------------------------------- C C C STORAGE ALLOCATION C C SINCE PDECOL IS A DYNAMICALLY DIMENSIONED PROGRAM, MOST OF ITS C WORKING STORAGE IS PROVIDED BY THE USER IN THE ARRAYS WORK AND IWORK. C THE FOLLOWING GIVES A LIST OF THE ARRAYS WHICH MAKE UP THE CONTENTS C WORK AND IWORK, THEIR LENGTHS, AND THEIR USES. WHEN MORE THAN ONE C NAME IS GIVEN, IT INDICATES THAT DIFFERENT NAMES ARE USED FOR THE C SAME ARRAY IN DIFFERENT PARTS OF THE PROGRAM. THE DIFFERENT NAMES C OCCUR BECAUSE PDECOL IS AN AMALGAMATION OF SEVERAL OTHER CODES C WRITTEN BY DIFFERENT PEOPLE AND WE HAVE TRIED TO LEAVE THE SEPARATE C PARTS AS UNCHANGED FROM THEIR ORIGINAL VERSIONS AS POSSIBLE. C C C NAMES LENGTH USE C --------- ------------ ------------------------------------- C C BC 4*NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK C C A 3*KORD*NCPTS BASIS FUNCTION VALUES AT COLLOCATION POINT C WORK(IW1) C C XT NCPTS + KORD BREAKPOINT SEQUENCE FOR GENERATION OF BASI C WORK(IW2) FUNCTION VALUES. C C XC NCPTS COLLOCATION POINTS. C WORK(IW3) C C CMAX NPDE*NCPTS VALUES USED IN ESTIMATING TIME C YMAX INTEGRATION ERRORS. C WORK(IW4) C C ERROR NPDE*NCPTS TIME INTEGRATION ERRORS. C WORK(IW5) C C SAVE1 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW6) METHOD. C C SAVE2 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW7) METHOD. C C SAVE3 NPDE*NCPTS WORKING STORAGE FOR THE TIME INTEGRATION C WORK(IW8) METHOD. C C UVAL 3*NPDE WORKING STORAGE FOR VALUES OF U, UX, AND C WORK(IW9) UXX AT ONE POINT. C C C NPDE*NCPTS* CURRENT BASIS FUNCTION COEFFICIENT VALUES C Y (MAXDER+1) AND THEIR SCALED TIME DERIVATIVES. C WORK(IW10) C C DFDU NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW11) JACOBIAN MATRIX. C C DFDUX NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW12) JACOBIAN MATRIX. C C DFDUXX NPDE**2 WORKING STORAGE USED TO COMPUTE THE C WORK(IW13) JACOBIAN MATRIX. C C DBDU NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK(IW14) C C DBDUX NPDE**2 BOUNDARY CONDITION INFORMATION. C WORK(IW15) C C DZDT NPDE BOUNDARY CONDITION INFORMATION. C WORK(IW16) C C PW NPDE*NCPTS* STORAGE AND PROCESSING OF THE JACOBIAN C WORK(IW17) (3*ML+1) MATRIX. C C ILEFT NCPTS POINTERS TO BREAKPOINT SEQUENCE FOR C IWORK GENERATION OF BASIS FUNCTION VALUES. C C IPIV NPDE*NCPTS PIVOT INFORMATION FOR THE LU DECOMPOSED C IWORK(IW18) JACOBIAN MATRIX PW. C C WHERE... C C NCPTS = KORD*NINT - NCC*(NINT-1) C ML = NPDE*(KORD+IQUAD-1) - 1 C IQUAD = 1 IF KORD = 3 AND A NULL BOUNDARY CONDITION EXISTS C IQUAD = 0 OTHERWISE C MAXDER = 5 UNLESS OTHERWISE SET BY THE USER INTO /OPTION/. C C THE COMMON BLOCK /OPTION/ CONTAINS THE TWO VARIABLES NOGAUS AND C MAXDER. NOGAUS IS SET .EQ. 0 IN THE BLOCK DATA. IT CAN BE CHANGED C TO BE SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE NOT C DESIRED WHEN NCC = 2 (SEE ABOVE AND COLPNT). MAXDER IS SET C .EQ. 5 IN THE BLOCK DATA AND ITS VALUE REPRESENTS THE C MAXIMUM ORDER OF TIME INTEGRATION FORMULA ALLOWED. ITS VALUE C AFFECTS THE STORAGE REQUIRED IN WORK AND MAY BE CHANGED IF C DESIRED. SEE COSET FOR RESTRICTIONS. THESE CHANGES MAY BE MADE BY C THE USER BY ACCESSING /OPTION/ IN HIS CALLING PROGRAM (BEFORE THE C FIRST CALL TO PDECOL) OR BY CHANGING THE DATA STATEMENT IN C THE BLOCK DATA. C C----------------------------------------------------------------------- C C C COMMUNICATION C C EACH SUBROUTINE IN THE PACKAGE CONTAINS A COMMUNICATION SUMMARY C AS INDICATED BELOW. C C PACKAGE ROUTINES CALLED.. EVAL,INITAL,INTERP,STIFIB C USER ROUTINES CALLED.. BNDRY C CALLED BY.. USERS MAIN PROGRAM C FORTRAN FUNCTIONS USED.. ABS,AMAX1,FLOAT,SQRT C----------------------------------------------------------------------- SAVE COMMON /GEAR0/ DTUSED,NQUSED,NSTEP,NFE,NJE COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W COMMON /OPTION/ NOGAUS,MAXDER COMMON /SIZES/ NIN,KOR,NC,NPD,NCPTS,NEQN,IQUAD COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10, * IW11,IW12,IW13,IW14,IW15,IW16,IW17,IW18 COMMON /IOUNIT/ LOUT DIMENSION WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))* * (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))), * IWORK((NPDE+1)*(KORD+(NINT-1)*(KORD-NCC))),XBKPT(NINT+1) IF (INDEX .EQ. 0) GO TO 60 IF (INDEX .EQ. 2) GO TO 70 IF (INDEX .EQ. 4) GO TO 80 IF (INDEX .EQ. 3) GO TO 90 C----------------------------------------------------------------------- C SEVERAL CHECKS ARE MADE HERE TO DETERMINE IF THE INPUT PARAMETERS C HAVE LEGAL VALUES. ERROR CHECKS ARE MADE ON INDEX, EPS, (T0-TOUT)*DT, C DT, NINT, KORD, NCC, NPDE, MF, WHETHER THE BREAKPOINT SEQUENCE IS C STRICTLY INCREASING, AND WHETHER THERE IS SUFFICIENT STORAGE C PROVIDED FOR WORK AND IWORK. PROBLEM DEPENDENT PARAMETERS ARE C CALCULATED AND PLACED IN COMMON. C----------------------------------------------------------------------- IERID = -6 IF (INDEX .NE. 1) GO TO 320 IERID = IERID - 1 IF (EPS .LE. 0.) GO TO 320 IERID = IERID - 1 IF ((T0-TOUT)*DT .GT. 0.) GO TO 320 IERID = IERID - 1 IF (DT .EQ. 0.0) GO TO 320 IERID = IERID - 1 NIN = NINT IF (NIN .LT. 1) GO TO 320 IERID = IERID - 1 KOR = KORD IF (KOR .LT. 3 .OR. KOR .GT. 20) GO TO 320 IERID = IERID - 1 NC = NCC IF (NCC .LT. 2 .OR. NCC .GE. KOR) GO TO 320 IERID = IERID - 1 NPD = NPDE NPDE2 = NPD*NPD IF (NPDE .LT. 1) GO TO 320 IERID = IERID - 1 IF (MF.NE.22.AND.MF.NE.21.AND.MF.NE.12.AND.MF.NE.11) GO TO 320 IERID = IERID - 1 DO 10 K=1,NIN IF(XBKPT(K) .GE. XBKPT(K+1)) GO TO 320 10 CONTINUE NCPTS = KOR + (NIN - 1) * (KOR - NCC) NEQN = NPDE * NCPTS ML = (KOR-1)*NPDE - 1 MU = ML MW = ML + ML + 1 N0W = NEQN*MW IWSAVE = IWORK(1) IISAVE = IWORK(2) IW1 = 4*NPDE2 + 1 IW2 = IW1 + 3*KORD*NCPTS IW3 = IW2 + NCPTS + KORD IW4 = IW3 + NCPTS IW5 = IW4 + NEQN IW6 = IW5 + NEQN IW7 = IW6 + NEQN IW8 = IW7 + NEQN IW9 = IW8 + NEQN IW10 = IW9 + 3*NPDE IW11 = IW10 + NEQN*(MAXDER+1) IW12 = IW11 + NPDE2 IW13 = IW12 + NPDE2 IW14 = IW13 + NPDE2 IW15 = IW14 + NPDE2 IW16 = IW15 + NPDE2 IW17 = IW16 + NPDE IW18 = NCPTS + 1 IERID = IERID - 1 IWSTOR = IW17 + NEQN*(3*ML+1) - 1 IISTOR = IW18 + NEQN - 1 IF ( IWSAVE .LT. IWSTOR .OR. IISAVE .LT. IISTOR ) GO TO 335 C----------------------------------------------------------------------- C PERFORM INITIALIZATION TASKS. IF KORD .EQ. 3 THEN CALCULATE THE BAND- C WIDTH OF THE ASSOCIATED MATRIX PROBLEM BY DETERMINING THE TYPE OF C BOUNDARY CONDITIONS, THEN CHECK FOR SUFFICIENT STORAGE AGAIN. C----------------------------------------------------------------------- CALL INITAL(KOR,WORK(IW1),WORK(IW6),XBKPT,WORK(IW2),WORK(IW3), * WORK(IW17),IWORK(IW18),IWORK) IF(IQUAD .NE. 0) GO TO 280 IF( KOR .NE. 3 ) GO TO 40 CALL EVAL(1,NPDE,WORK(IW6),WORK(IW9),WORK(IW1),IWORK) CALL BNDRY(T0,WORK(IW3),WORK(IW9),WORK(IW9+NPDE),WORK(IW14), * WORK(IW15),WORK(IW16),NPDE) DO 20 K=1,NPDE I = K + NPDE*(K-1) - 1 IF(WORK(IW14+I) .EQ. 0.0 .AND. WORK(IW15+I) .EQ. 0.0) * IQUAD = 1 20 CONTINUE CALL EVAL(NCPTS,NPDE,WORK(IW6),WORK(IW9),WORK(IW1),IWORK) CALL BNDRY(T0,WORK(IW3+NCPTS-1),WORK(IW9),WORK(IW9+NPDE), * WORK(IW14),WORK(IW15),WORK(IW16),NPDE) DO 30 K=1,NPDE I = K + NPDE*(K-1) - 1 IF(WORK(IW14+I) .EQ. 0.0 .AND. WORK(IW15+I) .EQ. 0.0) * IQUAD = 1 30 CONTINUE ML = ML + IQUAD*NPDE MU = ML MW = ML + ML + 1 N0W = NEQN*MW 40 CONTINUE IWSTOR = IW17 + NEQN*(3*ML+1) - 1 IF ( IWSAVE .LT. IWSTOR ) GO TO 335 C----------------------------------------------------------------------- C IF INITIAL VALUES OF CMAX OTHER THAN THOSE SET BELOW ARE DESIRED, C THEY SHOULD BE SET HERE. ALL CMAX(I) MUST BE POSITIVE. C HAVING PROPER VALUES OF CMAX FOR THE PROBLEM BEING SOLVED IS AS C IMPORTANT AS CHOOSING EPS (SEE ABOVE), SINCE ERRORS ARE C MEASURED RELATIVE TO CMAX. IF VALUES FOR DTMN OR DTMX, THE C BOUNDS ON ABS(DT), OTHER THAN THOSE BELOW ARE DESIRED, THEY C SHOULD BE SET BELOW. C----------------------------------------------------------------------- DO 50 I = 1,NEQN I1 = I - 1 WORK(IW4+I1) = ABS(WORK(IW6+I1)) IF (WORK(IW4+I1) .LT. 1.) WORK(IW4+I1) = 1. 50 WORK(IW10+I1) = WORK(IW6+I1) N = NEQN T = T0 DTC = DT DTMN = ABS(DT) DTUSED = 0. EPSC = EPS MFC = MF JSTART = 0 EPSJ = SQRT(UROUND) NM1 = NEQN - 1 N0ML = NEQN*ML NHCUT = 0 KFLAG = 0 TOUTP = T0 IF ( T0 .EQ. TOUT ) GO TO 360 60 DTMX = ABS(TOUT-TOUTP)*10. GO TO 140 C 70 DTMX = ABS(TOUT-TOUTP)*10. IF ((T-TOUT)*DTC .GE. 0.) GO TO 340 GO TO 150 C 80 IF ((T-TOUT)*DTC .GE. 0.) GO TO 300 JSTART = -1 EPSC = EPS MFC = MF GO TO 100 C 90 DTMX = DT 100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110) 110 FORMAT(36H WARNING.. T + DT = T ON NEXT STEP.) C----------------------------------------------------------------------- C TAKE A TIME STEP BY CALLING THE INTEGRATOR. C----------------------------------------------------------------------- CALL STIFIB (NEQN,WORK(IW10),WORK(IW4),WORK(IW5),WORK(IW6), * WORK(IW7),WORK(IW8),WORK(IW17),IWORK(IW18),WORK,IWORK) C KGO = 1 - KFLAG GO TO (120, 160, 220, 260, 280), KGO C KFLAG = 0, -1, -2, -3 -4 C 120 CONTINUE C----------------------------------------------------------------------- C NORMAL RETURN FROM INTEGRATOR. C C THE WEIGHTS CMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED, C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL C FOR THE MACHINE PRECISION. C C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY C STEP SHOULD BE INSERTED HERE. C C IF INDEX = 3, SAVE1 IS SET TO THE CURRENT C VALUES ON RETURN. C IF INDEX = 2, DT IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF C ERROR), AND THEN THE CURRENT C VALUES ARE PUT IN SAVE1 ON RETURN. C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF C ARE C COMPUTED AND STORED IN SAVE1 ON RETURN. C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 340 INSTEAD OF 360. C----------------------------------------------------------------------- D = 0. DO 130 I = 1,NEQN I1 = I - 1 AYI = ABS(WORK(IW10+I1)) WORK(IW4+I1) = AMAX1(WORK(IW4+I1), AYI) 130 D = D + (AYI/WORK(IW4+I1))**2 D = D*(UROUND/EPS)**2 IF (D .GT. FLOAT(NEQN)) GO TO 240 IF (INDEX .EQ. 3) GO TO 340 IF (INDEX .EQ. 2) GO TO 150 140 IF ((T-TOUT)*DTC .LT. 0.) GO TO 100 CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6)) GO TO 360 C 150 IF (((T+DTC)-TOUT)*DTC .LE. 0.) GO TO 100 IF (ABS(T-TOUT) .LE. 100.*UROUND*DTMX) GO TO 340 IF ((T-TOUT)*DTC .GE. 0.) GO TO 340 DTC = (TOUT - T)*(1. - 4.*UROUND) JSTART = -1 GO TO 100 C----------------------------------------------------------------------- C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF C KFLAG = -2 OR -4, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE. C TO RECOVER, DT AND DTMN ARE REDUCED BY A FACTOR OF .1 UP TO 10 C TIMES BEFORE GIVING UP. C----------------------------------------------------------------------- 160 WRITE (LOUT,170) T 170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/ * 40H ERROR TEST FAILED WITH ABS(DT) = DTMIN/) 180 IF (NHCUT .EQ. 10) GO TO 200 NHCUT = NHCUT + 1 DTMN = .1*DTMN DTC = .1*DTC WRITE (LOUT,190) DTC 190 FORMAT(25H DT HAS BEEN REDUCED TO ,E16.8, * 26H AND STEP WILL BE RETRIED//) JSTART = -1 GO TO 100 C 200 WRITE (LOUT,210) 210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) GO TO 340 C 220 WRITE (LOUT,230) T,DTC 230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H DT =, * E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) GO TO 340 C 240 WRITE (LOUT,250) T 250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/ * 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) KFLAG = -2 GO TO 340 C 260 WRITE (LOUT,270) T 270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/ * 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/) GO TO 180 C 280 WRITE (LOUT,290) 290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED, * 35H PROBABLY DUE TO STORAGE OVERWRITES//) KFLAG = -4 GO TO 340 C 300 WRITE(LOUT,310) T,TOUT,DTC 310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./ * 4H T =,E16.8,9H TOUT =,E16.8,8H DTC =,E16.8/ * 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./ * 41H DESIRED PARAMETER CHANGES WERE NOT MADE.) CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6)) INDEX = -5 RETURN C 320 WRITE(LOUT,330) IERID 330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//) INDEX = IERID RETURN C 335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE 336 FORMAT(//21H INSUFFICIENT STORAGE/24H WORK MUST BE OF LENGTH, * I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X, * 12HYOU PROVIDED,I10//) INDEX = IERID RETURN C 340 TOUT = T DO 350 I = 1,NEQN I1 = I - 1 350 WORK(IW6+I1) = WORK(IW10+I1) 360 INDEX = KFLAG TOUTP = TOUT DT = DTUSED IF (KFLAG .NE. 0) DT = DTC RETURN END SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK) C----------------------------------------------------------------------- C SUBROUTINE VALUES COMPUTES THE SOLUTION U AND THE FIRST NDERV C DERIVATIVES OF U AT THE NPTS POINTS X AND AT TIME TOUT AND RETURNS C THEM IN THE ARRAY USOL. THIS ROUTINE MUST BE USED TO OBTAIN C SOLUTION VALUES SINCE PDECOL DOES NOT RETURN ANY SOLUTION VALUES C TO THE USER. SEE PDECOL. C C THE CALLING PARAMETERS ARE... C X = AN ARBITRARY VECTOR OF SPATIAL POINTS OF LENGTH NPTS AT C WHICH THE SOLUTION AND THE FIRST NDERV DERIVATIVE VALUES C ARE TO BE CALCULATED. IF X .LT. XLEFT ( X .GT. XRIGHT ) C THEN THE PIECEWISE POLYNOMIAL OVER THE LEFTMOST ( RIGHT- C MOST ) INTERVAL IS EVALUATED TO CALCULATE THE SOLUTION C VALUES AT THIS UNUSUAL VALUE OF X. SEE PDECOL. C C USOL = AN ARRAY WHICH CONTAINS THE SOLUTION AND THE FIRST C NDERV DERIVATIVES OF THE SOLUTION AT ALL THE POINTS IN C THE INPUT VECTOR X. IN PARTICULAR, USOL(J,I,K) CONTAINS C THE VALUE OF THE (K-1)-ST DERIVATIVE OF THE J-TH PDE C COMPONENT AT THE I-TH POINT OF THE X VECTOR FOR C J = 1 TO NPDE, I = 1 TO NPTS, AND K = 1 TO NDERV+1. C C SCTCH = A USER SUPPLIED WORKING STORAGE ARRAY OF LENGTH AT LEAST C KORD*(NDERV+1). SEE BELOW AND PDECOL FOR DEFINITIONS OF C THESE PARAMETERS. C C NDIM1 = THE FIRST DIMENSION OF THE OUTPUT ARRAY USOL IN THE CALLING C PROGRAM. NDIM1 MUST BE .GE. NPDE. C C NDIM2 = THE SECOND DIMENSION OF THE OUTPUT ARRAY USOL IN THE C CALLING PROGRAM. NDIM2 MUST BE .GE. NPTS. C C NPTS = THE NUMBER OF POINTS IN THE X VECTOR. C C NDERV = THE NUMBER OF DERIVATIVE VALUES OF THE SOLUTION THAT ARE C TO BE CALCULATED. NDERV SHOULD BE LESS THAN KORD SINCE C THE KORD-TH DERIVATIVE OF A POLYNOMIAL OF DEGREE KORD-1 C IS EQUAL TO ZERO. SEE PDECOL. C C WORK = THE USERS WORKING STORAGE ARRAY WHICH IS USED IN THIS CASE C TO PROVIDE THE CURRENT BASIS FUNCTION COEFFICIENTS AND THE C PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE. C C PACKAGE ROUTINES CALLED.. BSPLVD,INTERV C USER ROUTINES CALLED.. NONE C CALLED BY.. USERS MAIN PROGRAM C FORTRAN FUNCTIONS USED.. NONE C C----------------------------------------------------------------------- SAVE ILEFT, MFLAG COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IDUM(12) COMMON /OPTION/ NOGAUS,MAXDER DIMENSION USOL(NDIM1,NDIM2,NDERV+1),X(NPTS),SCTCH(KORD*(NDERV+1)), * WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))* * (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))) DATA ILEFT/0/, MFLAG/0/ NDERV1 = NDERV + 1 DO 20 IPTS=1,NPTS CALL INTERV(WORK(IW2),NCPTS,X(IPTS),ILEFT,MFLAG) CALL BSPLVD(WORK(IW2),KORD,X(IPTS),ILEFT,SCTCH,NDERV1) IK = ILEFT - KORD DO 10 M=1,NDERV1 I1 = (M-1)*KORD DO 10 K=1,NPDE USOL(K,IPTS,M) = 0. DO 10 I=1,KORD I2 = (I+IK-1)*NPDE + IW6 - 1 USOL(K,IPTS,M) = USOL(K,IPTS,M) + WORK(I2+K) * SCTCH(I+I1) 10 CONTINUE 20 CONTINUE RETURN END BLOCK DATA C----------------------------------------------------------------------- C IN THE FOLLOWING DATA STATEMENT, SET.. C LOUT = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES DURING C THE INTEGRATION. C NOGAUS = SET .EQ. 1 IF THE GAUSS-LEGENDRE COLLOCATION POINTS ARE C NOT DESIRED WHEN NCC = 2 (SEE PDECOL AND COLPNT). C MAXDER = SET .EQ. 5. ITS VALUE REPRESENTS THE MAXIMUM ORDER OF C THE TIME INTEGRATION ALLOWED. ITS VALUE AFFECTS THE STOR- C AGE REQUIRED IN WORK AND MAY BE CHANGED IF DESIRED C (SEE COSET FOR RESTRICTIONS). C UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST C POSITIVE U SUCH THAT 1. + U .NE. 1. ON THE MACHINE. C----------------------------------------------------------------------- COMMON /GEAR1/ DUM(5),UROUND,IDUM(4) COMMON /OPTION/ NOGAUS,MAXDER COMMON /IOUNIT/ LOUT C*** C*** UROUND SET TO SINGLE PRECISION FOR A SUN SPARC2 C*** DATA LOUT,NOGAUS,MAXDER,UROUND/6,0,5,5.960464E-08/ END SUBROUTINE INITAL(K,A,RHS,X,XT,XC,PW,IPIV,ILEFT) C----------------------------------------------------------------------- C INITAL IS CALLED ONLY ONCE BY PDECOL TO PERFORM INITIALIZATION TASKS. C THESE TASKS INCLUDE - 1) DEFINING THE PIECEWISE POLYNOMIAL SPACE C BREAKPOINT SEQUENCE, 2) CALLING THE SUBROUTINE COLPNT TO DEFINE THE C REQUIRED COLLOCATION POINTS, 3) DEFING THE PIECEWISE POLYNOMIAL SPACE C BASIS FUNCTION VALUES (PLUS FIRST AND SECOND DERIVATIVE VALUES) AT C THE COLLOCATION POINTS, AND 4) DEFINING THE INITIAL BASIS FUNCTION C COEFFICIENTS WHICH DETERMINE THE PIECEWISE POLYNOMIAL WHICH C INTERPOLATES THE USER SUPPLIED (UINIT) INITIAL CONDITION FUNCTION(S) C AT THE COLLOCATION POINTS. C C K = ORDER OF PIECEWISE POLYNOMIAL SPACE. C A = BASIS FUNCTION VALUES GENERATED BY INITAL. C RHS = TEMPORARY STORAGE USED TO RETURN INITIAL CONDITION COEFFICIENT C VALUES. C X = USER DEFINED PIECEWISE POLYNOMIAL BREAKPOINTS. C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE GENERATED BY INITAL. C XC = COLLOCATION POINTS GENERATED BY INITAL. C PW = STORAGE FOR BAND MATRIX USED TO GENERATE INITIAL C COEFFICIENT VALUES. C IPIV = PIVOT INFORMATION FOR LINEAR EQUATION SOLVER DECB-SOLB. C ILEFT = POINTERS TO BREAKPOINT SEQUENCE GENERATED BY INITAL. C C PACKAGE ROUTINES CALLED.. BSPLVD,COLPNT,DECB,INTERV,SOLB C USER ROUTINES CALLED.. UINIT C CALLED BY.. PDECOL C FORTRAN FUNCTIONS USED.. MAX0,MIN0 C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IER COMMON /GEAR9/ EPSJ,R0,ML,MU,IDUM(3),N0W DIMENSION A(K,3,NCPTS),RHS(NEQN),X(NINT+1),XT(NCPTS+KORD), * XC(NCPTS),PW(NEQN*(3*ML+1)), * IPIV(NEQN),ILEFT(NCPTS) MFLAG = -2 IER = 0 C----------------------------------------------------------------------- C SET UP THE PIECEWISE POLYNOMIAL SPACE BREAKPOINT SEQUENCE. C----------------------------------------------------------------------- KRPT = KORD - NCC DO 10 I=1,KORD XT(NCPTS+I) = X(NINT+1) 10 XT(I) = X(1) DO 20 I=2,NINT I1 = (I-2)*KRPT + KORD DO 20 J=1,KRPT 20 XT(I1+J) = X(I) C----------------------------------------------------------------------- C SET UP COLLOCATION POINTS ARRAY XC. C----------------------------------------------------------------------- CALL COLPNT(X, XC, XT) C----------------------------------------------------------------------- C GENERATE THE ILEFT ARRAY. STORE THE BASIS FUNCTION VALUES IN THE C ARRAY A. THE ARRAY A IS DIMENSIONED A(KORD,3,NCPTS) AND A(K,J,I) C CONTAINS THE VALUE OF THE (J-1)-ST DERIVATIVE (J = 1,2,3) OF THE K-TH C NONZERO BASIS FUNCTION (K = 1, ... ,KORD) AT THE I-TH COLLOCATION C POINT (I = 1, ... ,NCPTS). SET UP RHS FOR INTERPOLATING THE INITIAL C CONDITIONS AT THE COLLOCATION POINTS. SET THE INTERPOLATION MATRIX C INTO THE BANDED MATRIX PW. C----------------------------------------------------------------------- DO 30 I=1,N0W 30 PW(I) = 0. DO 40 I=1,NCPTS CALL INTERV(XT,NCPTS,XC(I),ILEFT(I),MFLAG) CALL BSPLVD(XT,KORD,XC(I),ILEFT(I),A(1,1,I),3) I1 = NPDE * (I-1) CALL UINIT(XC(I),RHS(I1+1),NPDE) ICOL = ILEFT(I) - I - 1 JL = MAX0(1,I+2-NCPTS) JU = MIN0(KORD,KORD+I-2) DO 40 J=JL,JU J1 = I1 + NEQN * (NPDE * (ICOL + J) - 1) DO 40 JJ=1,NPDE 40 PW(JJ+J1) = A(J,1,I) C----------------------------------------------------------------------- C LU DECOMPOSE THE MATRIX PW. C----------------------------------------------------------------------- CALL DECB (NEQN,NEQN,ML,MU,PW,IPIV,IER) IF ( IER .NE. 0 ) RETURN C----------------------------------------------------------------------- C SOLVE THE LINEAR SYSTEM PW*Z = RHS. THIS GIVES THE BASIS FUNCTION C COEFFICIENTS FOR THE INITIAL CONDITIONS. C----------------------------------------------------------------------- CALL SOLB (NEQN,NEQN,ML,MU,PW,RHS,IPIV) RETURN END SUBROUTINE COLPNT(X, XC, XT) C----------------------------------------------------------------------- C COLPNT IS CALLED ONLY ONCE BY INITAL TO DEFINE THE REQUIRED COLLOCA- C TION POINTS WHICH ARE TO BE USED WITH THE USER SELECTED PIECEWISE C POLYNOMIAL SPACE. THE COLLOCATION POINTS ARE CHOSEN SUCH THAT THEY C ARE EITHER THE POINTS AT WHICH THE PIECEWISE POLYNOMIAL SPACE BASIS C FUNCTIONS ATTAIN THEIR UNIQUE MAXIMUM VALUES, OR, THE GAUSS-LEGENDRE C QUADRATURE POINTS WITHIN EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL, C DEPENDING UPON THE SPACE BEING USED AND THE DESIRE OF THE USER. C C X = USER DEFINED PIECEWISE POLYNOMIAL BREAKPOINTS. C XC = COLLOCATION POINTS DEFINED BY COLPNT. C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE. C C PACKAGE ROUTINES CALLED.. BSPLVD,INTERV C USER ROUTINES CALLED.. NONE C CALLED BY.. INITAL C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- SAVE ILEFT COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /OPTION/ NOGAUS,MAXDER DIMENSION RHO(40),X(NINT+1),XC(NCPTS),XT(NCPTS+KORD) DATA ILEFT/0/ C----------------------------------------------------------------------- C IF THE VARIABLE NOGAUS IN THE COMMON BLOCK /OPTION/ IS SET .EQ. 1, C THE USE OF THE GAUSS-LEGENDRE POINTS IS PROHIBITED FOR ALL CASES. C NOGAUS IS CURRENTLY SET .EQ. 0 BY A DATA STATEMENT IN THE BLOCK DATA. C THE USER MAY CHANGE THIS AS DESIRED. C----------------------------------------------------------------------- IF ( NCC .NE. 2 .OR. NOGAUS .EQ. 1 ) GO TO 200 C----------------------------------------------------------------------- C COMPUTE THE COLLOCATION POINTS TO BE AT THE GAUSS-LEGENDRE POINTS IN C EACH PIECEWISE POLYNOMIAL SPACE SUBINTERVAL. THE ARRAY RHO IS SET TO C CONTAIN THE GAUSS-LEGENDRE POINTS FOR THE STANDARD INTERVAL (-1,1). C----------------------------------------------------------------------- IPTS = KORD - 2 GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170, * 180),IPTS 10 RHO(1) = 0. GO TO 190 20 RHO(2) = .577350269189626E-00 RHO(1) = - RHO(2) GO TO 190 30 RHO(3) = .774596669241483E-00 RHO(1) = - RHO(3) RHO(2) = 0. GO TO 190 40 RHO(3) = .339981043584856E-00 RHO(2) = - RHO(3) RHO(4) = .861136311594053E-00 RHO(1) = - RHO(4) GO TO 190 50 RHO(4) = .538469310105683E-00 RHO(2) = - RHO(4) RHO(5) = .906179845938664E-00 RHO(1) = - RHO(5) RHO(3) = 0. GO TO 190 60 RHO(4) = .238619186083197E-00 RHO(3) = - RHO(4) RHO(5) = .661209386466265E-00 RHO(2) = - RHO(5) RHO(6) = .932469514203152E-00 RHO(1) = - RHO(6) GO TO 190 70 RHO(5) = .405845151377397E-00 RHO(3) = - RHO(5) RHO(6) = .741531185599394E-00 RHO(2) = - RHO(6) RHO(7) = .949107912342759E-00 RHO(1) = - RHO(7) RHO(4) = 0. GO TO 190 80 RHO(5) = .183434642495650E-00 RHO(4) = - RHO(5) RHO(6) = .525532409916329E-00 RHO(3) = - RHO(6) RHO(7) = .796666477413627E-00 RHO(2) = - RHO(7) RHO(8) = .960289856497536E-00 RHO(1) = - RHO(8) GO TO 190 90 RHO( 5) = .0 RHO( 6) = .324253423403809E-00 RHO( 7) = .613371432700590E-00 RHO( 8) = .836031107326636E-00 RHO( 9) = .968160239507626E-00 DO 95 I=1,4 95 RHO(I) = -RHO(10-I) GO TO 190 100 RHO( 6) = .148874338981631E-00 RHO( 7) = .433395394129247E-00 RHO( 8) = .679409568299024E-00 RHO( 9) = .865063366688984E-00 RHO(10) = .973906528517172E-00 DO 105 I=1,5 105 RHO(I) = -RHO(11-I) GO TO 190 110 RHO( 6) = .0 RHO( 7) = .269543155952345E-00 RHO( 8) = .519096129206812E-00 RHO( 9) = .730152005574049E-00 RHO(10) = .887062599768095E-00 RHO(11) = .978228658146057E-00 DO 115 I=1,5 115 RHO(I) = -RHO(12-I) GO TO 190 120 RHO( 7) = .125233408511469E-00 RHO( 8) = .367831498998180E-00 RHO( 9) = .587317954286617E-00 RHO(10) = .769902674194305E-00 RHO(11) = .904117256370475E-00 RHO(12) = .981560634246719E-00 DO 125 I=1,6 125 RHO(I) = -RHO(13-I) GO TO 190 130 RHO( 7) = .0 RHO( 8) = .230458315955135E-00 RHO( 9) = .448492751036447E-00 RHO(10) = .642349339440340E-00 RHO(11) = .801578090733310E-00 RHO(12) = .917598399222978E-00 RHO(13) = .984183054718588E-00 DO 135 I=1,6 135 RHO(I) = -RHO(14-I) GO TO 190 140 RHO( 8) = .108054948707344E-00 RHO( 9) = .319112368927890E-00 RHO(10) = .515248636358154E-00 RHO(11) = .687292904811685E-00 RHO(12) = .827201315069765E-00 RHO(13) = .928434883663574E-00 RHO(14) = .986283808696812E-00 DO 145 I=1,7 145 RHO(I) = -RHO(15-I) GO TO 190 150 RHO( 8) = .0 RHO( 9) = .201194093997435E-00 RHO(10) = .394151347077563E-00 RHO(11) = .570972172608539E-00 RHO(12) = .724417731360170E-00 RHO(13) = .848206583410427E-00 RHO(14) = .937273392400706E-00 RHO(15) = .987992518020485E-00 DO 155 I = 1,7 155 RHO(I) = -RHO(16-I) GO TO 190 160 RHO( 9) = .950125098376374E-01 RHO(10) = .281603550779259E-00 RHO(11) = .458016777657227E-00 RHO(12) = .617876244402644E-00 RHO(13) = .755404408355003E-00 RHO(14) = .865631202387832E-00 RHO(15) = .944575023073233E-00 RHO(16) = .989400934991650E-00 DO 165 I=1,8 165 RHO(I) = -RHO(17-I) GO TO 190 170 RHO( 9) = .0 RHO(10) = .178484181495848E-00 RHO(11) = .351231763453876E-00 RHO(12) = .512690537086477E-00 RHO(13) = .657671159216691E-00 RHO(14) = .781514003896801E-00 RHO(15) = .880239153726986E-00 RHO(16) = .950675521768768E-00 RHO(17) = .990575475314417E-00 DO 175 I=1,8 175 RHO(I) = -RHO(18-I) GO TO 190 180 RHO(10) = .847750130417353E-01 RHO(11) = .251886225691506E-00 RHO(12) = .411751161462843E-00 RHO(13) = .559770831073948E-00 RHO(14) = .691687043060353E-00 RHO(15) = .803704958972523E-00 RHO(16) = .892602466497556E-00 RHO(17) = .955823949571398E-00 RHO(18) = .991565168420931E-00 DO 185 I=1,9 185 RHO(I) = -RHO(19-I) C----------------------------------------------------------------------- C COMPUTE THE GAUSS-LEGENDRE COLLOCATION POINTS IN EACH SUBINTERVAL. C----------------------------------------------------------------------- 190 DO 195 I=1,NINT FAC = ( X(I+1) - X(I) ) * .5 DO 195 J = 1,IPTS KNOT = IPTS * (I-1) + J + 1 195 XC(KNOT) = X(I) + FAC * ( RHO(J) + 1. ) XC(1) = X(1) XC(NCPTS) = X(NINT+1) RETURN C----------------------------------------------------------------------- C COMPUTE THE COLLOCATION POINTS TO BE AT THE POINTS WHERE THE BASIS C FUNCTIONS ATTAIN THEIR MAXIMA. A BISECTION METHOD IS USED TO FIND C THE POINTS TO MACHINE PRECISION. THIS PROCESS COULD BE SPEEDED UP C BY USING A SECANT METHOD IF DESIRED. C----------------------------------------------------------------------- 200 ITOP = NCPTS - 1 MFLAG = -2 XC(1) = X(1) XC(NCPTS) = X(NINT+1) DO 240 I=2,ITOP XOLD = 1.E+20 XL = XT(I) XR = XT(I+KORD) 210 XNEW = .5 * (XL + XR) IF( XOLD .EQ. XNEW ) GO TO 240 CALL INTERV(XT,NCPTS,XNEW,ILEFT,MFLAG) CALL BSPLVD(XT,KORD,XNEW,ILEFT,RHO,2) DO 220 J=1,KORD IF( I .EQ. J + ILEFT - KORD ) GO TO 230 220 CONTINUE 230 XVAL = RHO(KORD+J) IF( XVAL .EQ. 0.0 ) XR = XNEW IF( XVAL .GT. 0.0 ) XL = XNEW IF( XVAL .LT. 0.0 ) XR = XNEW XOLD = XNEW GO TO 210 240 XC(I) = XR RETURN END SUBROUTINE BSPLVD ( XT, K, X, ILEFT, VNIKX, NDERIV ) C----------------------------------------------------------------------- C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE. C SEE REFERENCE BELOW. C C CALCULATES THE VALUE AND THE FIRST NDERIV-1 DERIVATIVES OF ALL C B-SPLINES WHICH DO NOT VANISH AT X. THE ROUTINE FILLS THE TWO- C DIMENSIONAL ARRAY VNIKX(J,IDERIV), J=IDERIV, ... ,K WITH NONZERO C VALUES OF B-SPLINES OF ORDER K+1-IDERIV, IDERIV=NDERIV, ... ,1, BY C REPEATED CALLS TO BSPLVN. C C XT = PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE. C K = ORDER OF THE PIECEWISE POLYNOMIAL SPACE. C X = POINT AT WHICH THE B-SPLINE IS TO BE EVALUATED. C ILEFT = POINTER TO THE BREAKPOINT SEQUENCE. C VNIKX = TABLE OF B-SPLINE VALUES AND DERIVATIVES. C NDERIV = DETERMINES NUMBER OF DERIVATIVES TO BE GENERATED. C C REFERENCE C C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J. C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472. C C PACKAGE ROUTINES CALLED.. BSPLVN C USER ROUTINES CALLED.. NONE C CALLED BY.. COLPNT,INITAL,VALUES C FORTRAN FUNCTIONS USED.. FLOAT,MAX0 C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD DIMENSION XT(NCPTS+KORD),VNIKX(K,NDERIV) DIMENSION A(20,20) KO = K + 1 - NDERIV CALL BSPLVN(XT,KO,1,X,ILEFT,VNIKX(NDERIV,NDERIV)) IF (NDERIV .LE. 1) GO TO 120 IDERIV = NDERIV DO 20 I=2,NDERIV IDERVM = IDERIV-1 DO 10 J=IDERIV,K 10 VNIKX(J-1,IDERVM) = VNIKX(J,IDERIV) IDERIV = IDERVM CALL BSPLVN(XT,0,2,X,ILEFT,VNIKX(IDERIV,IDERIV)) 20 CONTINUE DO 40 I=1,K DO 30 J=1,K 30 A(I,J) = 0. 40 A(I,I) = 1. KMD = K DO 110 M=2,NDERIV KMD = KMD - 1 FKMD = FLOAT(KMD) I = ILEFT J = K 50 JM1 = J-1 IPKMD = I + KMD DIFF = XT(IPKMD) -XT(I) IF (JM1 .EQ. 0) GO TO 80 IF (DIFF .EQ. 0.) GO TO 70 DO 60 L=1,J 60 A(L,J) = (A(L,J) - A(L,J-1))/DIFF*FKMD 70 J = JM1 I = I - 1 GO TO 50 80 IF (DIFF .EQ. 0.) GO TO 90 A(1,1) = A(1,1)/DIFF*FKMD 90 DO 110 I=1,K V = 0. JLOW = MAX0(I,M) DO 100 J=JLOW,K 100 V = A(I,J)*VNIKX(J,M) + V 110 VNIKX(I,M) = V 120 RETURN END SUBROUTINE BSPLVN ( XT, JHIGH, INDEX, X, ILEFT, VNIKX ) C----------------------------------------------------------------------- C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE. C SEE REFERENCE BELOW. C C CALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT THE C POINT X OF ORDER MAX(JHIGH,(J+1)(INDEX-1)) FOR THE BREAKPOINT SEQ- C UENCE XT. ASSUMING THAT XT(ILEFT) .LE. X .LE. XT(ILEFT+1), THE ROUT- C INE RETURNS THE B-SPLINE VALUES IN THE ONE DIMENSIONAL ARRAY VNIKX. C C FOR DEFINITIONS OF CALLING ARGUMENTS SEE ABOVE AND BSPLVD. C C REFERENCE C C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J. C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. BSPLVD C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- SAVE J,DELTAM,DELTAP COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD DIMENSION DELTAM(20),DELTAP(20) DIMENSION XT(NCPTS+KORD),VNIKX(*) DATA J/1/,DELTAM/20*0.E-00/,DELTAP/20*0.E-00/ GO TO (10,20),INDEX 10 J = 1 VNIKX(1) = 1. IF (J .GE. JHIGH) GO TO 40 20 IPJ = ILEFT+J DELTAP(J) = XT(IPJ) - X IMJP1 = ILEFT-J+1 DELTAM(J) = X - XT(IMJP1) VMPREV = 0. JP1 = J+1 DO 30 L=1,J JP1ML = JP1-L VM = VNIKX(L)/(DELTAP(L) + DELTAM(JP1ML)) VNIKX(L) = VM*DELTAP(L) + VMPREV 30 VMPREV = VM*DELTAM(JP1ML) VNIKX(JP1) = VMPREV J = JP1 IF (J .LT. JHIGH) GO TO 20 40 RETURN END SUBROUTINE INTERV ( XT, LXT, X, ILEFT, MFLAG ) C----------------------------------------------------------------------- C THIS SUBROUTINE IS PART OF THE B-SPLINE PACKAGE FOR THE STABLE C EVALUATION OF ANY B-SPLINE BASIS FUNCTION OR DERIVATIVE VALUE. C SEE REFERENCE BELOW. C C COMPUTES LARGEST ILEFT IN (1,LXT) SUCH THAT XT(ILEFT) .LE. X. THE C PROGRAM STARTS THE SEARCH FOR ILEFT WITH THE VALUE OF ILEFT THAT WAS C RETURNED AT THE PREVIOUS CALL (AND WAS SAVED IN THE LOCAL VARIABLE C ILO) TO MINIMIZE THE WORK IN THE COMMON CASE THAT THE VALUE OF X ON C THIS CALL IS CLOSE TO THE VALUE OF X ON THE PREVIOUS CALL. SHOULD C THIS ASSUMPTION NOT BE VALID, THEN THE PROGRAM LOCATES ILO AND IHI C SUCH THAT XT(ILO) .LE. X .LT. XT(IHI) AND, ONCE THEY ARE FOUND USES C BISECTION TO FIND THE CORRECT VALUE FOR ILEFT. MFLAG IS AN ERROR FLAG. C C FOR DEFINITIONS OF CALLING ARGUMENTS SEE ABOVE AND BSPLVD. C C REFERENCE C C DEBOOR, C., PACKAGE FOR CALCULATING WITH B-SPLINES, SIAM J. C NUMER. ANAL., VOL. 14, NO. 3, JUNE 1977, PP. 441-472. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. COLPNT,INITAL,VALUES C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- SAVE ILO DIMENSION XT(LXT) IF(MFLAG.EQ.-2) ILO = 1 IHI = ILO + 1 IF (IHI .LT. LXT) GO TO 20 IF (X .GE. XT(LXT)) GO TO 110 IF (LXT .LE. 1) GO TO 90 ILO = LXT - 1 GO TO 21 20 IF (X .GE. XT(IHI)) GO TO 40 21 IF (X .GE. XT(ILO)) GO TO 100 C----------------------------------------------------------------------- C NOW X .LT. XT(IHI). FIND LOWER BOUND. C----------------------------------------------------------------------- 30 ISTEP = 1 31 IHI = ILO ILO = IHI - ISTEP IF (ILO .LE. 1) GO TO 35 IF (X .GE. XT(ILO)) GO TO 50 ISTEP = ISTEP*2 GO TO 31 35 ILO = 1 IF (X .LT. XT(1)) GO TO 90 GO TO 50 C----------------------------------------------------------------------- C NOW X .GE. XT(ILO). FIND UPPER BOUND. C----------------------------------------------------------------------- 40 ISTEP = 1 41 ILO = IHI IHI = ILO + ISTEP IF (IHI .GE. LXT) GO TO 45 IF (X .LT. XT(IHI)) GO TO 50 ISTEP = ISTEP*2 GO TO 41 45 IF (X .GE. XT(LXT)) GO TO 110 IHI = LXT C----------------------------------------------------------------------- C NOW XT(ILO) .LE. X .LT. XT(IHI). NARROW THE INTERVAL. C----------------------------------------------------------------------- 50 MIDDLE = (ILO + IHI)/2 IF (MIDDLE .EQ. ILO) GO TO 100 C----------------------------------------------------------------------- C NOTE.. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1. C----------------------------------------------------------------------- IF (X .LT. XT(MIDDLE)) GO TO 53 ILO = MIDDLE GO TO 50 53 IHI = MIDDLE GO TO 50 C----------------------------------------------------------------------- C SET OUTPUT AND RETURN. C----------------------------------------------------------------------- 90 MFLAG = -1 ILEFT = 1 RETURN 100 MFLAG = 0 ILEFT = ILO RETURN 110 MFLAG = 1 ILEFT = LXT RETURN END SUBROUTINE STIFIB (N0,Y,YMAX,ERROR,SAVE1,SAVE2,SAVE3, * PW,IPIV,WORK,IWORK) C----------------------------------------------------------------------- C STIFIB PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM, C A(Y,T)*(DY/DT) = G(Y,T), WHERE Y = (Y(1),Y(2), ... ,Y(N)). C STIFIB IS FOR USE WHEN THE MATRICES A AND DG/DY HAVE BANDED OR NEARLY C BANDED FORM. THE DEPENDENCE OF A(Y,T) ON Y IS ASSUMED TO BE WEAK. C C REFERENCE C C HINDMARSH, A.C., PRELIMINARY DOCUMENTATION OF GEARIB.. SOLUTION C OF IMPLICIT SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS WITH C BANDED JACOBIANS, LAWRENCE LIVERMORE LAB, UCID-30130, FEBRUARY C 1976. C C COMMUNICATION WITH STIFIB IS DONE WITH THE FOLLOWING VARIABLES.. C C Y AN N0 BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES C AND THEIR SCALED DERIVATIVES. LMAX IS 13 FOR THE ADAMS C METHODS AND 6 FOR THE GEAR METHODS. LMAX - 1 = MAXDER C IS THE MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET. C Y(I,J+1) CONTAINS THE J-TH DERIVATIVE OF Y(I), SCALED BY C H**J/FACTORIAL(J) (J = 0,1,...,NQ). C N0 A CONSTANT INTEGER .GE. N, USED FOR DIMENSIONING PURPOSES. C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN. C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. C HMIN, THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY C TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE. C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN PDECOL. C UROUND THE UNIT ROUNDOFF OF THE MACHINE. C N THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. C MF THE METHOD FLAG. SEE DESCRIPTION IN PDECOL. C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C 0 THE STEP WAS SUCCESFUL. C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED C WITH ABS(H) = HMIN. C -2 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR ABS(H) = HMIN. C -4 SINGULAR A-MATRIX ENCOUNTERED. C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND C THE Y ARRAY ARE AS OF THE BEGINNING OF THE LAST C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. C JSTART AN INTEGER USED ON INPUT AND OUTPUT. C ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS.. C 0 PERFORM THE FIRST STEP. C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF C H, EPS, N, AND/OR MF. C ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD. C YMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL C ERRORS IN Y ARE COMPARED. C ERROR AN ARRAY OF N ELEMENTS. ERROR(I)/TQ(2) IS THE ESTIMATED C ONE-STEP ERROR IN Y(I). C SAVE1,SAVE2,SAVE3 THREE WORKING STORAGE ARRAYS, EACH OF LENGTH N. C PW A BLOCK OF LOCATIONS USED FOR THE CHORD ITERATION C MATRIX. SEE DESCRIPTION IN PDECOL. C IPIV AN INTEGER ARRAY OF LENGTH N FOR PIVOT INFORMATION. C ML,MU THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF C THE CHORD ITERATION MATRIX. SEE DESCRIPTION IN PDECOL. C WORK,IWORK WORKING ARRAYS WHICH ARE USED TO PASS APPROPRIATE C ARRAYS TO OTHER SUBROUTINES. C C PACKAGE ROUTINES CALLED.. COSET,DIFFUN,PSETIB,RES,SOLB C USER ROUTINES CALLED.. NONE C CALLED BY.. PDECOL C FORTRAN FUNCTIONS USED.. ABS,AMAX1,AMIN1,FLOAT C----------------------------------------------------------------------- C SAVE EL,OLDL0,TQ,IER,NQ,L,METH,MITER SAVE BND,CON,CRATE,D,D1,E,EDN,ENQ1,ENQ2,ENQ3,EPSOLD, * EUP,FN,HOLD,OLDL0,PR1,PR2,PR3,R1,RC,RH,RMAX,TOLD, * I,IDOUB,IER,IREDO,IRET,IWEVAL,J,J1,J2,L,LMAX,M,MEO,METH, * MFOLD,MIO,MITER,NEWQ,NOLD,NQ,NSTEPJ COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /ISTART/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11, * IW12,IW13,IW14,IW15,IW16,IW17,IW18 COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,N,MF,KFLAG,JSTART COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W COMMON /GEAR0/ HUSED,NQUSED,NSTEP,NFE,NJE COMMON /OPTION/ NOGAUS,MAXDER DIMENSION Y(NEQN,MAXDER+1),YMAX(NEQN),ERROR(NEQN),SAVE1(NEQN), * SAVE2(NEQN), * SAVE3(NEQN),PW(NEQN*(3*ML+1)),IPIV(NEQN), * IWORK((NPDE+1)*NCPTS), * WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))* * (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))) DIMENSION EL(13),TQ(4) DATA EL(2)/1./, OLDL0/1./, TQ(1)/0./, IER/0/ KFLAG = 0 TOLD = T IF (JSTART .GT. 0) GO TO 200 IF (JSTART .NE. 0) GO TO 120 C----------------------------------------------------------------------- C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL YDOT IS C CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 C FOR THE NEXT INCREASE. C----------------------------------------------------------------------- NQ = 1 IER = 0 CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK) IF ( IER .NE. 0 ) GO TO 685 DO 110 I = 1,N 110 Y(I,2) = H*SAVE1(I) METH = MF/10 MITER = MF - 10*METH L = 2 IDOUB = 3 RMAX = 1.E+04 RC = 0. CRATE = 1. EPSOLD = EPS HOLD = H MFOLD = MF NOLD = N NSTEP = 0 NSTEPJ = 0 NFE = 0 NJE = 1 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET C THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET. C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER. C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES. C IF THE CALLER HAS CHANGED H, Y MUST BE RESCALED. C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT C FURTHER CHANGES IN H FOR THAT MANY STEPS. C----------------------------------------------------------------------- 120 IF (MF .EQ. MFOLD) GO TO 150 MEO = METH MIO = MITER METH = MF/10 MITER = MF - 10*METH MFOLD = MF IF (MITER .NE. MIO) IWEVAL = MITER IF (METH .EQ. MEO) GO TO 150 IDOUB = L + 1 IRET = 1 130 CALL COSET (METH, NQ, EL, TQ) LMAX = MAXDER + 1 RC = RC*EL(1)/OLDL0 OLDL0 = EL(1) 140 FN = FLOAT(N) EDN = FN*(TQ(1)*EPS)**2 E = FN*(TQ(2)*EPS)**2 EUP = FN*(TQ(3)*EPS)**2 BND = FN*(TQ(4)*EPS)**2 GO TO (160, 170, 200), IRET 150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160 EPSOLD = EPS NOLD = N IRET = 1 GO TO 140 160 IF (H .EQ. HOLD) GO TO 200 RH = H/HOLD H = HOLD IREDO = 3 GO TO 175 170 RH = AMAX1(RH,HMIN/ ABS(H)) 175 RH = AMIN1(RH,HMAX/ ABS(H),RMAX) R1 = 1. DO 180 J = 2,L R1 = R1*RH DO 180 I = 1,N 180 Y(I,J) = Y(I,J)*R1 H = H*RH RC = RC*RH IDOUB = L + 1 IF (IREDO .EQ. 0) GO TO 690 C----------------------------------------------------------------------- C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY C MULTIPLYING THE Y ARRAY BY THE PASCAL TRIANGLE MATRIX. C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE PW TO BE UPDATED. C IN ANY CASE, PW IS UPDATED AT LEAST EVERY 40-TH STEP. C PW IS THE CHORD ITERATION MATRIX A - H*EL(1)*(DG/DY). C----------------------------------------------------------------------- 200 IF ( ABS(RC-1.) .GT. 0.3) IWEVAL = MITER IF (NSTEP .GE. NSTEPJ+40) IWEVAL = MITER T = T + H DO 210 J1 = 1,NQ DO 210 J2 = J1,NQ J = (NQ + J1) - J2 DO 210 I = 1,N 210 Y(I,J) = Y(I,J) + Y(I,J+1) C----------------------------------------------------------------------- C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH C IS DEPENDENT ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED C IN THE VECTOR ERROR(I). THE Y ARRAY IS NOT ALTERED IN THE CORRECTOR C LOOP. THE UPDATED Y VECTOR IS STORED TEMPORARILY IN SAVE1. C THE UPDATED H*YDOT IS STORED IN SAVE2. C----------------------------------------------------------------------- 220 DO 230 I = 1,N SAVE2(I) = Y(I,2) 230 ERROR(I) = 0. M = 0 CALL RES (T, H, Y, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK, * WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9)) NFE = NFE + 1 IF (IWEVAL .LE. 0) GO TO 350 C----------------------------------------------------------------------- C IF INDICATED, THE MATRIX PW IS REEVALUATED BEFORE STARTING THE C CORRECTOR ITERATION. IWEVAL IS SET TO 0 AS AN INDICATOR C THAT THIS HAS BEEN DONE. PW IS COMPUTED AND PROCESSED IN PSETIB. C----------------------------------------------------------------------- IWEVAL = 0 RC = 1. NJE = NJE + 1 NSTEPJ = NSTEP CON = -H*EL(1) CALL PSETIB (Y, PW, N0, CON, MITER, IER, WORK(IW1), IWORK, * WORK(IW3),WORK(IW9),SAVE2,IPIV,YMAX,WORK(IW11),WORK(IW12), * WORK(IW13),WORK(IW16),WORK(IW14),WORK(IW15),WORK,NPDE) IF (IER .NE. 0) GO TO 420 C----------------------------------------------------------------------- C COMPUTE THE CORRECTOR ERROR, R SUB M, AND SOLVE THE LINEAR SYSTEM C WITH THAT AS RIGHT-HAND SIDE AND PW AS COEFFICIENT MATRIX, C USING THE LU DECOMPOSITION OF PW. C----------------------------------------------------------------------- 350 CALL SOLB (N0, N, ML, MU, PW, SAVE3, IPIV) 370 D = 0. DO 380 I = 1,N ERROR(I) = ERROR(I) + SAVE3(I) D = D + (SAVE3(I)/YMAX(I))**2 SAVE1(I) = Y(I,1) + EL(1)*ERROR(I) 380 SAVE2(I) = Y(I,2) + ERROR(I) C----------------------------------------------------------------------- C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. C----------------------------------------------------------------------- 400 IF (M .NE. 0) CRATE = AMAX1(.9*CRATE,D/D1) IF ((D*AMIN1(1.E0,2.*CRATE)) .LE. BND) GO TO 450 D1 = D M = M + 1 IF (M .EQ. 3) GO TO 410 CALL RES(T, H, SAVE1, SAVE2, SAVE3, NPDE, NCPTS, WORK(IW1), IWORK, * WORK, WORK(IW14), WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9)) GO TO 350 C----------------------------------------------------------------------- C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. C IF THE MATRIX PW IS NOT UP TO DATE, IT IS REEVALUATED FOR THE C NEXT TRY. OTHERWISE THE Y ARRAY IS RETRACTED TO ITS VALUES C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A C NO-CONVERGENCE EXIT IS TAKEN. C----------------------------------------------------------------------- 410 NFE = NFE + 2 IF (IWEVAL .EQ. -1) GO TO 440 420 T = TOLD RMAX = 2. DO 430 J1 = 1,NQ DO 430 J2 = J1,NQ J = (NQ + J1) - J2 DO 430 I = 1,N 430 Y(I,J) = Y(I,J) - Y(I,J+1) IF ( ABS(H) .LE. HMIN*1.00001) GO TO 680 RH = .25 IREDO = 1 GO TO 170 440 IWEVAL = MITER GO TO 220 C----------------------------------------------------------------------- C THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -1 TO SIGNAL C THAT PW MAY NEED UPDATING ON SUBSEQUENT STEPS. THE ERROR TEST C IS MADE AND CONTROL PASSES TO STATEMENT 500 IF IT FAILS. C----------------------------------------------------------------------- 450 IWEVAL = -1 NFE = NFE + M D = 0. DO 460 I = 1,N 460 D = D + (ERROR(I)/YMAX(I))**2 IF (D .GT. E) GO TO 500 C----------------------------------------------------------------------- C AFTER A SUCCESSFUL STEP, UPDATE THE Y ARRAY. C CONSIDER CHANGING H IF IDOUB = 1. OTHERWISE DECREASE IDOUB BY 1. C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A C FACTOR OF AT LEAST 1.1. IF NOT, IDOUB IS SET TO 10 TO PREVENT C TESTING FOR THAT MANY STEPS. C----------------------------------------------------------------------- KFLAG = 0 IREDO = 0 NSTEP = NSTEP + 1 HUSED = H NQUSED = NQ DO 470 J = 1,L DO 470 I = 1,N 470 Y(I,J) = Y(I,J) + EL(J)*ERROR(I) IF (IDOUB .EQ. 1) GO TO 520 IDOUB = IDOUB - 1 IF (IDOUB .GT. 1) GO TO 700 IF (L .EQ. LMAX) GO TO 700 DO 490 I = 1,N 490 Y(I,LMAX) = ERROR(I) GO TO 700 C----------------------------------------------------------------------- C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. C RESTORE T AND THE Y ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR C ONE LOWER ORDER. C----------------------------------------------------------------------- 500 KFLAG = KFLAG - 1 T = TOLD DO 510 J1 = 1,NQ DO 510 J2 = J1,NQ J = (NQ + J1) - J2 DO 510 I = 1,N 510 Y(I,J) = Y(I,J) - Y(I,J+1) RMAX = 2. IF ( ABS(H) .LE. HMIN*1.00001) GO TO 660 IF (KFLAG .LE. -3) GO TO 640 IREDO = 2 PR3 = 1.E+20 GO TO 540 C----------------------------------------------------------------------- C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE. C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE C ADDITIONAL SCALED DERIVATIVE. C----------------------------------------------------------------------- 520 PR3 = 1.E+20 IF (L .EQ. LMAX) GO TO 540 D1 = 0. DO 530 I = 1,N 530 D1 = D1 + ((ERROR(I) - Y(I,LMAX))/YMAX(I))**2 ENQ3 = .5/ FLOAT(L+1) PR3 = ((D1/EUP)**ENQ3)*1.4 + 1.4E-06 540 ENQ2 = .5/ FLOAT(L) PR2 = ((D/E)**ENQ2)*1.2 + 1.2E-06 PR1 = 1.E+20 IF (NQ .EQ. 1) GO TO 560 D = 0. DO 550 I = 1,N 550 D = D + (Y(I,L)/YMAX(I))**2 ENQ1 = .5/ FLOAT(NQ) PR1 = ((D/EDN)**ENQ1)*1.3 + 1.3E-06 560 IF (PR2 .LE. PR3) GO TO 570 IF (PR3 .LT. PR1) GO TO 590 GO TO 580 570 IF (PR2 .GT. PR1) GO TO 580 NEWQ = NQ RH = 1./PR2 GO TO 620 580 NEWQ = NQ - 1 RH = 1./PR1 GO TO 620 590 NEWQ = L RH = 1./PR3 IF (RH .LT. 1.1) GO TO 610 DO 600 I = 1,N 600 Y(I,NEWQ+1) = ERROR(I)*EL(L)/ FLOAT(L) GO TO 630 610 IDOUB = 10 GO TO 700 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1)) GO TO 610 C----------------------------------------------------------------------- C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. C IN ANY CASE H IS RESET ACCORDING TO RH AND THE Y ARRAY IS RESCALED. C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. C----------------------------------------------------------------------- IF (NEWQ .EQ. NQ) GO TO 170 630 NQ = NEWQ L = NQ + 1 IRET = 2 GO TO 130 C----------------------------------------------------------------------- C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE C Y ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED. C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2. C----------------------------------------------------------------------- 640 IF (KFLAG .EQ. -7) GO TO 670 RH = .1 RH = AMAX1(HMIN/ ABS(H),RH) H = H*RH IER = 0 CALL DIFFUN (N, T, Y, SAVE1, IER, PW, IPIV, WORK, IWORK) IF ( IER .NE. 0 ) GO TO 685 NJE = NJE + 1 DO 650 I = 1,N 650 Y(I,2) = H*SAVE1(I) IWEVAL = MITER IDOUB = 10 IF (NQ .EQ. 1) GO TO 200 NQ = 1 L = 2 IRET = 3 GO TO 130 C----------------------------------------------------------------------- C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. C----------------------------------------------------------------------- 660 KFLAG = -1 GO TO 700 670 KFLAG = -2 GO TO 700 680 KFLAG = -3 GO TO 700 685 KFLAG = -4 GO TO 700 690 RMAX = 10. 700 HOLD = H JSTART = NQ RETURN END SUBROUTINE GFUN (T,C,UDOT,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT, * XC,UVAL,ILEFT) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL. C C SUBROUTINE GFUN COMPUTES THE FUNCTION UDOT=G(C,T), THE RIGHT- C HAND SIDE OF THE SEMI-DISCRETE APPROXIMATION TO THE ORIGINAL C SYSTEM OF PARTIAL DIFFERENTIAL EQUATIONS AND UPDATES THE BOUNDARY C CONDITION INFORMATION. C C PACKAGE ROUTINES CALLED.. EVAL C USER ROUTINES CALLED.. BNDRY,F C CALLED BY.. DIFFUN,PSETIB,RES C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPT,NEQN,IQUAD DIMENSION C(NPDE,NCPTS),UDOT(NPDE,NCPTS) DIMENSION A(NCPTS*3*KORD),BC(NPDE,NPDE,4), * XC(NCPTS),UVAL(NPDE,3),ILEFT(NCPTS) DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE) DO 10 K=1,4 DO 10 J=1,NPDE DO 10 I=1,NPDE BC(I,J,K) = 0.0 10 CONTINUE C----------------------------------------------------------------------- C UPDATE THE LEFT BOUNDARY VALUES. SAVE LEFT BOUNDARY CONDITION C INFORMATION IN THE FIRST 2*NPDE*NPDE LOCATIONS OF BC. C C NOTE.. UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K). C----------------------------------------------------------------------- CALL EVAL(1,NPDE,C,UVAL,A,ILEFT) CALL BNDRY(T,XC(1),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE) CALL F(T,XC(1),UVAL,UVAL(1,2),UVAL(1,3),UDOT,NPDE) ILIM = KORD + 2 DO 30 K=1,NPDE BC(K,K,1) = 1. IF( DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0 ) GO TO 30 UDOT(K,1) = DZDT(K) DO 20 J=1,NPDE BC(K,J,2) = A(ILIM) * DBDUX(K,J) BC(K,J,1) = DBDU(K,J) - BC(K,J,2) 20 CONTINUE 30 CONTINUE C----------------------------------------------------------------------- C MAIN LOOP TO FORM RIGHT SIDE OF ODES AT THE COLLOCATION POINTS. C----------------------------------------------------------------------- ILIM = NCPTS - 1 DO 40 I=2,ILIM CALL EVAL(I,NPDE,C,UVAL,A,ILEFT) CALL F(T,XC(I),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,I),NPDE) 40 CONTINUE C----------------------------------------------------------------------- C UPDATE THE RIGHT BOUNDARY VALUES. SAVE THE RIGHT BOUNDARY CONDITION C INFORMATION IN THE LAST 2*NPDE*NPDE LOCATIONS IN BC. C----------------------------------------------------------------------- CALL EVAL(NCPTS,NPDE,C,UVAL,A,ILEFT) CALL F(T,XC(NCPTS),UVAL,UVAL(1,2),UVAL(1,3),UDOT(1,NCPTS),NPDE) CALL BNDRY(T,XC(NCPTS),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE) ILIM = NCPTS * 3 * KORD - KORD - 1 DO 60 K=1,NPDE BC(K,K,4) = 1. IF( DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0 ) GO TO 60 UDOT(K,NCPTS) = DZDT(K) DO 50 J=1,NPDE BC(K,J,3) = A(ILIM) * DBDUX(K,J) BC(K,J,4) = DBDU(K,J) - BC(K,J,3) 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE EVAL(ICPT,NPDE,C,UVAL,A,ILEFT) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL. C C SUBROUTINE EVAL EVALUATES U(K), UX(K), AND UXX(K), K=1 TO NPDE, C AT THE COLLOCATION POINT WITH INDEX ICPT USING THE VALUES OF C THE BASIS FUNCTION COEFFICIENTS IN C AND THE BASIS FUNCTION VALUES C STORED IN A. THE RESULTS ARE STORED IN UVAL AS FOLLOWS.. C UVAL(K,1) = U(K), UVAL(K,2) = UX(K), AND UVAL(K,3) = UXX(K). C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. GFUN,PDECOL,PSETIB C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD DIMENSION C(NPDE,NCPTS),UVAL(NPDE,3),A(NCPTS*3*KORD),ILEFT(NCPTS) IK = ILEFT(ICPT) - KORD IC = 3*KORD*(ICPT-1) DO 10 M=1,3 ICC = IC + KORD*(M-1) DO 10 J=1,NPDE UVAL(J,M) = 0. DO 10 I=1,KORD UVAL(J,M) = UVAL(J,M) + C(J,I+IK)*A(I+ICC) 10 CONTINUE RETURN END SUBROUTINE DIFFUN (N, T, Y, YDOT, IER, PW, IPIV, WORK, IWORK) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL. C C THIS ROUTINE COMPUTES YDOT = A(Y,T)**-1 * G(Y,T) BY USE OF C THE ROUTINES GFUN, ADDA, DECB, AND SOLB. C C PACKAGE ROUTINES CALLED.. ADDA,DECB,GFUN,SOLB C USER ROUTINES CALLED.. NONE C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /ISTART/ IW1,IW2,IW3,IDUM(5),IW9,IW10,IW11,IW12,IW13,IW14, * IW15,IW16,IW17,IW18 COMMON /OPTION/ NOGAUS,MAXDER DIMENSION Y(NEQN),YDOT(NEQN),PW(NEQN*(3*ML+1)), * IPIV(NEQN),IWORK((NPDE+1)*NCPTS), * WORK(KORD+NPDE*(4+9*NPDE)+(KORD+(NINT-1)*(KORD-NCC))* * (3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4))) CALL GFUN (T, Y, YDOT, NPDE, NCPTS, WORK(IW1), WORK, WORK(IW14), * WORK(IW15), WORK(IW16), WORK(IW3), WORK(IW9), IWORK) DO 10 I = 1,N0W 10 PW(I) = 0. N0 = NM1 + 1 CALL ADDA (PW, N0, WORK(IW1), IWORK, WORK, NPDE) CALL DECB (N0, N, ML, MU, PW, IPIV, IER) IF ( IER .NE. 0 ) RETURN CALL SOLB (N0, N, ML, MU, PW, YDOT, IPIV) RETURN END SUBROUTINE ADDA(PW,N0,A,ILEFT,BC,NPDE) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL AND STIFIB. C C SUBROUTINE ADDA ADDS THE MATRIX A TO THE MATRIX STORED IN PW IN C BAND FORM. PW IS STORED BY DIAGONALS WITH THE LOWERMOST DIAGONAL C STORED IN THE FIRST COLUMN OF THE ARRAY. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. DIFFUN,PSETIB C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W DIMENSION PW(NEQN,3*ML+1),A(3*KORD*NCPTS), * ILEFT(NCPTS),BC(NPDE,NPDE,4) C----------------------------------------------------------------------- C ADD THE BOUNDARY CONDITION PORTIONS OF THE A MATRIX TO PW ( THE FIRST C AND LAST BLOCK ROWS). C----------------------------------------------------------------------- ICOL = (ILEFT(1) + IQUAD - 1) * NPDE DO 10 I=1,NPDE IBOT = NEQN - NPDE + I DO 10 J=1,NPDE IND = ICOL + J - I PW(I,IND) = PW(I,IND) + BC(I,J,1) PW(I,IND+NPDE) = PW(I,IND+NPDE) + BC(I,J,2) PW(IBOT,IND-NPDE) = PW(IBOT,IND-NPDE) + BC(I,J,3) PW(IBOT,IND) = PW(IBOT,IND) + BC(I,J,4) 10 CONTINUE C----------------------------------------------------------------------- C UPDATE THE REMAINING ROWS OF PW BY ADDING THE APPROPRIATE VALUES C IN A TO PW. C----------------------------------------------------------------------- IND = NCPTS - 1 DO 20 I=2,IND I1 = (I-1) * NPDE I2 = (I-1) * KORD * 3 ICOL = ILEFT(I) - I + IQUAD - 1 DO 20 J=1,KORD J1 = (ICOL+J) * NPDE J2 = I2 + J DO 20 JJ=1,NPDE 20 PW(I1+JJ,J1) = PW(I1+JJ,J1) + A(J2) RETURN END SUBROUTINE RES(T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT, * XC,UVAL) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL. C C SUBROUTINE RES COMPUTES THE RESIDUAL VECTOR R = H*G(C,T) - A(C,T)*V C WHERE H IS THE CURRENT TIME STEP SIZE, G IS A VECTOR, A IS A C MATRIX, V IS A VECTOR, AND T IS THE CURRENT TIME. C C PACKAGE ROUTINES CALLED.. GFUN C USER ROUTINES CALLED.. NONE C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- SAVE COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPT,NEQN,IQUAD DIMENSION C(NPDE,NCPTS),R(NPDE,NCPTS),V(NPDE,NCPTS) DIMENSION A(3*KORD*NCPTS),ILEFT(NCPTS),BC(NPDE,NPDE,4),XC(NCPTS), * UVAL(3*NPDE) DIMENSION DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE),DZDT(NPDE) C----------------------------------------------------------------------- C FORM G(C,T) AND STORE IN R. C----------------------------------------------------------------------- CALL GFUN(T,C,R,NPDE,NCPTS,A,BC,DBDU,DBDUX,DZDT,XC,UVAL,ILEFT) C----------------------------------------------------------------------- C FORM THE FIRST AND LAST BLOCK ROWS OF THE RESIDUAL VECTOR C WHICH ARE DEPENDENT ON THE BOUNDARY CONDITIONS. C----------------------------------------------------------------------- ILIM = NCPTS - 1 DO 20 I=1,NPDE SUM1 = 0.0 SUM2 = 0.0 DO 10 J=1,NPDE SUM1 = SUM1 + BC(I,J,1) * V(J,1) + BC(I,J,2) * V(J,2) SUM2 = SUM2 + BC(I,J,3) * V(J,ILIM) + BC(I,J,4) * V(J,NCPTS) 10 CONTINUE R(I,1) = H * R(I,1) - SUM1 R(I,NCPTS) = H * R(I,NCPTS) - SUM2 20 CONTINUE C----------------------------------------------------------------------- C FORM THE REMAINING COMPONENTS OF THE RESIDUAL VECTOR. C----------------------------------------------------------------------- DO 50 ICPTS=2,ILIM I2 = (ICPTS-1) * KORD * 3 ICOL = ILEFT(ICPTS) - KORD DO 40 JJ=1,NPDE SUM1 = 0. DO 30 J=1,KORD SUM1 = SUM1 + A(I2+J) * V(JJ,ICOL+J) 30 CONTINUE R(JJ,ICPTS) = H*R(JJ,ICPTS) - SUM1 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE PSETIB (C, PW, N0, CON, MITER, IER, A, ILEFT, XC, UVAL, * SAVE2,IPIV,CMAX,DFDU,DFDUX,DFDUXX,DZDT,DBDU,DBDUX,BC,NPDE) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL AND STIFIB. C C PSETIB IS CALLED BY STIFIB TO COMPUTE AND PROCESS THE MATRIX C PW = A - H*EL(1)*(DG/DC), WHERE A AND DG/DC ARE TREATED IN BAND C FORM. DG/DC IS COMPUTED, EITHER WITH THE AID OF THE USER-SUPPLIED C ROUTINE DERIVF IF MITER = 1, OR BY FINITE DIFFERENCING WITH THE AID C OF THE PACKAGE-SUPPLIED ROUTINE DIFFF IF MITER = 2. FINALLY, C PW IS SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER C SOLUTION OF LINEAR SYSTEMS WITH PW AS COEFFICIENT MATRIX. C SEE SUBROUTINES DECB AND SOLB. C C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION C WITH PSETIB USES THE FOLLOWING.. C EPSJ = SQRT(UROUND), USED IN THE NUMERICAL JACOBIAN INCREMENTS. C MW = ML + MU + 1. C NM1 = N0 - 1. C N0ML = N0*ML. C N0W = N0*MW. C C PACKAGE ROUTINES CALLED.. ADDA,DECB,DIFFF,EVAL,GFUN C USER ROUTINES CALLED.. BNDRY,DERIVF C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. ABS,FLOAT,MAX0,MIN0,SQRT C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD COMMON /GEAR1/ T,H,DUMMY(3),UROUND,N,IDUMMY(3) COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W DIMENSION PW(NEQN,3*ML+1),C(NEQN),CMAX(NEQN) DIMENSION A(3*KORD*NCPTS),ILEFT(NCPTS),BC(4*NPDE*NPDE), * XC(NCPTS),UVAL(NPDE,3),SAVE2(NEQN),IPIV(NEQN) DIMENSION DFDU(NPDE,NPDE),DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE) DIMENSION DZDT(NPDE),DBDU(NPDE,NPDE),DBDUX(NPDE,NPDE) DO 10 I=1,NEQN DO 5 J=1,MW 5 PW(I,J)=0.0E0 10 CONTINUE IF ( MITER .EQ. 1 ) GO TO 25 CALL GFUN (T, C, SAVE2, NPDE, NCPTS,A,BC,DBDU,DBDUX,DZDT,XC, * UVAL,ILEFT) D = 0. DO 20 I = 1,N 20 D = D + SAVE2(I)**2 R0 = ABS(H)* SQRT(D/FLOAT(N0))*1.E+03*UROUND C----------------------------------------------------------------------- C COMPUTE BLOCK ROWS OF JACOBIAN. C----------------------------------------------------------------------- 25 DO 30 I=1,NCPTS I1 = (I-1)*NPDE I2 = (I-1)*KORD*3 CALL EVAL(I,NPDE,C,UVAL,A,ILEFT) IF ( MITER .EQ. 1 ) * CALL DERIVF(T,XC(I),UVAL,UVAL(1,2),UVAL(1,3), * DFDU,DFDUX,DFDUXX,NPDE) IF ( MITER .EQ. 2 ) * CALL DIFFF(T,XC(I),I,UVAL,UVAL(1,2),UVAL(1,3), * DFDU,DFDUX,DFDUXX,NPDE,CMAX,SAVE2) ICOL = ILEFT(I) - I + IQUAD - 1 KLOW = MAX0(1,I+2-NCPTS) KUP = MIN0(KORD,KORD+I-2) DO 30 KBLK=KLOW,KUP J1 = (ICOL+KBLK)*NPDE J2 = I2 + KBLK J3 = J2 + KORD J4 = J3 + KORD DO 30 L=1,NPDE DO 30 K=1,NPDE PW(I1+K,J1-K+L) = DFDU(K,L)*A(J2) + DFDUX(K,L)*A(J3) * + DFDUXX(K,L)*A(J4) 30 CONTINUE C----------------------------------------------------------------------- C MODIFY THE LAST AND THE FIRST BLOCK ROWS FOR THE BOUNDARY CONDITIONS. C CURRENT INFORMATION FOR THE RIGHT BOUNDARY CONDITION IS ALREADY IN C THE ARRAYS DBDU, DBDUX AS A RESULT OF A PREVIOUS CALL TO GFUN. C----------------------------------------------------------------------- IROW = NEQN - NPDE DO 50 K=1,NPDE IROW = IROW + 1 IF(DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0) GO TO 50 DO 40 J=1,MW PW(IROW,J) = 0.0 40 CONTINUE 50 CONTINUE CALL EVAL(1,NPDE,C,UVAL,A,ILEFT) CALL BNDRY(T,XC(1),UVAL,UVAL(1,2),DBDU,DBDUX,DZDT,NPDE) DO 70 K=1,NPDE IF(DBDU(K,K) .EQ. 0.0 .AND. DBDUX(K,K) .EQ. 0.0) GO TO 70 DO 60 J=1,MW PW(K,J) = 0.0 60 CONTINUE 70 CONTINUE DO 80 I = 1,N0 DO 85 J=1,MW 85 PW(I,J)=PW(I,J)*CON 80 CONTINUE C----------------------------------------------------------------------- C ADD MATRIX A(C,T) TO PW. C----------------------------------------------------------------------- CALL ADDA (PW, N0, A, ILEFT, BC, NPDE) C----------------------------------------------------------------------- C DO LU DECOMPOSITION ON PW. C----------------------------------------------------------------------- CALL DECB (N0, N, ML, MU, PW, IPIV, IER) RETURN END SUBROUTINE DIFFF(T,X,IPT,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE,CMAX, * SAVE2) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN PDECOL. C C SUBROUTINE DIFFF IS USED IF MITER=2 TO PROVIDE FINITE DIFFERENCE C APPROXIMATIONS FOR THE PARTIAL DERIVATIVES OF THE K-TH USER DEFINED C FUNCTION IN THE F ROUTINE WITH RESPECT TO THE VARIABLES U, UX, AND C UXX. THESE PARTIALS WITH RESPECT TO U, UX, AND UXX ARE COMPUTED, C STORED, AND RETURNED IN THE NPDE BY NPDE ARRAYS DFDU, DFDUX, AND C DFDUXX, RESPECTIVELY, AT COLLOCATION POINT NUMBER IPT. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. F C CALLED BY.. PSETIB C FORTRAN FUNCTIONS USED.. AMAX1 C----------------------------------------------------------------------- COMMON /GEAR9/ EPSJ,R0,ML,MU,MW,NM1,N0ML,N0W COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD DIMENSION U(NPDE),UX(NPDE),UXX(NPDE),DFDU(NPDE,NPDE), * DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE),CMAX(NEQN),SAVE2(NEQN) ID = (IPT-1) * NPDE DO 40 J=1,NPDE UJ = U(J) R = EPSJ * CMAX(J) R = AMAX1(R,R0) U(J) = U(J) + R RINV = 1. / R CALL F(T,X,U,UX,UXX,DFDU(1,J),NPDE) DO 10 I=1,NPDE 10 DFDU(I,J) = ( DFDU(I,J) - SAVE2(I+ID) ) * RINV U(J) = UJ UJ = UX(J) UX(J) = UX(J) + R CALL F(T,X,U,UX,UXX,DFDUX(1,J),NPDE) DO 20 I=1,NPDE 20 DFDUX(I,J) = ( DFDUX(I,J) - SAVE2(I+ID) ) * RINV UX(J) = UJ UJ = UXX(J) UXX(J) = UXX(J) + R CALL F(T,X,U,UX,UXX,DFDUXX(1,J),NPDE) DO 30 I=1,NPDE 30 DFDUXX(I,J) = ( DFDUXX(I,J) - SAVE2(I+ID) ) * RINV UXX(J) = UJ 40 CONTINUE RETURN END SUBROUTINE INTERP (TOUT, Y, N0, Y0) C----------------------------------------------------------------------- C CALLING ARGUMENTS ARE DEFINED BELOW AND IN STIFIB C C SUBROUTINE INTERP COMPUTES INTERPOLATED VALUES OF THE DEPENDENT C VARIABLE Y AND STORES THEM IN Y0. THE INTERPOLATION IS TO THE C POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY Y, AS FOLLOWS.. C NQ C Y0(I) = SUM Y(I,J+1)*S**J , C J=0 C WHERE S = -(T-TOUT)/H. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. PDECOL C FORTRAN FUNCTIONS USED.. NONE C----------------------------------------------------------------------- COMMON /SIZES/ NINT,KORD,NCC,NPD,NCPTS,NEQN,IQUAD COMMON /OPTION/ NOGAUS,MAXDER COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART DIMENSION Y0(NEQN),Y(NEQN,MAXDER+1) DO 10 I = 1,N 10 Y0(I) = Y(I,1) L = JSTART + 1 S = (TOUT - T)/H S1 = 1. DO 30 J = 2,L S1 = S1*S DO 20 I = 1,N 20 Y0(I) = Y0(I) + S1*Y(I,J) 30 CONTINUE RETURN END SUBROUTINE COSET (METH, NQ, EL, TQ) C----------------------------------------------------------------------- C COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE. C THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD. C THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE C IN RELATION TO TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE C PERTST ARRAY. C C THE VECTORS EL AND TQ DEPEND ON METH AND NQ. C THE MAXIMUM ORDER, MAXDER, OF THE METHODS AVAILABLE IS CURRENTLY C 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS. MAXDER DEFAULTS C TO 5 UNLESS THE USER SETS MAXDER TO SOME OTHER LEGITIMATE VALUE C THROUGH THE COMMON BLOCK /OPTION/. SEE PDECOL FOR ADDITIONAL DETAILS. C LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE Y ARRAY (SEE STIFIB C AND THE VARIABLE C, Y, OR WORK(IW10) IN PDECOL. C C THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT C ONE PERCENT ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW C IS.. COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ, C COEFFICIENTS FOR ORDER NQ + 1. WITHIN EACH GROUP ARE THE C COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE C BDF METHODS. C C REFERENCE C C GEAR, C.W., NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY C DIFFERENTIAL EQUATIONS, PRENTICE-HALL, ENGLEWOOD CLIFFS, C N. J., 1971. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. STIFIB C FORTRAN FUNCTIONS USED.. FLOAT C----------------------------------------------------------------------- DIMENSION PERTST(12,2,3),EL(13),TQ(4) DATA PERTST / 1.,1.,2.,1.,.3158,.07407,.01391,.002182, * .0002945,.00003492,.000003692,.0000003524, * 1.,1.,.5,.1667,.04167,1.,1.,1.,1.,1.,1.,1., * 2.,12.,24.,37.89,53.33,70.08,87.97,106.9, * 126.7,147.4,168.8,191.0, * 2.0,4.5,7.333,10.42,13.7,1.,1.,1.,1.,1.,1.,1., * 12.0,24.0,37.89,53.33,70.08,87.97,106.9, * 126.7,147.4,168.8,191.0,1., * 3.0,6.0,9.167,12.5,1.,1.,1.,1.,1.,1.,1.,1. / C GO TO (1,2),METH 1 GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ 2 GO TO (201,202,203,204,205),NQ C----------------------------------------------------------------------- C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY. C FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE C GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I).. C L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ. C FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY C DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K, L(-1) = 0, C WHERE K = FACTORIAL(NQ-1). C FOR THE BDF METHODS, C L(T) = (T+1)*(T+2)* ... *(T+NQ)/K, C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ). C C THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS.. C IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12, C BDF METHODS OF ORDERS 1 TO 5. C----------------------------------------------------------------------- 101 EL(1) = 1.0E-00 GO TO 900 102 EL(1) = 0.5E-00 EL(3) = 0.5E-00 GO TO 900 103 EL(1) = 4.1666666666667E-01 EL(3) = 0.75E-00 EL(4) = 1.6666666666667E-01 GO TO 900 104 EL(1) = 0.375E-00 EL(3) = 9.1666666666667E-01 EL(4) = 3.3333333333333E-01 EL(5) = 4.1666666666667E-02 GO TO 900 105 EL(1) = 3.4861111111111E-01 EL(3) = 1.0416666666667E-00 EL(4) = 4.8611111111111E-01 EL(5) = 1.0416666666667E-01 EL(6) = 8.3333333333333E-03 GO TO 900 106 EL(1) = 3.2986111111111E-01 EL(3) = 1.1416666666667E-00 EL(4) = 0.625E-00 EL(5) = 1.7708333333333E-01 EL(6) = 0.025E-00 EL(7) = 1.3888888888889E-03 GO TO 900 107 EL(1) = 3.1559193121693E-01 EL(3) = 1.225E-00 EL(4) = 7.5185185185185E-01 EL(5) = 2.5520833333333E-01 EL(6) = 4.8611111111111E-02 EL(7) = 4.8611111111111E-03 EL(8) = 1.9841269841270E-04 GO TO 900 108 EL(1) = 3.0422453703704E-01 EL(3) = 1.2964285714286E-00 EL(4) = 8.6851851851852E-01 EL(5) = 3.3576388888889E-01 EL(6) = 7.7777777777778E-02 EL(7) = 1.0648148148148E-02 EL(8) = 7.9365079365079E-04 EL(9) = 2.4801587301587E-05 GO TO 900 109 EL(1) = 2.9486800044092E-01 EL(3) = 1.3589285714286E-00 EL(4) = 9.7655423280423E-01 EL(5) = 0.4171875E-00 EL(6) = 1.1135416666667E-01 EL(7) = 0.01875E-00 EL(8) = 1.9345238095238E-03 EL(9) = 1.1160714285714E-04 EL(10)= 2.7557319223986E-06 GO TO 900 110 EL(1) = 2.8697544642857E-01 EL(3) = 1.4144841269841E-00 EL(4) = 1.0772156084656E-00 EL(5) = 4.9856701940035E-01 EL(6) = 0.1484375E-00 EL(7) = 2.9060570987654E-02 EL(8) = 3.7202380952381E-03 EL(9) = 2.9968584656085E-04 EL(10)= 1.3778659611993E-05 EL(11)= 2.7557319223986E-07 GO TO 900 111 EL(1) = 2.8018959644394E-01 EL(3) = 1.4644841269841E-00 EL(4) = 1.1715145502646E-00 EL(5) = 5.7935819003527E-01 EL(6) = 1.8832286155203E-01 EL(7) = 4.1430362654321E-02 EL(8) = 6.2111441798942E-03 EL(9) = 6.2520667989418E-04 EL(10)= 4.0417401528513E-05 EL(11)= 1.5156525573192E-06 EL(12)= 2.5052108385442E-08 GO TO 900 112 EL(1) = 2.7426554003160E-01 EL(3) = 1.5099386724387E-00 EL(4) = 1.2602711640212E-00 EL(5) = 6.5923418209877E-01 EL(6) = 2.3045800264550E-01 EL(7) = 5.5697246105232E-02 EL(8) = 9.4394841269841E-03 EL(9) = 1.1192749669312E-03 EL(10)= 9.0939153439153E-05 EL(11)= 4.8225308641975E-06 EL(12)= 1.5031265031265E-07 EL(13)= 2.0876756987868E-09 GO TO 900 201 EL(1) = 1.0E-00 GO TO 900 202 EL(1) = 6.6666666666667E-01 EL(3) = 3.3333333333333E-01 GO TO 900 203 EL(1) = 5.4545454545455E-01 EL(3) = EL(1) EL(4) = 9.0909090909091E-02 GO TO 900 204 EL(1) = 0.48E-00 EL(3) = 0.7E-00 EL(4) = 0.2E-00 EL(5) = 0.02E-00 GO TO 900 205 EL(1) = 4.3795620437956E-01 EL(3) = 8.2116788321168E-01 EL(4) = 3.1021897810219E-01 EL(5) = 5.4744525547445E-02 EL(6) = 3.6496350364964E-03 C 900 DO 910 K = 1,3 910 TQ(K) = PERTST(NQ,METH,K) TQ(4) = .5E-00*TQ(2)/ FLOAT(NQ+2) RETURN END SUBROUTINE DECB (NDIM, N, ML, MU, B, IPIV, IER) C----------------------------------------------------------------------- C SUBROUTINES DECB AND SOLB FORM A TWO SUBROUTINE PACKAGE FOR THE C DIRECT SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IN WHICH THE C COEFFICIENT MATRIX IS REAL AND BANDED. C C LU DECOMPOSITION OF BAND MATRIX A.. L*U = P*A , WHERE P IS A C PERMUTATION MATRIX, L IS A UNIT LOWER TRIANGULAR MATRIX, C AND U IS AN UPPER TRIANGULAR MATRIX. C N = ORDER OF MATRIX. C B = N BY (2*ML+MU+1) ARRAY CONTAINING THE MATRIX A ON INPUT C AND ITS FACTORED FORM ON OUTPUT. C ON INPUT, B(I,K) (1.LE.I.LE.N) CONTAINS THE K-TH C DIAGONAL OF A, OR A(I,J) IS STORED IN B(I,J-I+ML+1). C ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH C U IN COLUMNS 1 TO ML+MU+1, AND L IN COLUMNS C ML+MU+2 TO 2*ML+MU+1. C ML,MU = WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, NOT C COUNTING THE MAIN DIAGONAL. TOTAL BANDWIDTH IS ML+MU+1. C NDIM = THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B. C NDIM MUST BE .GE. N. C IPIV = ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION. C IER = ERROR INDICATOR.. C = 0 IF NO ERROR, C = K IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR). C THE INPUT ARGUMENTS ARE NDIM, N, ML, MU, B. C THE OUTPUT ARGUMENTS ARE B, IPIV, IER. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. DIFFUN,INITAL,PSETIB C FORTRAN FUNCTIONS USED.. ABS,MIN0 C----------------------------------------------------------------------- DIMENSION B(NDIM,2*ML+MU+1),IPIV(N) IER = 0 IF (N .EQ. 1) GO TO 92 LL = ML + MU + 1 N1 = N - 1 IF (ML .EQ. 0) GO TO 32 DO 30 I = 1,ML II = MU + I K = ML + 1 - I DO 10 J = 1,II 10 B(I,J) = B(I,J+K) K = II + 1 DO 20 J = K,LL 20 B(I,J) = 0. 30 CONTINUE 32 LR = ML DO 90 NR = 1,N1 NP = NR + 1 IF (LR .NE. N) LR = LR + 1 MX = NR XM = ABS(B(NR,1)) IF (ML .EQ. 0) GO TO 42 DO 40 I = NP,LR IF ( ABS(B(I,1)) .LE. XM) GO TO 40 MX = I XM = ABS(B(I,1)) 40 CONTINUE 42 IPIV(NR) = MX IF (MX .EQ. NR) GO TO 60 DO 50 I = 1,LL XX = B(NR,I) B(NR,I) = B(MX,I) 50 B(MX,I) = XX 60 XM = B(NR,1) IF (XM .EQ. 0.) GO TO 100 B(NR,1) = 1./XM IF (ML .EQ. 0) GO TO 90 XM = -B(NR,1) KK = MIN0(N-NR,LL-1) DO 80 I = NP,LR J = LL + I - NR XX = B(I,1)*XM B(NR,J) = XX DO 70 II = 1,KK 70 B(I,II) = B(I,II+1) + XX*B(NR,II+1) 80 B(I,LL) = 0. 90 CONTINUE 92 NR = N IF (B(N,1) .EQ. 0.) GO TO 100 B(N,1) = 1./B(N,1) RETURN 100 IER = NR RETURN END SUBROUTINE SOLB (NDIM, N, ML, MU, B, Y, IPIV) C----------------------------------------------------------------------- C SUBROUTINES DECB AND SOLB FORM A TWO SUBROUTINE PACKAGE FOR THE C DIRECT SOLUTION OF A SYSTEM OF LINEAR EQUATIONS IN WHICH THE C COEFFICIENT MATRIX IS REAL AND BANDED. C C SOLUTION OF A*X = C GIVEN LU DECOMPOSITION OF A FROM DECB. C Y = RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT, C = SOLUTION VECTOR X ON OUTPUT. C ALL THE ARGUMENTS ARE INPUT ARGUMENTS. C THE OUTPUT ARGUMENT IS Y. C C PACKAGE ROUTINES CALLED.. NONE C USER ROUTINES CALLED.. NONE C CALLED BY.. DIFFUN,INITAL,STIFIB C FORTRAN FUNCTIONS USED.. MIN0 C----------------------------------------------------------------------- DIMENSION B(NDIM,2*ML+MU+1),Y(N),IPIV(N) IF (N .EQ. 1) GO TO 60 N1 = N - 1 LL = ML + MU + 1 IF (ML .EQ. 0) GO TO 32 DO 30 NR = 1,N1 IF (IPIV(NR) .EQ. NR) GO TO 10 J = IPIV(NR) XX = Y(NR) Y(NR) = Y(J) Y(J) = XX 10 KK = MIN0(N-NR,ML) DO 20 I = 1,KK 20 Y(NR+I) = Y(NR+I) + Y(NR)*B(NR,LL+I) 30 CONTINUE 32 LL = LL - 1 Y(N) = Y(N)*B(N,1) KK = 0 DO 50 NB = 1,N1 NR = N - NB IF (KK .NE. LL) KK = KK + 1 DP = 0. IF (LL .EQ. 0) GO TO 50 DO 40 I = 1,KK 40 DP = DP + B(NR,I+1)*Y(NR+I) 50 Y(NR) = (Y(NR) - DP)*B(NR,1) RETURN 60 Y(1) = Y(1)*B(1,1) RETURN END C C*** Example Driver 1 C C PROGRAM EXAMPL(UGOUT,TAPE3=UGOUT) INTEGER IOUT PARAMETER (IOUT = 6) COMMON /ENDPT/ XLEFT COMMON /GEAR0/ DTUSED,NQ,NSTEPS,NF,NJ DIMENSION U(2,31),XBKPT(31),SCTCH(10),WORK(5000),IWORK(500) C C INITIALIZE PARAMETERS AND PDECOL CALLING ARGUMENTS C NPDE = 2 NINT = 30 NPTS = NINT + 1 KORD = 4 NCC = 2 T0 = 0.0 TOUT = 1.E-3 DT = 1.E-7 EPS = 1.E-4 MF = 21 INDEX = 1 IWORK(1) = 5000 IWORK(2) = 500 DX = 1. / FLOAT(NPTS-1) DO 10 I=1,NPTS XBKPT(I) = FLOAT(I-1) * DX 10 CONTINUE XLEFT = XBKPT(1) C C CALL THE PACKAGE TO INTEGRATE TO TIME T = TOUT C 20 CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) C C CHECK FOR EXECUTION ERRORS C IF ( INDEX .NE. 0 ) GO TO 70 C C OUTPUT PERFORMANCE DATA AND COMPUTED SOLUTION VALUES C WRITE(IOUT,30) TOUT,DTUSED,NSTEPS 30 FORMAT(//10X,3HT= ,E10.3,7H DT= ,E10.3,15H TOTAL STEPS=,I5) CALL VALUES(XBKPT,U,SCTCH,NPDE,NPTS,NPTS,0,WORK) DO 60 K=1,NPDE WRITE(IOUT,40) K 40 FORMAT(/10X,13HPDE COMPONENT,I2/) WRITE(IOUT,50) (U(K,I),I=1,NPTS) 50 FORMAT(10X,5E12.4) 60 CONTINUE C C SET NEW OUTPUT TIME AND CONTINUE THE INTEGRATION IF TOUT .LT. 11. C OTHERWISE, TERMINATE THE PROBLEM C TOUT = TOUT * 10. IF ( TOUT .LT. 11. ) GO TO 20 70 WRITE(IOUT,80) INDEX 80 FORMAT(10X,7HINDEX= ,I3) STOP END SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) FVAL(1) = U(2)*U(2)*UXX(1) - U(1)*U(2) - U(1)**2 + 10. * + 2.*U(2)*UX(2)*UX(1) FVAL(2) = U(1)*U(1)*UXX(2) + U(1)*U(2) - U(2)**2 * + UXX(1) + 2.*U(1)*UX(1)*UX(2) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) COMMON /ENDPT/ XLEFT IF( X .NE. XLEFT ) GO TO 10 DBDU(1,1) = 1. DBDU(1,2) = 0. DBDU(2,1) = 0. DBDU(2,2) = 1. DBDUX(1,1) = 0. DBDUX(1,2) = 0. DBDUX(2,1) = 0. DBDUX(2,2) = 0. DZDT(1) = 0. DZDT(2) = 0. RETURN 10 DBDU(1,1) = U(2) * COS( U(1) * U(2) ) DBDU(1,2) = U(1) * COS( U(1) * U(2) ) DBDU(2,1) = U(2) * SIN( U(1) * U(2) ) DBDU(2,2) = U(1) * SIN( U(1) * U(2) ) DBDUX(1,1) = 1. DBDUX(1,2) = 0. DBDUX(2,1) = 0. DBDUX(2,2) = 1. DZDT(1) = 0. DZDT(2) = 0. RETURN END SUBROUTINE UINIT(X,U,NPDE) DIMENSION U(NPDE) C C SET INITIAL CONDITIONS. NOTE THAT PI = 4.*ATAN(1.) C U(1) = .5 * ( X + 1. ) U(2) = 4.*ATAN(1.) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) DFDU(1,1) = -U(2) - 2.*U(1) DFDU(1,2) = 2.*U(2)*UXX(1) - U(1) + 2.*UX(2)*UX(1) DFDU(2,1) = 2.*U(1)*UXX(2) + U(2) + 2.*UX(1)*UX(2) DFDU(2,2) = U(1) - 2.*U(2) C DFDUX(1,1) = 2.*U(2)*UX(2) DFDUX(1,2) = 2.*U(2)*UX(1) DFDUX(2,1) = 2.*U(1)*UX(2) DFDUX(2,2) = 2.*U(1)*UX(1) C DFDUXX(1,1) = U(2)*U(2) DFDUXX(1,2) = 0. DFDUXX(2,1) = 1. DFDUXX(2,2) = U(1)*U(1) RETURN END C C*** Example Driver 2 C C PROGRAM SINE(TSIN,TSTOUT,TAPE2=TSIN,TAPE3=TSTOUT) C THIS PROGRAM WAS WRITTEN TO BE RUN IN SINGLE PRECISION ON A CDC7600. INTEGER IN, IOUT PARAMETER (IN=2, IOUT=6) DIMENSION XPLT(25),UPLT(2,2,404),SCRTCH(10) DIMENSION ERRMX(2),UEXACT(2),XBKPT(404) DIMENSION WORK(30000),IWORK(1000) CHARACTER *80 TITLE COMMON /GEAR0/ HUSED,NQUSED,NS,NF,NJ COMMON /ISTART/ IW1,IW2,IW3,IDUM(15) COMMON /PARAMS/ PI C READ IN AND INITIALIZE PDECOL CALLING ARGUMENTS AND PARAMETERS. PI = 4. * ATAN(1.0) OPEN(IN,FILE='EG2.DAT',STATUS='OLD',FORM='FORMATTED') READ(IN,10) TITLE 10 FORMAT(A) READ(IN,20) NPDE,MF,TOUT TKEEP = TOUT 20 FORMAT(2I3,E10.3) WRITE(IOUT,10) TITLE WRITE(IOUT,30) NPDE,MF,TOUT 30 FORMAT(//8H NPDE = ,I2,5X,5HMF = ,I2, * 5X,26H RESULTS ARE FOR TIME T = ,E10.3// *27H KORD NINT EPS ,24HNSTEPS NRES NJAC NQ , *1X,30H ERROR(1) ERROR(2) ) 40 READ(IN,50) KORD,NINT,EPS 50 FORMAT(2I3,E10.3) IF(KORD.LE.0) STOP NPTS = NINT + 1 NCC = 2 IWORK(1) = 30000 IWORK(2) = 1000 T0 = 0.0 C SET TOUT EACH TIME ROUND SO IF THERE IS A FAILURE THE NEXT C SET OF DATA GOES TO THE RIGHT VALUE INSTAED OF STOPPING EARLY TOUT = TKEEP DT = EPS * .5 INDEX = 1 NCPTS1 = KORD*NINT - NCC*(NINT-1) - 1 DX = 1. / FLOAT(NINT) DO 60 I=1,NPTS XBKPT(I) = FLOAT(I-1)*DX 60 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) C CHECK FOR ERRORS AND THEN COMPUTE AN ESTIMATE FOR THE OVERALL C MAXIMUM ERROR FOR EACH PDE COMPONENT. OUTPUT RESULTS. IF(INDEX.NE.0) WRITE(IOUT,70) INDEX 70 FORMAT(7H INDEX=,I3) DO 80 I=1,NPDE 80 ERRMX(I) = 0. C WORK(IW3) TO WORK(IW3+NCPTS1) CONTAINS THE COLLOCATION POINTS. DO 100 J=1,NCPTS1 DX = (WORK(IW3+J)-WORK(IW3+J-1))/24. DO 90 I=2,24 90 XPLT(I) = FLOAT(I-1)*DX + WORK(IW3+J-1) XPLT(1) = WORK(IW3+J-1) XPLT(25) = WORK(IW3+J) DO 100 I=1,25 CALL VALUES(XPLT(I),UPLT(1,1,I),SCRTCH,1,1,1,0,WORK) CALL TRUSOL(XPLT(I),TOUT,UEXACT) DO 100 K=1,NPDE UPLT(K,2,I) = ABS( UPLT(K,1,I) - UEXACT(K) ) IF(UPLT(K,2,I) .GT. ERRMX(K)) ERRMX(K) = UPLT(K,2,I) 100 CONTINUE WRITE(IOUT,110) KORD,NINT,EPS,NS,NF,NJ,NQUSED,(ERRMX(I),I=1,NPDE) 110 FORMAT(I4,I6,1E14.3,1X,4I6,2E15.3) GO TO 40 END C SUBROUTINES FOR SINE PROBLEM. SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) COMMON /PARAMS/ PI FVAL(1) = UXX(1) + PI*PI*SIN(PI*X) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) DBDU(1,1) = 1. DBDUX(1,1) = 0. DZDT(1) = 0. RETURN END SUBROUTINE UINIT(X,U,NPDE) DIMENSION U(NPDE) CALL TRUSOL(X,0.,U(1)) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) DFDUXX(1,1) = 1. DFDUX(1,1) = 0. DFDU(1,1) = 0. RETURN END SUBROUTINE TRUSOL(X,T,UEXACT) DIMENSION UEXACT(1) COMMON /PARAMS/ PI UEXACT(1) = 1. + SIN(PI*X) * ( 1. - EXP(-PI*PI*T) ) RETURN END C C*** Example Driver 3 C C PROGRAM TWO(TSIN,TSTOUT,TAPE2=TSIN,TAPE3=TSTOUT) C THIS PROGRAM WAS WRITTEN TO BE RUN IN SINGLE PRECISION ON A CDC7600. INTEGER IN, IOUT PARAMETER (IN=2, IOUT=6) DIMENSION XPLT(25),UPLT(2,2,404),SCRTCH(10) DIMENSION ERRMX(2),UEXACT(2),XBKPT(404) DIMENSION WORK(30000),IWORK(1000) CHARACTER *80 TITLE COMMON /GEAR0/ HUSED,NQUSED,NS,NF,NJ COMMON /ISTART/ IW1,IW2,IW3,IDUM(15) COMMON /PARAMS/ PI C READ IN AND INITIALIZE PDECOL CALLING ARGUMENTS AND PARAMETERS. PI = 4. * ATAN(1.0) OPEN(IN,FILE='EG3.DAT',STATUS='OLD',FORM='FORMATTED') READ(IN,10) TITLE 10 FORMAT(A) READ(IN,20) NPDE,MF,TOUT 20 FORMAT(2I3,E10.3) WRITE(IOUT,10) TITLE WRITE(IOUT,30) NPDE,MF,TOUT TKEEP = TOUT 30 FORMAT(//8H NPDE = ,I2,5X,5HMF = ,I2, * 5X,26H RESULTS ARE FOR TIME T = ,E10.3// *27H KORD NINT EPS ,24HNSTEPS NRES NJAC NQ , *1X,30H ERROR(1) ERROR(2) ) 40 READ(IN,50) KORD,NINT,EPS 50 FORMAT(2I3,E10.3) IF(KORD.LE.0) STOP NPTS = NINT + 1 NCC = 2 TOUT = TKEEP IWORK(1) = 30000 IWORK(2) = 1000 T0 = 0.0 DT = EPS * .5 INDEX = 1 NCPTS1 = KORD*NINT - NCC*(NINT-1) - 1 DX = 1. / FLOAT(NINT) DO 60 I=1,NPTS XBKPT(I) = FLOAT(I-1)*DX 60 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) C CHECK FOR ERRORS AND THEN COMPUTE AN ESTIMATE FOR THE OVERALL C MAXIMUM ERROR FOR EACH PDE COMPONENT. OUTPUT RESULTS. IF(INDEX.NE.0) WRITE(IOUT,70) INDEX 70 FORMAT(7H INDEX=,I3) DO 80 I=1,NPDE 80 ERRMX(I) = 0. C WORK(IW3) TO WORK(IW3+NCPTS1) CONTAINS THE COLLOCATION POINTS. DO 100 J=1,NCPTS1 DX = (WORK(IW3+J)-WORK(IW3+J-1))/24. DO 90 I=2,24 90 XPLT(I) = FLOAT(I-1)*DX + WORK(IW3+J-1) XPLT(1) = WORK(IW3+J-1) XPLT(25) = WORK(IW3+J) DO 100 I=1,25 CALL VALUES(XPLT(I),UPLT(1,1,I),SCRTCH,2,1,1,0,WORK) CALL TRUSOL(XPLT(I),TOUT,UEXACT) DO 100 K=1,NPDE UPLT(K,2,I) = ABS( UPLT(K,1,I) - UEXACT(K) ) IF(UPLT(K,2,I) .GT. ERRMX(K)) ERRMX(K) = UPLT(K,2,I) 100 CONTINUE WRITE(IOUT,110) KORD,NINT,EPS,NS,NF,NJ,NQUSED,(ERRMX(I),I=1,NPDE) 110 FORMAT(I4,I6,1E14.3,1X,4I6,2E15.3) GO TO 40 END C SUBROUTINES FOR TWO PDE PROBLEM. SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) DIMENSION U(NPDE),UX(NPDE),UXX(NPDE),FVAL(NPDE) TEMP = EXP(-4.*X) FVAL(1) = (U(2)-1.0)*UXX(1) + UX(1)*UX(2) + (16.*X*T-2.*T * - 16.*(U(2)-1.0))*(U(1)-1.0) + 10.*X*TEMP FVAL(2) = UXX(2) + UX(1) + 4.0*U(1) - 4.0 + X*X - 2.*T * - 10.*T*TEMP RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) DIMENSION U(NPDE),UX(NPDE),DZDT(NPDE),DBDU(NPDE,NPDE), * DBDUX(NPDE,NPDE) IF( X.NE. 0. ) GO TO 10 DZDT(1) = 0. DZDT(2) = 0. DBDU(1,1) = 1. DBDU(1,2) = 0. DBDU(2,1) = 0. DBDU(2,2) = 1. DBDUX(1,1) = 0. DBDUX(1,2) = 0. DBDUX(2,1) = 0. DBDUX(2,2) = 0. RETURN 10 CONTINUE DZDT(1) = 0. DZDT(2) = 0. DBDU(1,1) = 3. DBDU(1,2) = 0. DBDU(2,1) = -5.4598150033142E+01 DBDU(2,2) = 0. DBDUX(1,1) = 1. DBDUX(1,2) = 0. DBDUX(2,1) = 0. DBDUX(2,2) = 5. RETURN END SUBROUTINE UINIT(X,U,NPDE) DIMENSION U(NPDE) CALL TRUSOL(X,0.,U) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) DIMENSION U(NPDE),UX(NPDE),UXX(NPDE),DFDU(NPDE,NPDE) DIMENSION DFDUX(NPDE,NPDE),DFDUXX(NPDE,NPDE) DFDU(1,1) = 16.*X*T - 2.*T - 16.*(U(2)-1.0) DFDU(1,2) = UXX(1) - 16.*(U(1)-1.0) DFDU(2,1) = 4. DFDU(2,2) = 0. DFDUX(1,1) = UX(2) DFDUX(1,2) = UX(1) DFDUX(2,1) = 1. DFDUX(2,2) = 0. DFDUXX(1,1) = U(2) - 1.0 DFDUXX(1,2) = 0. DFDUXX(2,1) = 0. DFDUXX(2,2) = 1. RETURN END SUBROUTINE TRUSOL(X,T,UEXACT) DIMENSION UEXACT(2) UEXACT(1) = 10.*X*T*EXP(-4.*X) + 1.0 UEXACT(2) = T*X*X + 1.0 RETURN END C C*** Example Driver 4 C C PROGRAM BURG(TSIN,TSTOUT,TAPE2=TSIN,TAPE3=TSTOUT) C THIS PROGRAM WAS WRITTEN TO BE RUN IN SINGLE PRECISION ON A CDC7600. INTEGER IN, IOUT PARAMETER (IN=2, IOUT=6) DIMENSION XPLT(25),UPLT(2,2,404),SCRTCH(10) DIMENSION ERRMX(2),UEXACT(2),XBKPT(404) DIMENSION WORK(30000),IWORK(1000) CHARACTER *80 TITLE COMMON /GEAR0/ HUSED,NQUSED,NS,NF,NJ COMMON /ISTART/ IW1,IW2,IW3,IDUM(15) COMMON /PARAMS/ PI C READ IN AND INITIALIZE PDECOL CALLING ARGUMENTS AND PARAMETERS. PI = 4. * ATAN(1.0) OPEN(IN,FILE='EG4.DAT',STATUS='OLD',FORM='FORMATTED') READ(IN,10) TITLE 10 FORMAT(A) READ(IN,20) NPDE,MF,TOUT 20 FORMAT(2I3,E10.3) WRITE(IOUT,10) TITLE WRITE(IOUT,30) NPDE,MF,TOUT 30 FORMAT(//8H NPDE = ,I2,5X,5HMF = ,I2, * 5X,26H RESULTS ARE FOR TIME T = ,E10.3// *27H KORD NINT EPS ,24HNSTEPS NRES NJAC NQ , *1X,30H ERROR(1) ERROR(2) ) 40 READ(IN,50) KORD,NINT,EPS 50 FORMAT(2I3,E10.3) IF(KORD.LE.0) STOP NPTS = NINT + 1 NCC = 2 IWORK(1) = 30000 IWORK(2) = 1000 T0 = 0.0 DT = EPS * .5 INDEX = 1 DT = EPS * .5 INDEX = 1 NCPTS1 = KORD*NINT - NCC*(NINT-1) - 1 DX = 1. / FLOAT(NINT) DO 60 I=1,NPTS XBKPT(I) = FLOAT(I-1)*DX 60 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) C CHECK FOR ERRORS AND THEN COMPUTE AN ESTIMATE FOR THE OVERALL C MAXIMUM ERROR FOR EACH PDE COMPONENT. OUTPUT RESULTS. IF(INDEX.NE.0) WRITE(IOUT,70) INDEX 70 FORMAT(7H INDEX=,I3) DO 80 I=1,NPDE 80 ERRMX(I) = 0. C WORK(IW3) TO WORK(IW3+NCPTS1) CONTAINS THE COLLOCATION POINTS. DO 100 J=1,NCPTS1 DX = (WORK(IW3+J)-WORK(IW3+J-1))/24. DO 90 I=2,24 90 XPLT(I) = FLOAT(I-1)*DX + WORK(IW3+J-1) XPLT(1) = WORK(IW3+J-1) XPLT(25) = WORK(IW3+J) DO 100 I=1,25 CALL VALUES(XPLT(I),UPLT(1,1,I),SCRTCH,1,1,1,0,WORK) CALL TRUSOL(XPLT(I),TOUT,UEXACT) DO 100 K=1,NPDE UPLT(K,2,I) = ABS( UPLT(K,1,I) - UEXACT(K) ) IF(UPLT(K,2,I) .GT. ERRMX(K)) ERRMX(K) = UPLT(K,2,I) 100 CONTINUE WRITE(IOUT,110) KORD,NINT,EPS,NS,NF,NJ,NQUSED,(ERRMX(I),I=1,NPDE) 110 FORMAT(I4,I6,1E14.3,1X,4I6,2E15.3) GO TO 40 END C SUBROUTINES FOR BURGERS EQUATIONS PROBLEM. SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) FVAL(1) = .003*UXX(1) - U(1)*UX(1) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) IF( X.NE. 0. ) GO TO 10 DZDT(1) = 0. DBDU(1,1) = 1. DBDUX(1,1) = 0. RETURN 10 CONTINUE DZDT(1) = 0. DBDU(1,1) = 0. DBDUX(1,1) = 0. RETURN END SUBROUTINE UINIT(X,U,NPDE) DIMENSION U(NPDE) CALL TRUSOL(X,0.,U) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) DFDUXX(1,1) = .003 DFDUX(1,1) = -U(1) DFDU(1,1) = -UX(1) RETURN END SUBROUTINE TRUSOL(X,T,UEXACT) DIMENSION UEXACT(1) ANUIN = 1./.006 A = -(X-.5 + 4.95*T)*ANUIN*.1 B = -(X-.5 + 0.75*T)*ANUIN*.5 C = -(X-.375) * ANUIN A = EXP(A) B = EXP(B) C = EXP(C) UEXACT(1) = (.1*A + .5*B + C) / (A + B + C) RETURN END C C*** Example Driver 5 C C PROGRAM WAVE(TSIN,TSTOUT,TAPE2=TSIN,TAPE3=TSTOUT) C THIS PROGRAM WAS WRITTEN TO BE RUN IN SINGLE PRECISION ON A CDC7600. INTEGER IN, IOUT PARAMETER (IN=2, IOUT=6) DIMENSION XPLT(25),UPLT(2,2,404),SCRTCH(10) DIMENSION ERRMX(2),UEXACT(2),XBKPT(404) DIMENSION WORK(30000),IWORK(1000) CHARACTER *80 TITLE COMMON /GEAR0/ HUSED,NQUSED,NS,NF,NJ COMMON /ISTART/ IW1,IW2,IW3,IDUM(15) COMMON /PARAMS/ PI C READ IN AND INITIALIZE PDECOL CALLING ARGUMENTS AND PARAMETERS. PI = 4. * ATAN(1.0) OPEN(IN,FILE='EG5.DAT',STATUS='OLD',FORM='FORMATTED') READ(IN,10) TITLE 10 FORMAT(A) READ(IN,20) NPDE,MF,TOUT 20 FORMAT(2I3,E10.3) WRITE(IOUT,10) TITLE WRITE(IOUT,30) NPDE,MF,TOUT 30 FORMAT(//8H NPDE = ,I2,5X,5HMF = ,I2, * 5X,26H RESULTS ARE FOR TIME T = ,E10.3// *27H KORD NINT EPS ,24HNSTEPS NRES NJAC NQ , *1X,30H ERROR(1) ERROR(2) ) TKEEP = TOUT 40 READ(IN,50) KORD,NINT,EPS 50 FORMAT(2I3,E10.3) IF(KORD.LE.0) STOP NPTS = NINT + 1 NCC = 2 TOUT=TKEEP IWORK(1) = 30000 IWORK(2) = 1000 T0 = 0.0 DT = EPS * .5 INDEX = 1 NCPTS1 = KORD*NINT - NCC*(NINT-1) - 1 DX = 1. / FLOAT(NINT) DO 60 I=1,NPTS XBKPT(I) = FLOAT(I-1)*DX 60 CONTINUE CALL PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF,INDEX, * WORK,IWORK) C CHECK FOR ERRORS AND THEN COMPUTE AN ESTIMATE FOR THE OVERALL C MAXIMUM ERROR FOR EACH PDE COMPONENT. OUTPUT RESULTS. IF(INDEX.NE.0) WRITE(IOUT,70) INDEX 70 FORMAT(7H INDEX=,I3) DO 80 I=1,NPDE 80 ERRMX(I) = 0. C WORK(IW3) TO WORK(IW3+NCPTS1) CONTAINS THE COLLOCATION POINTS. DO 100 J=1,NCPTS1 DX = (WORK(IW3+J)-WORK(IW3+J-1))/24. DO 90 I=2,24 90 XPLT(I) = FLOAT(I-1)*DX + WORK(IW3+J-1) XPLT(1) = WORK(IW3+J-1) XPLT(25) = WORK(IW3+J) DO 100 I=1,25 CALL VALUES(XPLT(I),UPLT(1,1,I),SCRTCH,1,1,1,0,WORK) CALL TRUSOL(XPLT(I),TOUT,UEXACT) DO 100 K=1,NPDE UPLT(K,2,I) = ABS( UPLT(K,1,I) - UEXACT(K) ) IF(UPLT(K,2,I) .GT. ERRMX(K)) ERRMX(K) = UPLT(K,2,I) 100 CONTINUE WRITE(IOUT,110) KORD,NINT,EPS,NS,NF,NJ,NQUSED,(ERRMX(I),I=1,NPDE) 110 FORMAT(I4,I6,1E14.3,1X,4I6,2E15.3) GO TO 40 END C SUBROUTINES FOR WAVE PROPAGATION PROBLEM. SUBROUTINE F(T,X,U,UX,UXX,FVAL,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) FVAL(1) = -UX(1) RETURN END SUBROUTINE BNDRY(T,X,U,UX,DBDU,DBDUX,DZDT,NPDE) DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) COMMON /PARAMS/ PI IF( X.NE. 0. ) GO TO 10 DZDT(1) = -2.*PI*COS(-4.*PI*T) DBDU(1,1) = 1. DBDUX(1,1) = 0. RETURN 10 CONTINUE DZDT(1) = 0. DBDU(1,1) = 0. DBDUX(1,1) = 0. RETURN END SUBROUTINE UINIT(X,U,NPDE) DIMENSION U(NPDE) CALL TRUSOL(X,0.,U) RETURN END SUBROUTINE DERIVF(T,X,U,UX,UXX,DFDU,DFDUX,DFDUXX,NPDE) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) DFDU(1,1) = 0. DFDUX(1,1) = -1. DFDUXX(1,1) = 0. RETURN END SUBROUTINE TRUSOL(X,T,UEXACT) DIMENSION UEXACT(1) COMMON /PARAMS/ PI UEXACT(1) = 1.0 + .5*SIN(4.*PI*(X-T)) RETURN END C*** Example 1 Data SINE PROBLEM INPUT DATA 1 21 .1 3 1 .01 3 3 .001 3 10 .00003 3 20 .00001 3 60 .000001 3100 .000001 4 1 .01 4 2 .001 4 4 .00001 4 8 .000001 4 16 .0000001 5 1 .0001 5 2 .00001 5 4 .000001 5 8 .00000001 6 1 .0001 6 2 .000001 6 4 .00000001 -1 C*** Example 2 Data SINE PROBLEM INPUT DATA 1 21 .1 3 1 .01 3 3 .001 3 10 .00003 3 20 .00001 3 60 .000001 3100 .000001 4 1 .01 4 2 .001 4 4 .00001 4 8 .000001 4 16 .0000001 5 1 .0001 5 2 .00001 5 4 .000001 5 8 .00000001 6 1 .0001 6 2 .000001 6 4 .00000001 -1 C*** Example 3 Data TWO PDE PROBLEM INPUT DATA 2 21 1.0 3 3 .0027 3 6 .00064 3 12 .00016 3 24 .00004 3 48 .00001 3 96 .0000025 4 1 .015 4 2 .0013 4 4 .00007 4 8 .000004 4 16 .00000025 4 32.000000015 4 64.000000001 5 1 .0035 5 2 .0001 5 4 .000003 5 8 .00000006 5 16.000000001 6 1 .00046 6 2 .000013 6 4 .00000013 6 8.000000001 -1 C*** Example 4 Data BURGERS EQUATION PROBLEM INPUT DATA 1 21 .8 3 50 .003 3100 .0003 3200 .0001 3400 .00004 4 20 .004 4 40 .001 4 80 .0001 4160 .00001 5 15 .004 5 30 .001 5 60 .00003 5100 .000003 6 10 .01 6 20 .001 6 40 .00003 6 80 .000001 -1 C*** Example 5 Data WAVE PROPAGATION PROBLEM INPUT DATA 1 11 .5 3 20 .001 3 40 .001 3 80 .00001 3160 .000001 3320 .0000001 4 5 .001 4 10 .0001 4 20 .00001 4 40 .000001 4 80 .00000001 41601.0000E-10 5 3 .001 5 6 .0004 5 12 .000002 5 24 .00000001 5 481.0000E-10 6 2 .001 6 4 .0001 6 8 .0000001 6 16 .00000001 6 321.0000E-10 -1 *** Example 1 Results T= 0.100E-02 DT= 0.340E-03 TOTAL STEPS= 11 PDE COMPONENT 1 0.5000E+00 0.5191E+00 0.5376E+00 0.5555E+00 0.5729E+00 0.5901E+00 0.6071E+00 0.6239E+00 0.6407E+00 0.6574E+00 0.6740E+00 0.6906E+00 0.7073E+00 0.7239E+00 0.7405E+00 0.7571E+00 0.7736E+00 0.7902E+00 0.8068E+00 0.8234E+00 0.8400E+00 0.8566E+00 0.8732E+00 0.8898E+00 0.9064E+00 0.9231E+00 0.9398E+00 0.9565E+00 0.9733E+00 0.9903E+00 0.1007E+01 PDE COMPONENT 2 0.3142E+01 0.3133E+01 0.3133E+01 0.3133E+01 0.3133E+01 0.3133E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3134E+01 0.3135E+01 0.3135E+01 0.3135E+01 0.3135E+01 0.3135E+01 0.3135E+01 0.3135E+01 0.3135E+01 T= 0.100E-01 DT= 0.205E-02 TOTAL STEPS= 21 PDE COMPONENT 1 0.5000E+00 0.5248E+00 0.5492E+00 0.5730E+00 0.5961E+00 0.6185E+00 0.6402E+00 0.6613E+00 0.6819E+00 0.7020E+00 0.7216E+00 0.7410E+00 0.7600E+00 0.7788E+00 0.7974E+00 0.8159E+00 0.8343E+00 0.8528E+00 0.8713E+00 0.8901E+00 0.9090E+00 0.9283E+00 0.9480E+00 0.9682E+00 0.9890E+00 0.1011E+01 0.1033E+01 0.1056E+01 0.1081E+01 0.1107E+01 0.1134E+01 PDE COMPONENT 2 0.3142E+01 0.3094E+01 0.3072E+01 0.3063E+01 0.3061E+01 0.3060E+01 0.3061E+01 0.3062E+01 0.3063E+01 0.3064E+01 0.3065E+01 0.3066E+01 0.3067E+01 0.3068E+01 0.3068E+01 0.3069E+01 0.3070E+01 0.3071E+01 0.3072E+01 0.3072E+01 0.3073E+01 0.3074E+01 0.3075E+01 0.3076E+01 0.3077E+01 0.3079E+01 0.3080E+01 0.3081E+01 0.3083E+01 0.3084E+01 0.3086E+01 T= 0.100E+00 DT= 0.567E-02 TOTAL STEPS= 51 PDE COMPONENT 1 0.5000E+00 0.5472E+00 0.5970E+00 0.6487E+00 0.7015E+00 0.7546E+00 0.8075E+00 0.8598E+00 0.9112E+00 0.9615E+00 0.1010E+01 0.1058E+01 0.1104E+01 0.1148E+01 0.1191E+01 0.1232E+01 0.1271E+01 0.1309E+01 0.1346E+01 0.1380E+01 0.1414E+01 0.1446E+01 0.1477E+01 0.1506E+01 0.1534E+01 0.1561E+01 0.1587E+01 0.1612E+01 0.1635E+01 0.1658E+01 0.1680E+01 PDE COMPONENT 2 0.3142E+01 0.3006E+01 0.2909E+01 0.2841E+01 0.2794E+01 0.2765E+01 0.2748E+01 0.2741E+01 0.2743E+01 0.2752E+01 0.2766E+01 0.2786E+01 0.2809E+01 0.2836E+01 0.2866E+01 0.2899E+01 0.2935E+01 0.2973E+01 0.3014E+01 0.3057E+01 0.3102E+01 0.3149E+01 0.3198E+01 0.3249E+01 0.3303E+01 0.3358E+01 0.3415E+01 0.3475E+01 0.3536E+01 0.3599E+01 0.3665E+01 T= 0.100E+01 DT= 0.381E+00 TOTAL STEPS= 76 PDE COMPONENT 1 0.5000E+00 0.5471E+00 0.5969E+00 0.6487E+00 0.7014E+00 0.7545E+00 0.8073E+00 0.8593E+00 0.9102E+00 0.9597E+00 0.1008E+01 0.1054E+01 0.1099E+01 0.1142E+01 0.1183E+01 0.1222E+01 0.1260E+01 0.1296E+01 0.1331E+01 0.1364E+01 0.1396E+01 0.1426E+01 0.1455E+01 0.1483E+01 0.1509E+01 0.1535E+01 0.1559E+01 0.1583E+01 0.1605E+01 0.1627E+01 0.1647E+01 PDE COMPONENT 2 0.3142E+01 0.3003E+01 0.2904E+01 0.2835E+01 0.2790E+01 0.2762E+01 0.2748E+01 0.2746E+01 0.2753E+01 0.2767E+01 0.2788E+01 0.2813E+01 0.2842E+01 0.2875E+01 0.2910E+01 0.2949E+01 0.2990E+01 0.3033E+01 0.3078E+01 0.3125E+01 0.3174E+01 0.3224E+01 0.3276E+01 0.3330E+01 0.3386E+01 0.3443E+01 0.3501E+01 0.3562E+01 0.3624E+01 0.3688E+01 0.3754E+01 T= 0.100E+02 DT= 0.415E+01 TOTAL STEPS= 83 PDE COMPONENT 1 0.5000E+00 0.5471E+00 0.5969E+00 0.6487E+00 0.7014E+00 0.7545E+00 0.8073E+00 0.8593E+00 0.9102E+00 0.9597E+00 0.1008E+01 0.1054E+01 0.1099E+01 0.1142E+01 0.1183E+01 0.1222E+01 0.1260E+01 0.1296E+01 0.1331E+01 0.1364E+01 0.1396E+01 0.1426E+01 0.1455E+01 0.1483E+01 0.1509E+01 0.1535E+01 0.1559E+01 0.1583E+01 0.1605E+01 0.1627E+01 0.1647E+01 PDE COMPONENT 2 0.3142E+01 0.3003E+01 0.2904E+01 0.2835E+01 0.2790E+01 0.2762E+01 0.2748E+01 0.2746E+01 0.2753E+01 0.2767E+01 0.2788E+01 0.2813E+01 0.2842E+01 0.2875E+01 0.2910E+01 0.2949E+01 0.2990E+01 0.3033E+01 0.3078E+01 0.3125E+01 0.3174E+01 0.3224E+01 0.3276E+01 0.3330E+01 0.3386E+01 0.3443E+01 0.3501E+01 0.3562E+01 0.3624E+01 0.3688E+01 0.3754E+01 INDEX= 0 *** Example 2 Results SINE PROBLEM INPUT DATA 28670 NPDE = 1 MF = 21 RESULTS ARE FOR TIME T = 0.100E+00 KORD NINT EPS NSTEPS NRES NJAC NQ ERROR(1) ERROR(2) 3 1 0.100E-01 7 9 3 2 0.687E-01 3 3 0.100E-02 9 14 4 3 0.163E-01 3 10 0.300E-04 19 26 6 4 0.116E-02 3 20 0.100E-04 22 39 7 4 0.300E-03 3 60 0.100E-05 47 73 10 3 0.339E-04 3 100 0.100E-05 56 83 10 3 0.870E-05 4 1 0.100E-01 7 9 3 2 0.916E-01 4 2 0.100E-02 10 16 4 3 0.492E-02 4 4 0.100E-04 22 38 7 4 0.395E-03 4 8 0.100E-05 36 59 9 4 0.274E-04 4 16 0.100E-06 161 286 11 2 0.691E-05 5 1 0.100E-03 14 20 4 4 0.273E-02 5 2 0.100E-04 22 38 7 4 0.348E-03 5 4 0.100E-05 33 46 8 4 0.132E-04 INTEGRATION HALTED BY DRIVER AT T = 0.50000000E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 5 8 0.100E-07 1 1 2 1 0.238E-06 6 1 0.100E-03 16 23 5 3 0.177E-02 6 2 0.100E-05 30 45 8 4 0.201E-04 INTEGRATION HALTED BY DRIVER AT T = 0.50000000E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 6 4 0.100E-07 1 1 2 1 0.238E-06 *** Example 3 Results TWO PDE PROBLEM INPUT DATA 30140 NPDE = 2 MF = 21 RESULTS ARE FOR TIME T = 0.100E+01 KORD NINT EPS NSTEPS NRES NJAC NQ ERROR(1) ERROR(2) 3 3 0.270E-02 10 14 5 1 0.133E+00 0.261E-01 3 6 0.640E-03 9 13 4 1 0.130E-01 0.637E-02 3 12 0.160E-03 11 22 5 1 0.125E-01 0.159E-02 3 24 0.400E-04 15 25 6 2 0.493E-02 0.404E-03 3 48 0.100E-04 19 34 8 2 0.149E-02 0.103E-03 3 96 0.250E-05 218 417 117 1 0.655E-03 0.523E-04 4 1 0.150E-01 10 17 5 1 0.262E+00 0.163E+00 4 2 0.130E-02 10 13 5 1 0.528E-01 0.131E-01 4 4 0.700E-04 11 16 5 1 0.398E-02 0.677E-03 4 8 0.400E-05 18 24 7 1 0.502E-03 0.409E-04 4 16 0.250E-06 400 812 204 1 0.401E-04 0.894E-05 INTEGRATION HALTED BY DRIVER AT T = 0.74999997E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 4 32 0.150E-07 1 1 2 1 0.238E-06 0.238E-06 INTEGRATION HALTED BY DRIVER AT T = 0.49999999E-09 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 4 64 0.100E-08 1 1 2 1 0.238E-06 0.238E-06 5 1 0.350E-02 11 18 5 1 0.998E-01 0.358E-01 5 2 0.100E-03 13 23 6 2 0.715E-02 0.116E-02 5 4 0.300E-05 19 32 7 2 0.362E-03 0.308E-04 5 8 0.600E-07 1795 4105 841 1 0.196E-04 0.296E-04 INTEGRATION HALTED BY DRIVER AT T = 0.49999999E-09 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 5 16 0.100E-08 1 1 2 1 0.238E-06 0.238E-06 6 1 0.460E-03 13 25 5 2 0.190E-01 0.524E-02 6 2 0.130E-04 15 21 6 2 0.681E-03 0.118E-03 6 4 0.130E-06 218 506 109 1 0.225E-04 0.215E-05 INTEGRATION HALTED BY DRIVER AT T = 0.49999999E-09 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 6 8 0.100E-08 1 1 2 1 0.238E-06 0.238E-06 *** Example 4 Results BURGERS EQUATION PROBLEM INPUT DATA 31410 NPDE = 1 MF = 21 RESULTS ARE FOR TIME T = 0.800E+00 KORD NINT EPS NSTEPS NRES NJAC NQ ERROR(1) ERROR(2) 3 50 0.300E-02 60 134 6 2 0.130E+00 3 100 0.300E-03 117 204 5 3 0.226E-01 3 200 0.100E-03 148 257 7 3 0.558E-02 3 400 0.400E-04 163 291 8 4 0.804E-03 4 20 0.400E-02 60 129 4 3 0.990E-01 4 40 0.100E-02 94 200 5 3 0.270E-01 4 80 0.100E-03 159 292 7 3 0.326E-02 4 160 0.100E-04 210 352 9 4 0.277E-03 5 15 0.400E-02 72 158 6 2 0.935E-01 5 30 0.100E-02 110 225 6 2 0.206E-01 5 60 0.300E-04 182 334 9 4 0.574E-03 5 100 0.300E-05 265 409 12 5 0.832E-03 6 10 0.100E-01 55 143 7 3 0.122E+00 6 20 0.100E-02 116 238 7 3 0.145E-01 6 40 0.300E-04 191 364 9 4 0.381E-03 KFLAG = -1 FROM INTEGRATOR AT T = 0.74923939E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000001E-07 AND STEP WILL BE RETRIED KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000000E-08 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.49999999E-09 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000001E-10 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000000E-11 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000000E-12 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.49999999E-13 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.49999999E-14 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000000E-15 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN DT HAS BEEN REDUCED TO 0.50000001E-16 AND STEP WILL BE RETRIED WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. WARNING.. T + DT = T ON NEXT STEP. KFLAG = -1 FROM INTEGRATOR AT T = 0.74924850E+00 ERROR TEST FAILED WITH ABS(DT) = DTMIN PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT INDEX= -1 6 80 0.100E-05 1009 1280 101 5 0.763E+12 *** Example 5 Results WAVE PROPAGATION PROBLEM INPUT DATA 32550 NPDE = 1 MF = 11 RESULTS ARE FOR TIME T = 0.500E+00 KORD NINT EPS NSTEPS NRES NJAC NQ ERROR(1) ERROR(2) 3 20 0.100E-02 19 35 5 4 0.587E-01 3 40 0.100E-02 19 35 5 4 0.139E-01 3 80 0.100E-04 49 92 7 4 0.346E-02 3 160 0.100E-05 117 176 9 4 0.947E-03 3 320 0.100E-06 319 375 14 3 0.275E-03 4 5 0.100E-02 20 37 5 4 0.709E-01 4 10 0.100E-03 28 52 5 4 0.529E-02 4 20 0.100E-04 36 60 6 5 0.485E-03 4 40 0.100E-05 56 101 7 5 0.586E-04 INTEGRATION HALTED BY DRIVER AT T = 0.50000000E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 4 80 0.100E-07 1 1 2 1 0.715E-06 INTEGRATION HALTED BY DRIVER AT T = 0.50000001E-10 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 4 160 0.100E-09 1 1 2 1 0.477E-06 5 3 0.100E-02 20 37 5 4 0.552E-01 5 6 0.400E-03 23 43 5 4 0.105E-01 5 12 0.200E-05 56 92 7 4 0.122E-03 INTEGRATION HALTED BY DRIVER AT T = 0.50000000E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 5 24 0.100E-07 1 1 2 1 0.107E-05 INTEGRATION HALTED BY DRIVER AT T = 0.50000001E-10 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 5 48 0.100E-09 1 1 2 1 0.596E-06 6 2 0.100E-02 21 39 5 4 0.185E+00 6 4 0.100E-03 30 56 5 4 0.227E-02 6 8 0.100E-06 82 122 8 5 0.386E-04 INTEGRATION HALTED BY DRIVER AT T = 0.50000000E-08 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 6 16 0.100E-07 1 1 2 1 0.596E-06 INTEGRATION HALTED BY DRIVER AT T = 0.50000001E-10 EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION INDEX= -2 6 32 0.100E-09 1 1 2 1 0.536E-06