----------------------------------------------------------------------- P D E C O L ----------------------------------------------------------------------- This is the March 24, 1978 version of PDECOL. PDECOL is a sophisticated package of subroutines which is designed to solve the general system of NPDE nonlinear partial differential equations of at most second order on the interval (XLEFT,XRIGHT) for t .GT. T0 which is of the form.... du/dt = f( t, x, u, ux, uxx ) where u = ( u(1), u(2), ... , u(NPDE) ) ux = ( ux(1), ux(2), ... , ux(NPDE) ) uxx = (uxx(1),uxx(2), ... ,uxx(NPDE) ) . Each u(k) is a function of the scalar quantities t and x. Ux(k) represents the first partial derivative of u(k) with respect to the variable x, uxx(k) represents the second partial derivative of u(k) with respect to the variable x, and du/dt is the vector of partial derivatives of u with respect to the time variable t. F represents an arbitrary vector valued function whose NPDE components define the respective partial differential equations of the PDE system. Boundary conditions Depending on the type of PDE(s), 0, 1, or 2 boundary conditions are required for each PDE in the system. These are imposed at XLEFT and/or XRIGHT and each must be of the form.... b(u,ux) = z(t) where b and z are arbitrary vector valued functions with NPDE components and u, ux, and t are as above. These boundary conditions must be consistent with the initial conditions which are described next. Initial conditions Each solution component u(k) is assumed to be a known (user provided) function of x at the initial time t = t0. The initial condition functions must be consistent with the boundary conditions above, i.E. The initial condition functions must satisfy the boundary conditions for t = t0. ----------------------------------------------------------------------- References 1. Madsen, N.K. and R.F. Sincovec, PDECOL - collocation software for partial differential equations, ACM-TOMS, vol. 5, No. 3, Sept. 1979, pp. 326-351. 2. Sincovec, R.F. and N.K. Madsen, Software for nonlinear partial differential equations, ACM-TOMS, vol. 1, No. 3, September 1975, pp. 232-260. 3. Hindmarsh, A.C., Preliminary documentation of GEARIB.. Solution of implicit systems of ordinary differential equations with banded Jacobians, Lawrence Livermore Lab, ucid-30130, February 1976. 4. DeBoor, C., Package for calculating with b-splines, SIAM J. Numer. Anal., Vol. 14, No. 3, June 1977, pp. 441-472. ---------------------------------------------------------------------- DETAILS ON USE OF PDECOL ---------------------------------------------------------------------- SUBROUTINE PDECOL(T0,TOUT,DT,XBKPT,EPS,NINT,KORD,NCC,NPDE,MF, * INDEX,WORK,IWORK) ----------------------------------------------------------------------- REQUIRED USER SUPPLIED SUBROUTINES THE USER IS REQUIRED TO CONSTRUCT THREE SUBPROGRAMS AND A MAIN PROGRAM WHICH DEFINE THE PDE PROBLEM WHOSE SOLUTION IS TO BE ATTEMPTED. THE THREE SUBPROGRAMS ARE... 1) SUBROUTINE F( T, X, U, UX, UXX, FVAL, NPDE ) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) THIS ROUTINE DEFINES THE DESIRED PARTIAL DIFFERENTIAL EQUATIONS TO BE SOLVED. THE PACKAGE PROVIDES VALUES OF THE INPUT SCALARS T AND X AND INPUT ARRAYS (LENGTH NPDE) U, UX, AND UXX, AND THE USER MUST CONSTRUCT THIS ROUTINE TO COMPUTE THE OUTPUT ARRAY FVAL (LENGTH NPDE) WHICH CONTAINS THE CORRESPONDING VALUES OF THE RIGHT HAND SIDES OF THE DESIRED PARTIAL DIFFERENTIAL EQUATIONS, I.E. FVAL(K) = THE VALUE OF THE RIGHT HAND SIDE OF THE K-TH PDE IN THE PDE SYSTEM ABOVE, FOR K = 1 TO NPDE. THE INCOMING VALUE OF THE SCALAR QUANTITY X WILL BE A COLLOCATION POINT VALUE (SEE INITAL AND COLPNT) AND THE INCOMING VALUES IN THE ARRAYS U, UX AND UXX CORRESPOND TO THIS POINT X AND TIME T. RETURN END 2) SUBROUTINE BNDRY( T, X, U, UX, DBDU, DBDUX, DZDT, NPDE ) DIMENSION U(NPDE), UX(NPDE), DZDT(NPDE) DIMENSION DBDU(NPDE,NPDE), DBDUX(NPDE,NPDE) THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH NEEDED INFORMATION ABOUT THE BOUNDARY CONDITION FUNCTIONS B AND Z ABOVE. THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, AND UX, AND THE USER IS TO DEFINE THE CORRESPONDING OUTPUT VALUES OF THE DERIVATIVES OF THE FUNCTIONS B AND Z WHERE.... DBDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO THE J-TH VARIABLE U(J). DBDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR FUNCTION B(U,UX) ABOVE WITH RESPECT TO THE J-TH VARIABLE UX(J). DZDT(K) = DERIVATIVE OF THE K-TH COMPONENT OF THE VECTOR FUNCTION Z(T) ABOVE WITH RESPECT TO THE VARIABLE T. NOTE... THE INCOMING VALUE OF X WILL BE EITHER XLEFT OR XRIGHT. IF NO BOUNDARY CONDITION IS DESIRED FOR SAY THE K-TH PDE AT ONE OR BOTH OF THE ENDPOINTS XLEFT OR XRIGHT, THEN DBDU(K,K) AND DBDUX(K,K) SHOULD BOTH BE SET TO ZERO WHEN BNDRY IS CALLED FOR THAT POINT. WE REFER TO THIS AS A NULL BOUNDARY CONDITION. THIS ROUTINE CAN BE STRUCTURED AS FOLLOWS... THE COMMON BLOCK /ENDPT/ IS NOT A PART OF PDECOL AND MUST BE SUPPLIED AND DEFINED BY THE USER. COMMON /ENDPT/ XLEFT IF( X .NE. XLEFT ) GO TO 10 HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE LEFT BOUNDARY POINT X = XLEFT. RETURN 10 CONTINUE HERE DEFINE AND SET PROPER VALUES FOR DBDU(K,J), DBDUX(K,J), AND DZDT(K) FOR K,J = 1 TO NPDE FOR THE RIGHT BOUNDARY POINT X = XRIGHT. RETURN END 3) SUBROUTINE UINIT( X, U, NPDE ) DIMENSION U(NPDE) THIS ROUTINE IS USED TO PROVIDE THE PDE PACKAGE WITH THE NEEDED INITIAL CONDITION FUNCTION VALUES. THE PACKAGE PROVIDES A VALUE OF THE INPUT VARIABLE X, AND THE USER IS TO DEFINE THE PROPER INITIAL VALUES (AT T = T0) FOR ALL OF THE PDE COMPONENTS, I.E. U(K) = DESIRED INITIAL VALUE OF PDE COMPONENT U(K) AT X AND T = T0 FOR K = 1 TO NPDE. NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT VALUE. THE INITIAL CONDITIONS AND BOUNDARY CONDITIONS MUST BE CONSISTENT (SEE ABOVE). RETURN END ----------------------------------------------------------------------- OPTIONAL USER SUPPLIED SUBROUTINE IF THE USER DESIRES TO USE THE MF = 11 OR 21 OPTION IN ORDER TO SAVE ABOUT 10-20 PERCENT IN EXECUTION TIME (SEE BELOW), THEN THE USER MUST PROVIDE THE FOLLOWING SUBROUTINE WHICH PROVIDES INFORMATION ABOUT THE DERIVATIVES OF THE FUNCTION F ABOVE. THIS PROVIDES FOR MORE EFFICIENT JACOBIAN MATRIX GENERATION. ON MOST COMPUTER SYSTEMS, THE USER WILL BE REQUIRED TO SUPPLY THIS SUBROUTINE AS A DUMMY SUBROUTINE IF THE OPTIONS MF = 12 OR 22 ARE USED (SEE BELOW). 1) 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) THE PACKAGE PROVIDES VALUES OF THE INPUT VARIABLES T, X, U, UX, AND UXX, AND THE USER SHOULD CONSTRUCT THIS ROUTINE TO PROVIDE THE FOLLOWING CORRESPONDING VALUES OF THE OUTPUT ARRAYS DFDU, DFDUX, AND DFDUXX FOR K,J = 1 TO NPDE... DFDU(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE PDE DEFINING FUNCTION F WITH RESPECT TO THE VARIABLE U(J). DFDUX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE PDE DEFINING FUNCTION F WITH RESPECT TO THE VARIABLE UX(J). DFDUXX(K,J) = PARTIAL DERIVATIVE OF THE K-TH COMPONENT OF THE PDE DEFINING FUNCTION F WITH RESPECT TO THE VARIABLE UXX(J). NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT VALUE. RETURN END ----------------------------------------------------------------------- METHODS USED THE PACKAGE PDECOL IS BASED ON THE METHOD OF LINES AND USES A FINITE ELEMENT COLLOCATION PROCEDURE (WITH PIECEWISE POLYNOMIALS AS THE TRIAL SPACE) FOR THE DISCRETIZATION OF THE SPATIAL VARIABLE X. THE COLLOCATION PROCEDURE REDUCES THE PDE SYSTEM TO A SEMI- DISCRETE SYSTEM WHICH THEN DEPENDS ONLY ON THE TIME VARIABLE T. THE TIME INTEGRATION IS THEN ACCOMPLISHED BY USE OF SLIGHTLY MODIFIED STANDARD TECHNIQUES (SEE REFS. 1,2). PIECEWISE POLYNOMIALS THE USER IS REQUIRED TO SELECT THE PIECEWISE POLYNOMIAL SPACE WHICH IS TO BE USED TO COMPUTE HIS APPROXIMATE SOLUTION. FIRST, THE ORDER, KORD, OF THE POLYNOMIALS TO BE USED MUST BE SPECIFIED (KORD = POLYNOMIAL DEGREE + 1). NEXT, THE NUMBER OF PIECES (INTERVALS), NINT, INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS TO BE DIVIDED, IS CHOSEN. THE NINT + 1 DISTINCT BREAKPOINTS OF THE DOMAIN MUST BE DEFINED AND SET INTO THE ARRAY XBKPT IN STRICTLY INCREASING ORDER, I.E. XLEFT=XBKPT(1) .LT. XBKPT(2) .LT. ... .LT. XBKPT(NINT+1)=XRIGHT. THE APPROXIMATE SOLUTION AT ANY TIME T WILL BE A POLYNOMIAL OF ORDER KORD OVER EACH SUBINTERVAL (XBKPT(I),XBKPT(I+1)). THE NUMBER OF CONTINUITY CONDITIONS, NCC, TO BE IMPOSED ACROSS ALL OF THE BREAKPOINTS IS THE LAST PIECE OF USER SUPPLIED DATA WHICH IS REQUIRED TO UNIQUELY DETERMINE THE DESIRED PIECEWISE POLYNOMIAL SPACE. FOR EXAMPLE, NCC = 2 WOULD REQUIRE THAT THE APPROXIMATE SOLUTION (MADE UP OF THE SEPARATE POLYNOMIAL PIECES) AND ITS FIRST SPATIAL DERIVATIVE BE CONTINUOUS AT THE BREAKPOINTS AND HENCE ON THE ENTIRE DOMAIN (XLEFT,XRIGHT). NCC = 3 WOULD REQUIRE THAT THE APPROXIMATE SOLUTION AND ITS FIRST AND SECOND SPATIAL DERIVATIVES BE CONTINUOUS AT THE BREAKPOINTS, ETC. THE DIMENSION OF THIS LINEAR SPACE IS KNOWN AND FINITE AND IS NCPTS = KORD*NINT - NCC*(NINT-1). THE WELL-KNOWN B-SPLINE BASIS (SEE REF. 3) FOR THIS SPACE IS USED BY PDECOL AND IT CONSISTS OF NCPTS KNOWN PIECEWISE POLYNOMIAL FUNCTIONS BF(I,X), FOR I=1 TO NCPTS, WHICH DO NOT DEPEND ON THE TIME VARIABLE T. WE WISH TO EMPHASIZE THAT THE PIECEWISE POLYNOMIAL SPACE USED IN PDECOL (WHICH IS SELECTED BY THE USER) WILL DETERMINE THE MAGNITUDE OF THE SPATIAL DISCRETIZATION ERRORS IN THE COMPUTED APPROXIMATE SOLUTION. THE PACKAGE HAS NO CONTROL OVER ERRORS INTRODUCED BY THE USERS CHOICE OF THIS SPACE. SEE INPUT PARAMETERS BELOW. COLLOCATION OVER PIECEWISE POLYNOMIALS THE BASIC ASSUMPTION MADE IS THAT THE APPROXIMATE SOLUTION SATISFIES NCPTS U(T,X) = SUM C(I,T) * BF(I,X) I=1 WHERE THE UNKNOWN COEFFICIENTS C DEPEND ONLY ON THE TIME T AND THE KNOWN BASIS FUNCTIONS DEPEND ONLY ON X (WE HAVE ASSUMED THAT NPDE = 1 FOR CONVENIENCE). SO, AT ANY GIVEN TIME T THE APPROX- IMATE SOLUTION IS A PIECEWISE POLYNOMIAL IN THE USER CHOSEN SPACE. THE SEMI-DISCRETE EQUATIONS (ACTUALLY ORDINARY DIFFERENTIAL EQUATIONS) WHICH DETERMINE THE COEFFICIENTS C ARE OBTAINED BY REQUIRING THAT THE ABOVE APPROXIMATE U(T,X) SATISFY THE PDE AND BOUNDARY CONDITIONS EXACTLY AT A SET OF NCPTS COLLOCATION POINTS (SEE COLPNT). THUS, PDECOL ACTUALLY COMPUTES THE BASIS FUNCTION COEFFICIENTS RATHER THAN SPECIFIC APPROXIMATE SOLUTION VALUES. ----------------------------------------------------------------------- USE OF PDECOL PDECOL IS CALLED ONCE FOR EACH DESIRED OUTPUT VALUE (TOUT) OF THE TIME T, AND IT IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, STIFIB, WHICH ADVANCES THE TIME BY TAKING SINGLE STEPS UNTIL T .GE. TOUT. INTERPOLATION TO THE EXACT TIME TOUT IS THEN DONE. SEE TOUT BELOW. SUMMARY OF SUGGESTED INPUT VALUES IT IS OF COURSE IMPOSSIBLE TO SUGGEST INPUT PARAMETER VALUES WHICH ARE APPROPRIATE FOR ALL PROBLEMS. THE FOLLOWING SUGGESTIONS ARE TO BE USED ONLY IF YOU HAVE NO IDEA OF BETTER VALUES FOR YOUR PROBLEM. DT = 1.E-10 XBKPT = CHOOSE NINT+1 EQUALLY SPACED VALUES SUCH THAT XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT. EPS = 1.E-4 NINT = ENOUGH SO THAT ANY FINE STRUCTURE OF THE PROBLEM MAY BE RESOLVED. KORD = 4 NCC = 2 MF = 22 INDEX = 1 (ON FIRST CALL ONLY, THEN 0 THEREAFTER). THE INPUT PARAMETERS ARE.. T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE (USED ONLY ON FIRST CALL). TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT. SINCE THE PACKAGE CHOOSES ITS OWN TIME STEP SIZES, THE INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT AND THE PACKAGE WILL INTERPOLATE TO T = TOUT. DT = THE INITIAL STEP SIZE IN T, IF INDEX = 1, OR, THE MAXIMUM STEP SIZE ALLOWED (MUST BE .GT. 0), IF INDEX = 3. USED FOR INPUT ONLY WHEN INDEX = 1 OR 3. SEE BELOW. XBKPT = THE ARRAY OF PIECEWISE POLYNOMIAL BREAKPOINTS. THE NINT+1 VALUES MUST BE STRICTLY INCREASING WITH XBKPT(1) = XLEFT AND XBKPT(NINT+1) = XRIGHT (USED ONLY ON FIRST CALL). EPS = THE RELATIVE TIME ERROR BOUND (USED ONLY ON THE FIRST CALL, UNLESS INDEX = 4). SINGLE STEP ERROR ESTIMATES DIVIDED BY CMAX(I) WILL BE KEPT LESS THAN EPS IN ROOT-MEAN-SQUARE NORM. THE VECTOR CMAX OF WEIGHTS IS COMPUTED IN PDECOL. INITIALLY CMAX(I) IS SET TO ABS(C(I)), WITH A DEFAULT VALUE OF 1 IF ABS(C(I)) .LT. 1. THEREAFTER, CMAX(I) IS THE LARGEST VALUE OF ABS(C(I)) SEEN SO FAR, OR THE INITIAL CMAX(I) IF THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT STATEMENTS 50 AND 130 BELOW. THE USER SHOULD EXERCISE SOME DISCRETION IN CHOOSING EPS. IN GENERAL, THE OVERALL RUNNING TIME FOR A PROBLEM WILL BE GREATER IF EPS IS CHOSEN SMALLER. THERE IS USUALLY LITTLE REASON TO CHOOSE EPS MUCH SMALLER THAN THE ERRORS WHICH ARE BEING INTRODUCED BY THE USERS CHOICE OF THE POLYNOMIAL SPACE. SEE RELATED COMMENTS CONCERNING CMAX BELOW STATEMENT 40. NINT = THE NUMBER OF SUBINTERVALS INTO WHICH THE SPATIAL DOMAIN (XLEFT,XRIGHT) IS TO BE DIVIDED (MUST BE .GE. 1) (USED ONLY ON FIRST CALL). KORD = THE ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED. ITS VALUE MUST BE GREATER THAN 2 AND LESS THAN 21. FOR FIRST ATTEMPTS WE RECOMMEND KORD = 4. IF THE SOLUTION IS SMOOTH AND MUCH ACCURACY IS DESIRED, HIGHER VALUES MAY PROVE TO BE MORE EFFICIENT. WE HAVE SELDOM USED VALUES OF KORD IN EXCESS OF 8 OR 9, THOUGH THEY ARE AVAILABLE FOR USE IN PDECOL (USED ONLY ON FIRST CALL). NCC = THE NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED ON THE APPROXIMATE SOLUTION AT THE BREAKPOINTS IN XBKPT. NCC MUST BE GREATER THAN 1 AND LESS THAN KORD. WE RECOMMEND THE USE OF NCC = 2 SINCE THEORY PREDICTS THAT DRAMATICALLY MORE ACCURATE RESULTS CAN OFTEN BE OBTAINED USING THIS CHOICE (USED ONLY ON FIRST CALL). NPDE = THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS IN THE SYSTEM TO BE SOLVED (USED ONLY ON FIRST CALL). MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS INDEX = 4). ALLOWED VALUES ARE 11, 12, 21, 22. FOR FIRST ATTEMPTS WE RECOMMEND THE USE OF MF = 22. MF HAS TWO DECIMAL DIGITS, METH AND MITER (MF = 10*METH + MITER). METH IS THE BASIC METHOD INDICATOR.. METH = 1 MEANS THE ADAMS METHODS (GENERALIZATIONS OF CRANK-NICOLSON). METH = 2 MEANS THE BACKWARD DIFFERENTIATION FORMULAS (BDF), OR STIFF METHODS OF GEAR. MITER IS THE ITERATION METHOD INDICATOR AND DETERMINES HOW THE JACOBIAN MATRIX IS TO BE COMPUTED.. MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN. FOR THIS USER SUPPLIES SUBROUTINE DERIVF. SEE DESCRIPTION ABOVE. MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED INTERNALLY BY FINITE DIFFERENCES. SEE SUBROUTINES PSETIB AND DIFFF. INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL, WITH THE FOLLOWING VALUES AND MEANINGS.. 1 THIS IS THE FIRST CALL FOR THIS PROBLEM. 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM, AND INTEGRATION IS TO CONTINUE. 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT EXACTLY (NO INTERPOLATION IS DONE). SEE NOTE BELOW. ASSUMES TOUT .GE. THE CURRENT T. IF TOUT IS .LT. THE CURRENT TIME, THEN TOUT IS RESET TO THE CURRENT TIME AND CONTROL IS RETURNED TO THE USER. A CALL TO VALUES WILL PRODUCE ANSWERS FOR THE NEW VALUE OF TOUT. 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING PROGRAM AFTER ONE STEP. TOUT IS IGNORED AND DT MUST BE SET .GT. 0 TO A MAXIMUM ALLOWED DT VALUE. SEE ABOVE. 4 THIS IS NOT THE FIRST CALL FOR THE PROBLEM, AND THE USER HAS RESET EPS AND/OR MF. SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0, IT NEED NOT BE RESET FOR NORMAL CONTINUATION. NOTE.. THE PACKAGE MUST HAVE TAKEN AT LEAST ONE SUCCESSFUL TIME STEP BEFORE A CALL WITH INDEX = 2 OR 4 IS ALLOWED. AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN. ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL. A CHANGE OF PARAMETERS WITH INDEX = 4 CAN BE MADE AFTER EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN PROVIDED AT LEAST ONE SUCCESSFUL TIME STEP HAS BEEN MADE. WORK = FLOATING POINT WORKING ARRAY FOR PDECOL. WE RECOMMEND THAT IT BE INITIALIZED TO ZERO BEFORE THE FIRST CALL TO PDECOL. ITS TOTAL LENGTH MUST BE AT LEAST KORD + 4*NPDE + 9*NPDE**2 + NCPTS*(3*KORD + 2) + NPDE*NCPTS*(3*ML + MAXDER + 7) WHERE... NCPTS = KORD*NINT - NCC*(NINT-1) ML = NPDE*(KORD+IQUAD-1) - 1 IQUAD = 1 IF KORD = 3 AND A NULL BOUNDARY CONDITION EXISTS IQUAD = 0 OTHERWISE MAXDER = 5 UNLESS OTHERWISE SET BY THE USER INTO /OPTION/. IWORK = INTEGER WORKING ARRAY FOR PDECOL. THE FIRST TWO LOCATIONS MUST BE DEFINED AS FOLLOWS... IWORK(1) = LENGTH OF USERS ARRAY WORK IWORK(2) = LENGTH OF USERS ARRAY IWORK THE TOTAL LENGTH OF IWORK MUST BE AT LEAST NCPTS*(NPDE + 1). OUTPUT THE SOLUTION VALUES ARE NOT RETURNED DIRECTLY TO THE USER BY PDECOL. THE METHODS USED IN PDECOL COMPUTE BASIS FUNCTION COEFFICIENTS, SO THE USER (AFTER A RETURN FROM PDECOL) MUST CALL THE PACKAGE ROUTINE VALUES TO OBTAIN HIS APPROXIMATE SOLUTION VALUES AT ANY DESIRED SPACE POINTS X AT THE TIME T = TOUT. SEE THE COMMENTS IN SUBROUTINE VALUES FOR DETAILS ON HOW TO PROPERLY MAKE THE CALL. EXECUTION ERROR MESSAGES WILL BE PRINTED BY PDECOL ON LOGICAL UNIT LOUT WHICH IS THE ONLY VARIABLE IN THE COMMON BLOCK /IOUNIT/. A DEFAULT OF LOUT = 3 IS SET IN THE BLOCK DATA. THE COMMON BLOCK /GEAR0/ CONTAINS THE VARIABLES DTUSED, NQUSED, NSTEP, NFE, AND NJE AND CAN BE ACCESSED EXTERNALLY BY THE USER IF DESIRED. RESPECTIVELY, IT CONTAINS THE STEP SIZE LAST USED (SUCCESS- FULLY), THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN SO FAR, THE NUMBER OF RESIDUAL EVALUATIONS (RES CALLS) SO FAR, AND THE NUMBER OF MATRIX EVALUATIONS (PSETIB CALLS) SO FAR. DIFFUN CALLS ARE COUNTED IN WITH RESIDUAL EVALUATIONS. THE OUTPUT PARAMETERS ARE.. DT = THE STEP SIZE USED LAST, WHETHER SUCCESSFULLY OR NOT. TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL, AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT IS THE CURRENT VALUE OF T TO WHICH THE INTEGRATION HAS BEEN COMPLETED. INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS, WITH THE FOLLOWING VALUES AND MEANINGS.. 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND. -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE ERROR TEST EVEN AFTER REDUCING DT BY A FACTOR OF 1.E10 FROM ITS INITIAL VALUE. -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED. -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE CORRECTOR CONVERGENCE EVEN AFTER REDUCING DT BY A FACTOR OF 1.E10 FROM ITS INITIAL VALUE. -4 SINGULAR MATRIX ENCOUNTERED. PROBABLY DUE TO STORAGE OVERWRITES. -5 INDEX WAS 4 ON INPUT, BUT THE DESIRED CHANGES OF PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT WAS NOT BEYOND T. INTERPOLATION TO T = TOUT WAS PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN, SIMPLY CALL AGAIN WITH INDEX = 4 AND A NEW TOUT. -6 ILLEGAL INDEX VALUE. -7 ILLEGAL EPS VALUE. -8 AN ATTEMPT TO INTEGRATE IN THE WRONG DIRECTION. THE SIGN OF DT IS WRONG RELATIVE TO T0 AND TOUT. -9 DT .EQ. 0.0. -10 ILLEGAL NINT VALUE. -11 ILLEGAL KORD VALUE. -12 ILLEGAL NCC VALUE. -13 ILLEGAL NPDE VALUE. -14 ILLEGAL MF VALUE. -15 ILLEGAL BREAKPOINTS - NOT STRICTLY INCREASING. -16 INSUFFICIENT STORAGE FOR WORK OR IWORK. ----------------------------------------------------------------------- SUBROUTINE VALUES(X,USOL,SCTCH,NDIM1,NDIM2,NPTS,NDERV,WORK) SUBROUTINE VALUES COMPUTES THE SOLUTION U AND THE FIRST NDERV DERIVATIVES OF U AT THE NPTS POINTS X AND AT TIME TOUT AND RETURNS THEM IN THE ARRAY USOL. THIS ROUTINE MUST BE USED TO OBTAIN SOLUTION VALUES SINCE PDECOL DOES NOT RETURN ANY SOLUTION VALUES TO THE USER. SEE PDECOL. THE CALLING PARAMETERS ARE... X = AN ARBITRARY VECTOR OF SPATIAL POINTS OF LENGTH NPTS AT WHICH THE SOLUTION AND THE FIRST NDERV DERIVATIVE VALUES ARE TO BE CALCULATED. IF X .LT. XLEFT ( X .GT. XRIGHT ) THEN THE PIECEWISE POLYNOMIAL OVER THE LEFTMOST ( RIGHT- MOST ) INTERVAL IS EVALUATED TO CALCULATE THE SOLUTION VALUES AT THIS UNUSUAL VALUE OF X. SEE PDECOL. USOL = AN ARRAY WHICH CONTAINS THE SOLUTION AND THE FIRST NDERV DERIVATIVES OF THE SOLUTION AT ALL THE POINTS IN THE INPUT VECTOR X. IN PARTICULAR, USOL(J,I,K) CONTAINS THE VALUE OF THE (K-1)-ST DERIVATIVE OF THE J-TH PDE COMPONENT AT THE I-TH POINT OF THE X VECTOR FOR J = 1 TO NPDE, I = 1 TO NPTS, AND K = 1 TO NDERV+1. SCTCH = A USER SUPPLIED WORKING STORAGE ARRAY OF LENGTH AT LEAST KORD*(NDERV+1). SEE BELOW AND PDECOL FOR DEFINITIONS OF THESE PARAMETERS. NDIM1 = THE FIRST DIMENSION OF THE OUTPUT ARRAY USOL IN THE CALLING PROGRAM. NDIM1 MUST BE .GE. NPDE. NDIM2 = THE SECOND DIMENSION OF THE OUTPUT ARRAY USOL IN THE CALLING PROGRAM. NDIM2 MUST BE .GE. NPTS. NPTS = THE NUMBER OF POINTS IN THE X VECTOR. NDERV = THE NUMBER OF DERIVATIVE VALUES OF THE SOLUTION THAT ARE TO BE CALCULATED. NDERV SHOULD BE LESS THAN KORD SINCE THE KORD-TH DERIVATIVE OF A POLYNOMIAL OF DEGREE KORD-1 IS EQUAL TO ZERO. SEE PDECOL. WORK = THE USERS WORKING STORAGE ARRAY WHICH IS USED IN THIS CASE TO PROVIDE THE CURRENT BASIS FUNCTION COEFFICIENTS AND THE PIECEWISE POLYNOMIAL BREAKPOINT SEQUENCE. PACKAGE ROUTINES CALLED.. BSPLVD,INTERV USER ROUTINES CALLED.. NONE CALLED BY.. USERS MAIN PROGRAM FORTRAN FUNCTIONS USED.. NONE ----------------------------------------------------------------------- ----------------------------------------------------------------------- THIS IS THE MARCH 24, 1978 VERSION OF PDECOL. THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY ANSI-FORTRAN RULES AS WAS CONVENIENTLY POSSIBLE. THE FORTRAN USED VIOLATES ANSI STANDARDS IN THE TWO WAYS LISTED BELOW.... 1. SUBSCRIPTS OF THE GENERAL FORM C*V1 + V2 + V3 ARE USED (POSSIBLY IN A PERMUTED ORDER), WHERE C IS AN INTEGER CONSTANT AND V1, V2, AND V3 ARE INTEGER VARIABLES. 2. ARRAY NAMES APPEAR SINGLY IN DATA STATEMENTS IN THE ROUTINES BSPLVN AND CSTCOL. ----------------------------------------------------------------------- -----------------------------------------------------------------------