*DECK DSS002 SUBROUTINE DSS002(XL,XU,N,U,UX) C... C... SUBROUTINE DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU C... C... ARGUMENT LIST C... C... XL LOWER BOUNDARY VALUE OF X (INPUT) C... C... XU UPPER BOUNDARY VALUE OF X (INPUT) C... C... N NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE C... BOUNDARY POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT C... THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS C... TO BE COMPUTED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL C... VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS C... (OUTPUT) C... C... SUBROUTINE DSS002 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM THE C... CLASSICAL THREE-POINT, SECOND-ORDER FINITE DIFFERENCE APPROXI- C... TIONS C... C... 2 C... U1 = (1/2DX)(-3U1 + 4U2 - U3) + O(DX ) (LEFT BOUNDARY, (1) C... X X = XL) C... C... 2 C... UI = (1/2DX)(UI+1 - UI-1) + O(DX ) (INTERIOR POINT, (2) C... X X NE XL, XU) C... C... 2 C... UN = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX ) (RIGHT BOUNDARY, (3) C... X X = XU) C... C... EQUATIONS (1) TO (3) APPLY OVER A GRID IN X WITH CORRESPONDING C... VALUES OF THE FUNCTION U(X) REPRESENTED AS C... C... U1 U2 U3 UI UN-2 UN-1 UN C... C... X=XL X=XL+DX X=XL+2DX ... X=XI ... X=XU-2DX X=XU-DX X=XU C... C... THE ORIGIN OF EQUATIONS (1) TO (3) IS OUTLINED BELOW. C... C... CONSIDER THE FOLLOWING POLYNOMIAL IN X OF ARBITRARY ORDER C... C... 2 3 C... U(X) = A0 + A1(X - X0) + A2(X - X0) + A3(X - X0) + .... (4) C... C... WE SEEK THE VALUES OF THE COEFFICIENTS A0, A1, A2, ... FOR A C... PARTICULAR FUNCTION U(X). IF X = X0 IS SUBSTITUTED IN EQUATION C... (4), WE HAVE IMMEDIATELY A0 = U(X0). NEXT, IF EQUATION (4) IS C... DIFFERENTIATED WITH RESPECT TO X, C... C... 2 C... DU(X)/DX = U (X) = A1 + 2A2(X - X0) + 3A3(X - X0) + ... (5) C... X C... C... AGAIN, WITH X = X0, A1 = DU(X0)/DX = U (X0). DIFFERENTIATION C... X C... OF EQUATION (5) IN TURN GIVES C... C... D2U(X)/DX2 = U (X) = 2A2 + 6A3(X - X0) + ... C... 2X C... C... AND FOR X = X0, A2 = U (X0)/2F (2F = 1*2, I.E., 2 FACTORIAL). C... 2X C... C... WE CAN CONTINUE THIS PROCESS OF DIFFERENTIATION FOLLOWED BY THE C... SUBSTITUTION X = X0 TO OBTAIN THE SUCCESSIVE COEFFICIENTS IN C... EQUATION (4), A3, A4, ... FINALLY, SUBSTITUTION OF THESE CO- C... EFFICIENTS IN EQUATION (4) GIVES C... C... 2 C... U(X) = U(X0) + U (X0)(X - X0) + U (X0)(X - X0) + C... X 1F 2X 2F C... (6) C... 3 4 C... U (X0)(X - X0) + U (X0)(X - X0) + ... C... 3X 3F 4X 4F C... C... THE CORRESPONDENCE BETWEEN EQUATION (6) AND THE WELL-KNOWN C... TAYLOR SERIES SHOULD BE CLEAR. THUS THE EXPANSION OF A C... FUNCTION, U(X), AROUND A NEIGHBORING POINT X0 IN TERMS OF U(X0) C... AND THE DERIVATIVES OF U(X) AT X = X0 IS EQUIVALENT TO APPROXI- C... MATING U(X) NEAR X0 BY A POLYNOMIAL. C... C... EQUATION (6) IS THE STARTING POINT FOR THE DERIVATION OF THE C... CLASSICAL FINITE DIFFERENCE APPROXIMATIONS OF DERIVATIVES SUCH C... AS THE THREE-POINT FORMULAS OF EQUATIONS (1), (2) AND (3). WE C... WILL NOW CONSIDER THE DERIVATION OF THESE THREE-POINT FORMULAS C... IN A STANDARD FORMAT WHICH CAN THEN BE EXTENDED TO HIGHER C... MULTI-POINT FORMULAS IN OTHER SUBROUTINES, E.G., FIVE-POINT C... FORMULAS IN SUBROUTINE DSS004. C... C... THREE-POINT FORMULAS C... C... (1) LEFT END, POINT I = 1 C... C... IF EQUATION (6) IS WRITTEN AROUND THE POINT X = XL FOR X = XL + C... DX AND X = XL + 2DX, FOR WHICH THE CORRESPONDING VALUES OF U(X) C... ARE U1, U2 AND U3 (U1 AND U2 ARE SEPARATED WITH RESPECT TO X BY C... DISTANCE DX AS ARE U2 AND U3, I.E., WE ASSUME A UNIFORM GRID C... SPACING, DX, FOR INDEPENDENT VARIABLE X) C... C... 2 3 C... U2 = U1 + U1 ( DX) + U1 ( DX) + U1 ( DX) + ... (7) C... X 1F 2X 2F 3X 3F C... C... 2 3 C... U3 = U1 + U1 (2DX) + U1 (2DX) + U1 (2DX) + ... (8) C... X 1F 2X 2F 3X 3F C... C... WE CAN NOW TAKE A LINEAR COMBINATION OF EQUATIONS (7) AND (8) C... BY FIRST MULTIPLYING EQUATION (7) BY A CONSTANT, A, AND EQUA- C... TION (8) BY CONSTANT B C... C... 2 3 C... A(U2 = U1 + U1 ( DX) + U1 ( DX) + U1 ( DX) + ...) (9) C... X 1F 2X 2F 3X 3F C... C... 2 3 C... B(U3 = U1 + U1 (2DX) + U1 (2DX) + U1 (2DX) + ...) (10) C... X 1F 2X 2F 3X 3F C... C... CONSTANTS A AND B ARE THEN SELECTED SO THAT THE COEFFICIENTS OF C... THE U1 TERMS SUM TO ONE (SINCE WE ARE INTERESTED IN OBTAINING C... X C... A FINITE DIFFERENCE APPROXIMATION FOR THIS FIRST DERIVATIVE). C... ALSO, WE SELECT A AND B SO THAT THE COEFFICIENTS OF THE U1 C... 2X C... TERMS SUM TO ZERO IN ORDER TO DROP OUT THE CONTRIBUTION OF THIS C... SECOND DERIVATIVE (THE BASIC IDEA IS TO DROP OUT AS MANY OF THE C... DERIVATIVES AS POSSIBLE IN THE TAYLOR SERIES BEYOND THE DERI- C... VATIVE OF INTEREST, IN THIS CASE U1 , IN ORDER TO PRODUCE A C... X C... FINITE DIFFERENCE APPROXIMATION FOR THE DERIVATIVE OF MAXIMUM C... ACCURACY). IN THIS CASE WE HAVE ONLY TWO CONSTANTS, A AND B, C... TO SELECT SO WE CAN DROP OUT ONLY THE SECOND DERIVATIVE, U1 , C... 2X C... IN THE TAYLOR SERIES (IN ADDITION TO RETAINING THE FIRST DERI- C... VATIVE). THIS PROCEDURE LEADS TO TWO LINEAR ALGEBRAIC EQUA- C... TIONS IN THE TWO CONSTANTS C... C... A + 2B = 1 C... C... A + 4B = 0 C... C... SOLUTION OF THESE EQUATIONS FOR A AND B GIVES C... C... A = 2, B = -1/2 C... C... SOLUTION OF EQUATIONS (9) AND (10) FOR U1 WITH THESE VALUES OF C... A AND B GIVES EQUATION (1) X C... C... 2 C... U1 = (1/2DX)(-3U1 + 4U2 - U3) + O(DX ) (1) C... X C... 2 C... THE TERM O(DX ) INDICATES A PRINCIPAL ERROR TERM DUE TO TRUNCA- C... 2 C... TION OF THE TAYLOR SERIES WHICH IS OF ORDER DX . THIS TERM IN C... 2 C... FACT EQUALS U1 DX /3F, WHICH IS EASILY OBTAINED IN DERIVING C... 3X C... EQUATION (1). C... C... THIS SAME BASIC PROCEDURE CAN NOW BE APPLIED TO THE DERIVATION C... OF EQUATIONS (2) AND (3). C... C... (2) INTERIOR POINT I C... C... 2 3 C... A(UI-1 = UI + UI (-DX) + UI (-DX) + UI (-DX) + ...) C... X 1F 2X 2F 3X 3F C... C... 2 3 C... B(UI+1 = UI + UI ( DX) + UI ( DX) + UI ( DX) + ...) C... X 1F 2X 2F 3X 3F C... C... -A + B = 1 C... C... A + B = 0 C... C... A = -1/2, B = 1/2 C... 2 C... UI = (1/2DX)(UI+1 - UI-1) + O(DX ) (2) C... X C... C... (3) RIGHT END, POINT I = N C... C... 2 3 C... A(UN-2 = UN + UN (-2DX) + UN (-2DX) + UN (-2DX) + ...) C... X 1F 2X 2F 3X 3F C... C... 2 3 C... B(UN-1 = UN + UN ( -DX) + UN ( -DX) + UN ( -DX) + ...) C... X 1F 2X 2F 3X 3F C... C... -2A - B = 1 C... C... 4A + B = 0 C... C... A = -2, B = 1/2 C... 2 C... UN = (1/2DX)(3UN - 4UN-1 + UN-2) + O(DX ) (3) C... X C... C... THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1), (2) AND (3) CAN C... BE SUMMARIZED AS C... C... -3 4 -1 C... C... 1/2 -1 0 1 C... C... 1 -4 3 C... C... WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 2, M = C... 1, P = 0, 1, 2 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER- C... ENTIATION, MATH. GAZ., VOL. 25, 1941). C... C... EQUATIONS (1), (2) AND (3) CAN NOW BE PROGRAMMED TO GENERATE C... THE DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF C... X C... SUBROUTINE DSS002, RESPECTIVELY). C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R2FDX, U, UX, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT DX=(XU-XL)/DFLOAT(N-1) R2FDX=1.D+00/(2.D+00*DX) NM1=N-1 C... C... EQUATION (1) (NOTE - THE RHS OF THE FINITE DIFFERENCE APPROXI- C... TIONS, EQUATIONS (1), (2) AND (3) HAVE BEEN FORMATTED SO THAT C... THE NUMERICAL WEIGHTING COEFFICIENTS CAN BE MORE EASILY ASSOCI- C... ATED WITH THE BICKLEY MATRIX LISTED ABOVE) UX(1)=R2FDX* 1( -3.D+00*U( 1) +4.D+00*U( 2) -1.D+00*U( 3)) C... C... EQUATION (2) DO 1 I=2,NM1 UX(I)=R2FDX* 1( -1.D+00*U(I-1) +0.D+00*U( I) +1.D+00*U(I+1)) 1 CONTINUE C... C... EQUATION (3) UX(N)=R2FDX* 1( 1.D+00*U(N-2) -4.D+00*U(N-1) +3.D+00*U( N)) RETURN END *DECK DSS004 SUBROUTINE DSS004(XL,XU,N,U,UX) C... C... SUBROUTINE DSS004 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL C... FIVE-POINT, FOURTH-ORDER FINITE DIFFERENCE APPROXIMATIONS C... C... ARGUMENT LIST C... C... XL LOWER BOUNDARY VALUE OF X (INPUT) C... C... XU UPPER BOUNDARY VALUE OF X (INPUT) C... C... N NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE C... BOUNDARY POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT C... THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS C... TO BE COMPUTED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL C... VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS C... (OUTPUT) C... C... THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR C... POLYNOMIALS) ARE GIVEN IN SUBROUTINE DSS002. C... C... FIVE-POINT FORMULAS C... C... (1) LEFT END, POINT I = 1 C... C... 2 3 4 C... A(U2 = U1 + U1 ( DX) + U1 ( DX) + U1 ( DX) + U1 ( DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 ( DX) + U1 ( DX) + U1 ( DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... B(U3 = U1 + U1 (2DX) + U1 (2DX) + U1 (2DX) + U1 (2DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (2DX) + U1 (2DX) + U1 (2DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... C(U4 = U1 + U1 (3DX) + U1 (3DX) + U1 (3DX) + U1 (3DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (3DX) + U1 (3DX) + U1 (3DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... D(U5 = U1 + U1 (4DX) + U1 (4DX) + U1 (4DX) + U1 (4DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (4DX) + U1 (4DX) + U1 (4DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... CONSTANTS A, B, C AND D ARE SELECTED SO THAT THE COEFFICIENTS C... OF THE U1 TERMS SUM TO ONE AND THE COEFFICIENTS OF THE U1 , C... X 2X C... U1 AND U1 TERMS SUM TO ZERO C... 3X 4X C... C... A + 2B + 3C + 4D = 1 C... C... A + 4B + 9C + 16D = 0 C... C... A + 8B + 27C + 64D = 0 C... C... A + 16B + 81C + 256D = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU- C... TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U C... 4X C... TERMS, FOR U1 GIVES THE FOLLOWING FIVE-POINT APPROXIMATION C... X C... 4 C... U1 = (1/12DX)(-25U1 + 48U2 - 36U3 + 16U4 - 3U5) + O(DX ) (1) C... X C... C... (2) INTERIOR POINT, I = 2 C... C... 2 3 4 C... A(U1 = U2 + U2 (-DX) + U2 (-DX) + U2 (-DX) + U2 (-DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U2 (-DX) + U2 (-DX) + U2 (-DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... B(U3 = U2 + U2 ( DX) + U2 ( DX) + U2 ( DX) + U2 ( DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U2 ( DX) + U2 ( DX) + U2 ( DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... C(U4 = U2 + U2 (2DX) + U2 (2DX) + U2 (2DX) + U2 (2DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U2 (2DX) + U2 (2DX) + U2 (2DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... D(U5 = U2 + U2 (3DX) + U2 (3DX) + U2 (3DX) + U2 (3DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U2 (3DX) + U2 (3DX) + U2 (3DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... -A + B + 2C + 3D = 1 C... C... A + B + 4C + 9D = 0 C... C... -A + B + 8C + 27D = 0 C... C... A + B + 16C + 81D = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU- C... TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U C... 4X C... TERMS, FOR U1 GIVES THE FOLLOWING FIVE-POINT APPROXIMATION C... X C... 4 C... U2 = (1/12DX)(-3U1 - 10U2 + 18U3 - 6U4 + U5) + O(DX ) (2) C... X C... C... (3) INTERIOR POINT I, I NE 2, N-1 C... C... 2 3 C... A(UI-2 = UI + UI (-2DX) + UI (-2DX) + UI (-2DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UI (-2DX) + UI (-2DX) + UI (-2DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... B(UI-1 = UI + UI ( -DX) + UI ( -DX) + UI ( -DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UI ( -DX) + UI ( -DX) + UI ( -DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... C(UI+1 = UI + UI ( DX) + UI ( DX) + UI ( DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UI ( DX) + UI ( DX) + UI ( DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... D(UI+2 = UI + UI ( 2DX) + UI ( 2DX) + UI ( 2DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UI ( 2DX) + UI ( 2DX) + UI ( 2DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... -2A - B + C + 2D = 1 C... C... 4A + B + C + 4D = 0 C... C... -8A - B + C + 8D = 0 C... C... 16A + B + C + 16D = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU- C... TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U C... 4X C... TERMS, FOR U1 GIVES THE FOLLOWING FIVE-POINT APPROXIMATION C... X C... 4 C... UI = (1/12DX)(UI-2 - 8UI-1 + 0UI + 8UI+1 - UI+2) + O(DX ) (3) C... X C... C... (4) INTERIOR POINT, I = N-1 C... C... 2 3 C... A(UN-4 = UN-1 + UN-1 (-3DX) + UN-1 (-3DX) + UN-1 (-3DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN-1 (-3DX) + UN-1 (-3DX) + UN-1 (-3DX) + ... C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... B(UN-3 = UN-1 + UN-1 (-2DX) + UN-1 (-2DX) + UN-1 (-2DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN-1 (-2DX) + UN-1 (-2DX) + UN-1 (-2DX) + ... C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... C(UN-2 = UN-1 + UN-1 ( -DX) + UN-1 (- -X) + UN-1 ( -DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN-1 ( -DX) + UN-1 ( -DX) + UN-1 ( -DX) + ... C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... D(UN = UN-1 + UN-1 ( DX) + UN-1 ( DX) + UN-1 ( DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN-1 ( DX) + UN-1 ( DX) + UN-1 ( DX) + ... C... 4X 4F 5X 5F 6X 6F C... C... -3A - 2B - C + D = 1 C... C... 9A + 4B + C + D = 0 C... C... -27A - 8B - C + D = 0 C... C... 81A + 16B + C + D = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU- C... TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U C... 4X C... TERMS, FOR U1 GIVES THE FOLLOWING FIVE-POINT APPROXIMATION C... X C... 4 C... UN-1 = (1/12DX)(-UN-4 + 6UN-3 - 18UN-2 + 10UN-1 + 3UN) + O(DX ) C... X C... (4) C... C... (5) RIGHT END, POINT I = N C... C... 2 3 C... A(UN-4 = UN + UN (-4DX) + UN (-4DX) + UN (-4DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN (-4DX) + UN (-4DX) + UN (-4DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... B(UN-3 = UN + UN (-3DX) + UN (-3DX) + UN (-3DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN (-3DX) + UN (-3DX) + UN (-3DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... C(UN-2 = UN + UN (-2DX) + UN (-2DX) + UN (-2DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN (-2DX) + UN (-2DX) + UN (-2DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... 2 3 C... D(UN-1 = UN + UN ( -DX) + UN ( -DX) + UN ( -DX) C... X 1F 2X 2F 3X 3F C... C... 4 5 6 C... + UN ( -DX) + UN ( -DX) + UN ( -DX) + ...) C... 4X 4F 5X 5F 6X 6F C... C... -4A - 3B - 2C - D = 1 C... C... 16A + 9B + 4C + D = 0 C... C... -64A - 27B - 8C - D = 0 C... C... 256A + 81B + 16C + D = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C AND D FOLLOWED BY THE SOLU- C... TION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE U C... 4X C... TERMS, FOR U1 GIVES THE FOLLOWING FIVE-POINT APPROXIMATION C... X C... 4 C... UN = (1/12DX)(3UN-4 - 16UN-3 + 36UN-2 - 48UN-1 + 25UN) + O(DX ) C... X C... (5) C... C... THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1) TO (5) CAN BE C... SUMMARIZED AS C... C... -25 48 -36 16 -3 C... C... -3 -10 18 -6 1 C... C... 1/12 1 -8 0 8 -1 C... C... -1 6 -18 10 3 C... C... 3 -16 36 -48 25 C... C... WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 4, M = C... 1, P = 0, 1, 2, 3, 4 (BICKLEY, W. G., FORMULAE FOR NUMERICAL C... DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941. NOTE - THE BICKLEY C... COEFFICIENTS HAVE BEEN DIVIDED BY A COMMON FACTOR OF TWO). C... C... EQUATIONS (1) TO (5) CAN NOW BE PROGRAMMED TO GENERATE THE C... DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF SUB- C... X C... ROUTINE DSS004 RESPECTIVELY). C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R4FDX, U, UX, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT DX=(XU-XL)/DFLOAT(N-1) R4FDX=1.D+00/(12.D+00*DX) NM2=N-2 C... C... EQUATION (1) (NOTE - THE RHS OF EQUATIONS (1), (2), (3), (4) C... AND (5) HAVE BEEN FORMATTED SO THAT THE NUMERICAL WEIGHTING C... COEFFICIENTS CAN BE MORE EASILY ASSOCIATED WITH THE BICKLEY C... MATRIX ABOVE) UX( 1)=R4FDX* 1(-25.D+00 *U( 1) 2 +48.D+00 *U( 2) 3 -36.D+00 *U( 3) 4 +16.D+00 *U( 4) 5 -3.D+00 *U( 5)) C... C... EQUATION (2) UX( 2)=R4FDX* 1( -3.D+00 *U( 1) 2 -10.D+00 *U( 2) 3 +18.D+00 *U( 3) 4 -6.D+00 *U( 4) 5 +1.D+00 *U( 5)) C... C... EQUATION (3) DO 1 I=3,NM2 UX( I)=R4FDX* 1( +1.D+00 *U(I-2) 2 -8.D+00 *U(I-1) 3 +0.D+00 *U( I) 4 +8.D+00 *U(I+1) 5 -1.D+00 *U(I+2)) 1 CONTINUE C... C... EQUATION (4) UX(N-1)=R4FDX* 1( -1.D+00 *U(N-4) 2 +6.D+00 *U(N-3) 3 -18.D+00 *U(N-2) 4 +10.D+00 *U(N-1) 5 +3.D+00 *U(N )) C... C... EQUATION (5) UX( N)=R4FDX* 1( +3.D+00 *U(N-4) 2 -16.D+00 *U(N-3) 3 +36.D+00 *U(N-2) 4 -48.D+00 *U(N-1) 5 +25.D+00 *U(N )) RETURN END *DECK DSS006 SUBROUTINE DSS006(XL,XU,N,U,UX) C... C... SUBROUTINE DSS006 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL C... SEVEN-POINT, SIXTH-ORDER FINITE DIFFERENCE APPROXIMATIONS C... C... ARGUMENT LIST C... C... XL LOWER BOUNDARY VALUE OF X (INPUT) C... C... XU UPPER BOUNDARY VALUE OF X (INPUT) C... C... N NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE C... BOUNDARY POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT C... THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS C... TO BE COMPUTED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL C... VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS C... (OUTPUT) C... C... THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR C... POLYNOMIALS) ARE GIVEN IN SUBROUTINES DSS002 AND DSS004. C... C... SEVEN-POINT FORMULAS C... C... (1) LEFT END, POINT I = 1 C... C... 2 3 4 C... A(U2 = U1 + U1 ( DX) + U1 ( DX) + U1 ( DX) + U1 ( DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 ( DX) + U1 ( DX) + U1 ( DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... B(U3 = U1 + U1 (2DX) + U1 (2DX) + U1 (2DX) + U1 (2DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (2DX) + U1 (2DX) + U1 (2DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... C(U4 = U1 + U1 (3DX) + U1 (3DX) + U1 (3DX) + U1 (3DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (3DX) + U1 (3DX) + U1 (3DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... D(U5 = U1 + U1 (4DX) + U1 (4DX) + U1 (4DX) + U1 (4DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (4DX) + U1 (4DX) + U1 (4DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... E(U6 = U1 + U1 (5DX) + U1 (5DX) + U1 (5DX) + U1 (5DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (5DX) + U1 (5DX) + U1 (5DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... F(U7 = U1 + U1 (6DX) + U1 (6DX) + U1 (6DX) + U1 (6DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (6DX) + U1 (6DX) + U1 (6DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... CONSTANTS A, B, C, D, E AND F ARE SELECTED SO THAT THE COEFFI- C... CIENTS OF THE U1 TERMS SUM TO ONE AND THE COEFFICIENTS OF C... X C... THE U1 , U1 , U1 , U1 AND U1 TERMS SUM TO ZERO. C... 2X 3X 4X 5X 6X C... C... 1 1 1 1 1 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 1 C... C... 2 2 2 2 2 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 0 C... C... 3 3 3 3 3 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 0 C... C... 4 4 4 4 4 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 0 C... C... 5 5 5 5 5 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 0 C... C... 6 6 6 6 6 C... A + 2 B + 3 C + 4 D + 5 E + 6 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR U1 GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... U1 = (1/6F)(-1764*U1 + 4320U2 - 5400U3 + 4800U4 C... X 6 (1) C... + 2700U5 + 864U6 - 120U7) + O(DX ) C... C... (2) INTERIOR POINT, I = 2 C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = 2 C... AS C... C... A(U1 = U2 + ...) C... C... B(U3 = U2 + ...) C... C... C(U4 = U2 + ...) C... C... D(U5 = U2 + ...) C... C... E(U6 = U2 + ...) C... C... F(U7 = U2 + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 1 C... C... 2 2 2 2 2 2 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 0 C... C... 3 3 3 3 3 3 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 0 C... C... 4 4 4 4 4 4 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 0 C... C... 5 5 5 5 5 5 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 0 C... C... 6 6 6 6 6 6 C... -1 A + 1 B + 2 C + 3 D + 4 E + 5 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR U2 GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... U2 = (1/6F)(-120U1 - 924U2 + 1800U3 - 1200U4 C... X 6 (2) C... + 600U5 - 180U6 + 24U7) + O(DX ) C... C... (3) INTERIOR POINT, I = 3 C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = 3 C... AS C... C... A(U1 = U3 + ...) C... C... B(U2 = U3 + ...) C... C... C(U4 = U3 + ...) C... C... D(U5 = U3 + ...) C... C... E(U6 = U3 + ...) C... C... F(U7 = U3 + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 1 C... C... 2 2 2 2 2 2 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 0 C... C... 3 3 3 3 3 3 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 0 C... C... 4 4 4 4 4 4 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 0 C... C... 5 5 5 5 5 5 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 0 C... C... 6 6 6 6 6 6 C... -2 A + -1 B + 1 C + 2 D + 3 E + 4 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR U3 GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... U3 = (1/6F)(24U1 - 288U2 - 420U3 + 960U4 C... X 6 (3) C... - 360U5 + 96U6 - 12U7) + O(DX ) C... C... (4) INTERIOR POINT, I NE 2, 3, N-2, N-1 C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = I C... AS C... C... A(UI-3 = UI + ...) C... C... B(UI-2 = UI + ...) C... C... C(UI-1 = UI + ...) C... C... D(UI+1 = UI + ...) C... C... E(UI+2 = UI + ...) C... C... F(UI+3 = UI + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 1 C... C... 2 2 2 2 2 2 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 0 C... C... 3 3 3 3 3 3 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 0 C... C... 4 4 4 4 4 4 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 0 C... C... 5 5 5 5 5 5 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 0 C... C... 6 6 6 6 6 6 C... -3 A + -2 B + -1 C + 1 D + 2 E + 3 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR UI GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... UI = (1/6F)(-12UI-3 + 108UI-2 - 540UI-1 + 0UI C... X 6 (4) C... + 540UI+1 - 108UI+2 + 12UI+3) + O(DX ) C... C... (5) INTERIOR POINT, I = N-2 C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N-2 C... AS C... C... A(UN-6 = UN-2 + ...) C... C... B(UN-5 = UN-2 + ...) C... C... C(UN-4 = UN-2 + ...) C... C... D(UN-3 = UN-2 + ...) C... C... E(UN-1 = UN-2 + ...) C... C... F(UN = UN-2 + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 1 C... C... 2 2 2 2 2 2 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 0 C... C... 3 3 3 3 3 3 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 0 C... C... 4 4 4 4 4 4 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 0 C... C... 5 5 5 5 5 5 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 0 C... C... 6 6 6 6 6 6 C... -4 A + -3 B + -2 C + -1 D + 1 E + 2 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR UN-2 GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... UN-2 = (1/6F)(12UN-6 - 96UN-5 + 360UN-4 - 960UN-3 C... X 6 (5) C... + 420UN-2 + 288UN-1 - 24UN) + O(DX ) C... C... (6) INTERIOR POINT, I = N-1 C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N-1 C... AS C... C... A(UN-6 = UN-1 + ...) C... C... B(UN-5 = UN-1 + ...) C... C... C(UN-4 = UN-1 + ...) C... C... D(UN-3 = UN-1 + ...) C... C... E(UN-2 = UN-1 + ...) C... C... F(UN = UN-1 + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 1 C... C... 2 2 2 2 2 2 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 0 C... C... 3 3 3 3 3 3 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 0 C... C... 4 4 4 4 4 4 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 0 C... C... 5 5 5 5 5 5 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 0 C... C... 6 6 6 6 6 6 C... -5 A + -4 B + -3 C + -2 D + -1 E + 1 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR UN-1 GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... UN-1 = (1/6F)(-24UN-6 + 180UN-5 - 600UN-4 + 1200UN-3 C... X 6 (6) C... - 1800UN-2 + 924UN-1 + 120UN) + O(DX ) C... C... (7) RIGHT END, POINT I = N C... C... THE PRECEDING TAYLOR SERIES CAN BE SUMMARIZED FOR POINT I = N C... AS C... C... A(UN-6 = UN + ...) C... C... B(UN-5 = UN + ...) C... C... C(UN-4 = UN + ...) C... C... D(UN-3 = UN + ...) C... C... E(UN-2 = UN + ...) C... C... F(UN-1 = UN + ...) C... C... THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... 1 1 1 1 1 1 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 1 C... C... 2 2 2 2 2 2 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0 C... C... 3 3 3 3 3 3 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0 C... C... 4 4 4 4 4 4 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0 C... C... 5 5 5 5 5 5 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0 C... C... 6 6 6 6 6 6 C... -6 A + -5 B + -4 C + -3 D + -2 E + -1 F = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E AND F FOLLOWED BY THE C... SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR UN GIVES THE FOLLOWING SEVEN-POINT APPROXI- C... 6X X C... TION C... C... UN = (1/6F)(120UN-6 - 864UN-5 + 2700UN-4 - 4800UN-3 C... X 6 (7) C... + 5400UN-2 - 4320UN-1 + 1764UN) + O(DX ) C... C... THE WEIGHTING COEFFICIENTS FOR EQUATIONS (1) TO (7) CAN BE C... SUMMARIZED AS C... C... -1764 4320 -5400 4800 -2700 864 -120 C... C... -120 -924 1800 -1200 600 -180 24 C... C... 24 -288 -420 960 -360 96 -12 C... C... 1/6F -12 108 -540 0 540 -108 12 C... C... 12 -96 360 -960 420 288 -24 C... C... -24 180 -600 1200 -1800 924 120 C... C... 120 -864 2700 -4800 5400 -4320 1764 C... C... WHICH ARE THE COEFFICIENTS REPORTED BY BICKLEY FOR N = 6, M = C... 1, P = 0 TO 6 (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFER- C... ENTIATION, MATH. GAZ., VOL. 25, 1941). C... C... EQUATIONS (1) TO (7) CAN NOW BE PROGRAMMED TO GENERATE THE C... DERIVATIVE U (X) OF FUNCTION U(X) (ARGUMENTS U AND UX OF SUB- C... X C... ROUTINE DSS006, RESPECTIVELY). C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R6FDX, U, UX, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT DX=(XU-XL)/DFLOAT(N-1) R6FDX=1.D+00/(720.D+00*DX) NM3=N-3 C... C... EQUATION (1) UX( 1)=R6FDX* 1( -1764.D+00 *U( 1) 2 +4320.D+00 *U( 2) 3 -5400.D+00 *U( 3) 4 +4800.D+00 *U( 4) 5 -2700.D+00 *U( 5) 6 +864.D+00 *U( 6) 7 -120.D+00 *U( 7)) C... C... EQUATION (2) UX( 2)=R6FDX* 1( -120.D+00 *U( 1) 2 -924.D+00 *U( 2) 3 +1800.D+00 *U( 3) 4 -1200.D+00 *U( 4) 5 +600.D+00 *U( 5) 6 -180.D+00 *U( 6) 7 +24.D+00 *U( 7)) C... C... EQUATION (3) UX( 3)=R6FDX* 1( +24.D+00 *U( 1) 2 -288.D+00 *U( 2) 3 -420.D+00 *U( 3) 4 +960.D+00 *U( 4) 5 -360.D+00 *U( 5) 6 +96.D+00 *U( 6) 7 -12.D+00 *U( 7)) C... C... EQUATION (4) DO 1 I=4,NM3 UX( I)=R6FDX* 1( -12.D+00 *U(I-3) 2 +108.D+00 *U(I-2) 3 -540.D+00 *U(I-1) 4 +0.D+00 *U(I ) 5 +540.D+00 *U(I+1) 6 -108.D+00 *U(I+2) 7 +12.D+00 *U(I+3)) 1 CONTINUE C... C... EQUATION (5) UX(N-2)=R6FDX* 1( +12.D+00 *U(N-6) 2 -96.D+00 *U(N-5) 3 +360.D+00 *U(N-4) 4 -960.D+00 *U(N-3) 5 +420.D+00 *U(N-2) 6 +288.D+00 *U(N-1) 7 -24.D+00 *U(N )) C... C... EQUATION (6) UX(N-1)=R6FDX* 1( -24.D+00 *U(N-6) 2 +180.D+00 *U(N-5) 3 -600.D+00 *U(N-4) 4 +1200.D+00 *U(N-3) 5 -1800.D+00 *U(N-2) 6 +924.D+00 *U(N-1) 7 +120.D+00 *U(N )) C... C... EQUATION (7) UX( N)=R6FDX* 1( +120.D+00 *U(N-6) 2 -864.D+00 *U(N-5) 3 +2700.D+00 *U(N-4) 4 -4800.D+00 *U(N-3) 5 +5400.D+00 *U(N-2) 6 -4320.D+00 *U(N-1) 7 +1764.D+00 *U(N )) RETURN END *DECK DSS008 SUBROUTINE DSS008(XL,XU,N,U,UX) C... C... SUBROUTINE DSS008 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL C... NINE-POINT, EIGHTH-ORDER FINITE DIFFERENCE APPROXIMATIONS C... C... ARGUMENT LIST C... C... XL LOWER BOUNDARY VALUE OF X (INPUT) C... C... XU UPPER BOUNDARY VALUE OF X (INPUT) C... C... N NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE C... BOUNDARY POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT C... THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS C... TO BE COMPUTED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL C... VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS C... (OUTPUT) C... C... THE MATHEMATICAL DETAILS OF THE FOLLOWING TAYLOR SERIES (OR C... POLYNOMIALS) ARE GIVEN IN SUBROUTINES DSS002, DSS004 AND C... DSS006. C... C... NINE-POINT FORMULAS C... C... (1) LEFT END, POINT I = 1 C... C... 2 3 4 C... A(U2 = U1 + U1 ( DX) + U1 ( DX) + U1 ( DX) + U1 ( DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 ( DX) + U1 ( DX) + U1 ( DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... B(U3 = U1 + U1 (2DX) + U1 (2DX) + U1 (2DX) + U1 (2DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (2DX) + U1 (2DX) + U1 (2DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... C(U4 = U1 + U1 (3DX) + U1 (3DX) + U1 (3DX) + U1 (3DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (3DX) + U1 (3DX) + U1 (3DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... D(U5 = U1 + U1 (4DX) + U1 (4DX) + U1 (4DX) + U1 (4DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (4DX) + U1 (4DX) + U1 (4DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... E(U6 = U1 + U1 (5DX) + U1 (5DX) + U1 (5DX) + U1 (5DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (5DX) + U1 (5DX) + U1 (5DX) + ...) C... 5X 5F 6X 6F 7X 7F C... 2 3 4 C... F(U7 = U1 + U1 (6DX) + U1 (6DX) + U1 (6DX) + U1 (6DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (6DX) + U1 (6DX) + U1 (6DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... G(U8 = U1 + U1 (7DX) + U1 (7DX) + U1 (7DX) + U1 (7DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (7DX) + U1 (7DX) + U1 (7DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... 2 3 4 C... H(U9 = U1 + U1 (8DX) + U1 (8DX) + U1 (8DX) + U1 (8DX) C... X 1F 2X 2F 3X 3F 4X 4F C... C... 5 6 7 C... + U1 (8DX) + U1 (8DX) + U1 (8DX) + ...) C... 5X 5F 6X 6F 7X 7F C... C... CONSTANTS A, B, C, D, E, F, G AND H ARE SELECTED SO THAT THE C... COEFFICIENTS OF THE U1 TERMS SUM TO ONE AND THE COEFFICIENTS OF C... X C... THE U1 , U1 , U1 , U1 , U1 , U1 AND U1 SUM TO ZERO. C... 2X 3X 4X 5X 6X, 7X 8X C... C... 1 1 1 1 1 1 1 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 1 C... C... 2 2 2 2 2 2 2 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 3 3 3 3 3 3 3 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 4 4 4 4 4 4 4 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 5 5 5 5 5 5 5 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 6 6 6 6 6 6 6 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 7 7 7 7 7 7 7 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... 8 8 8 8 8 8 8 C... A + 2 B + 3 C + 4 D + 5 E + 6 F + 7 G + 8 H = 0 C... C... SIMULTANEOUS SOLUTION FOR A, B, C, D, E, F, G AND H FOLLOWED BY C... THE SOLUTION OF THE PRECEDING TAYLOR SERIES, TRUNCATED AFTER THE C... U TERMS, FOR U1 GIVES THE FOLLOWING NINE-POINT APPROXIMATION C... 8X X C... C... U1 = (1/8F)(-109584U1 + 322560U2 - 564480U3 + 752640U4 C... X C... -705600U5 + 451584U6 - 188160U7 + 46080U8 C... 8 C... - 5040U9) + O(DX ) (1) C... C... THE PRECEDING ANALYSIS CAN BE REPEATED TO PRODUCE NINE-POINT C... APPROXIMATIONS FOR THE FIRST DERIVATIVES U2 , U3 , U4 , UI , C... X X X X C... UN-3 , UN-2 , UN-1 AND UN . THE RESULTS CAN BE SUMMARIZED BY C... X X X X C... THE FOLLOWING BICKLEY MATRIX FOR N = 8, M = 1, P = 0 TO 8, C... (BICKLEY, W. G., FORMULAE FOR NUMERICAL DIFFERENTIATION, MATH. C... GAZ., VOL. 25, 1941) C... C... -109584 322560 -564480 752640 -705600 C... C... -5040 -64224 141120 -141120 117600 C... C... 720 -11520 -38304 80640 -50400 C... C... -240 2880 -20160 -18144 50400 C... C... 1/8F 144 -1536 8064 -32256 0 C... C... -144 1440 -6720 20160 -50400 C... C... 240 -2304 10080 -26880 50400 C... C... -720 6720 -28224 70560 -117600 C... C... 5040 -46080 188160 -451584 705600 C... C... FROM THIS BICKLEY MATRIX, THE FINITE DIFFERENCE APPROXIMATION OF C... THE FIRST DERIVATIVE CAN BE PROGRAMMED FOR EACH OF THE GRID POINTS C... 1, 2, 3, 4,..., I,..., N-3, N-2, N-1, N (TAKING INTO ACCOUNT C... THE SYMMETRY PROPERTIES OF THE BICKLEY MATRIX). C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R8FDX, U, UX, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT DX=(XU-XL)/DFLOAT(N-1) R8FDX=1.D+00/(40320.D+00*DX) NM4=N-4 C... C... GRID POINT 1 UX( 1)=R8FDX* 1 (-109584.D+00 *U( 1) 2 +322560.D+00 *U( 2) 3 -564480.D+00 *U( 3) 4 +752640.D+00 *U( 4) 5 -705600.D+00 *U( 5) 6 +451584.D+00 *U( 6) 7 -188160.D+00 *U( 7) 8 +46080.D+00 *U( 8) 9 -5040.D+00 *U( 9)) C... C... GRID POINT 2 UX( 2)=R8FDX* 1 (-5040.D+00 *U( 1) 2 -64224.D+00 *U( 2) 3 +141120.D+00 *U( 3) 4 -141120.D+00 *U( 4) 5 +117600.D+00 *U( 5) 6 -70560.D+00 *U( 6) 7 +28224.D+00 *U( 7) 8 -6720.D+00 *U( 8) 9 +720.D+00 *U( 9)) C... C... GRID POINT 3 UX( 3)=R8FDX* 1 (+720.D+00 *U( 1) 2 -11520.D+00 *U( 2) 3 -38304.D+00 *U( 3) 4 +80640.D+00 *U( 4) 5 -50400.D+00 *U( 5) 6 +26880.D+00 *U( 6) 7 -10080.D+00 *U( 7) 8 +2304.D+00 *U( 8) 9 -240.D+00 *U( 9)) C... C... GRID POINT 4 UX( 4)=R8FDX* 1 (-240.D+00 *U( 1) 2 +2880.D+00 *U( 2) 3 -20160.D+00 *U( 3) 4 -18144.D+00 *U( 4) 5 +50400.D+00 *U( 5) 6 -20160.D+00 *U( 6) 7 +6720.D+00 *U( 7) 8 -1440.D+00 *U( 8) 9 +144.D+00 *U( 9)) C... C... GRID POINT I DO 1 I=5,NM4 UX( I)=R8FDX* 1 (+144.D+00 *U(I-4) 2 -1536.D+00 *U(I-3) 3 +8064.D+00 *U(I-2) 4 -32256.D+00 *U(I-1) 5 +0.D+00 *U(I ) 6 +32256.D+00 *U(I+1) 7 -8064.D+00 *U(I+2) 8 +1536.D+00 *U(I+3) 9 -144.D+00 *U(I+4)) 1 CONTINUE C... C... GRID POINT N-3 UX(N-3)=R8FDX* 1 (-144.D+00 *U(N-8) 2 +1440.D+00 *U(N-7) 3 -6720.D+00 *U(N-6) 4 +20160.D+00 *U(N-5) 5 -50400.D+00 *U(N-4) 6 +18144.D+00 *U(N-3) 7 +20160.D+00 *U(N-2) 8 -2880.D+00 *U(N-1) 9 +240.D+00 *U(N )) C... C... GRID POINT N-2 UX(N-2)=R8FDX* 1 (+240.D+00 *U(N-8) 2 -2304.D+00 *U(N-7) 3 +10080.D+00 *U(N-6) 4 -26880.D+00 *U(N-5) 5 +50400.D+00 *U(N-4) 6 -80640.D+00 *U(N-3) 7 +38304.D+00 *U(N-2) 8 +11520.D+00 *U(N-1) 9 -720.D+00 *U(N )) C... C... GRID POINT N-1 UX(N-1)=R8FDX* 1 (-720.D+00 *U(N-8) 2 +6720.D+00 *U(N-7) 3 -28224.D+00 *U(N-6) 4 +70560.D+00 *U(N-5) 5 -117600.D+00 *U(N-4) 6 +141120.D+00 *U(N-3) 7 -141120.D+00 *U(N-2) 8 +64224.D+00 *U(N-1) 9 +5040.D+00 *U(N )) C... C... GRID POINT N UX(N )=R8FDX* 1 (+5040.D+00 *U(N-8) 2 -46080.D+00 *U(N-7) 3 +188160.D+00 *U(N-6) 4 -451584.D+00 *U(N-5) 5 +705600.D+00 *U(N-4) 6 -752640.D+00 *U(N-3) 7 +564480.D+00 *U(N-2) 8 -322560.D+00 *U(N-1) 9 +109584.D+00 *U(N )) RETURN END *DECK DSS010 SUBROUTINE DSS010(XL,XU,N,U,UX) C... C... SUBROUTINE DSS010 COMPUTES THE FIRST DERIVATIVE, U , OF A C... X C... VARIABLE U OVER THE SPATIAL DOMAIN XL LE X LE XU FROM CLASSICAL C... ELEVEN-POINT, TENTH-ORDER FINITE DIFFERENCE APPROXIMATIONS C... C... ARGUMENT LIST C... C... XL LOWER BOUNDARY VALUE OF X (INPUT) C... C... XU UPPER BOUNDARY VALUE OF X (INPUT) C... C... N NUMBER OF GRID POINTS IN THE X DOMAIN INCLUDING THE C... BOUNDARY POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF U AT C... THE N GRID POINT POINTS FOR WHICH THE DERIVATIVE IS C... TO BE COMPUTED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE NUMERICAL C... VALUES OF THE DERIVATIVES OF U AT THE N GRID POINTS C... (OUTPUT) C... C... THE MATHEMATICAL DETAILS OF THE FINITE DIFFERENCE APPROXIMA- C... TIONS ARE GIVEN IN SUBROUTINES DSS002, DSS004, DSS006 AND C... DSS008, AND CAN BE SUMMARIZED BY THE FOLLOWING BICKLEY MATRIX C... FOR N = 10, M = 1, P = 0 TO 10 (BICKLEY, W. G., FORMULAE FOR C... NUMERICAL DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941) C... C... -10626840 36288000 -81648000 C... 145152000 -190512000 182891520 C... C... -362880 -6636960 16329600 C... -21772800 25401600 -22861440 C... C... 40320 -806400 -4419360 C... 9676800 -8467200 6773760 C... C... -10080 151200 -1360800 C... -2756160 6350400 -3810240 C... C... 4320 -57600 388800 C... -2073600 -1330560 4354560 C... C... -2880 36000 -216000 C... 864000 -3024000 0 C... C... 2880 -34560 194400 C... -691200 1814400 -4354560 C... C... -4320 50400 -272160 C... 907200 -2116800 3810240 C... C... 10080 -115200 604800 C... -1935360 4233600 -6773760 C... C... -40320 453600 -2332800 C... 7257600 -15240960 22861440 C... C... 362880 -4032000 20412000 C... -62208000 127008000 -182891520 C... C... EACH ENTRY IN THIS TABLE SHOULD BE MULTIPLIED BY 1/10F TO C... OBTAIN THE FINAL WEIGHTING COEFFICIENTS. FROM THIS BICKLEY C... MATRIX, THE FINITE DIFFERENCE APPROXIMATION OF THE FIRST C... DERIVATIVE CAN BE PROGRAMMED FOR EACH OF THE GRID POINTS 1, 2, C... 3, 4, 5,..., I,..., N-4, N-3, N-2, N-1 AND N (TAKING INTO C... ACCOUNT THE SYMMETRY PROPERTIES OF THE MATRIX). C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R10FDX, U, UX, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT DX=(XU-XL)/DFLOAT(N-1) R10FDX=1.D+00/(3628800.D+00*DX) NM5=N-5 C... C... GRID POINT 1 UX( 1)=R10FDX* 1 (-10628640.D+00 *U( 1) 2 +36288000.D+00 *U( 2) 3 -81648000.D+00 *U( 3) 4 +145152000.D+00 *U( 4) 5 -190512000.D+00 *U( 5) 6 +182891520.D+00 *U( 6) 7 -127008000.D+00 *U( 7) 8 +62208000.D+00 *U( 8) 9 -20412000.D+00 *U( 9) A +4032000.D+00 *U( 10) B -362880.D+00 *U( 11)) C... C... GRID POINT 2 UX( 2)=R10FDX* 1 (-362880.D+00 *U( 1) 2 -6636960.D+00 *U( 2) 3 +16329600.D+00 *U( 3) 4 -21772800.D+00 *U( 4) 5 +25401600.D+00 *U( 5) 6 -22861440.D+00 *U( 6) 7 +15240960.D+00 *U( 7) 8 -7257600.D+00 *U( 8) 9 +2332800.D+00 *U( 9) A -453600.D+00 *U( 10) B +40320.D+00 *U( 11)) C... C... GRID POINT 3 UX( 3)=R10FDX* 1 (+40320.D+00 *U( 1) 2 -806400.D+00 *U( 2) 3 -4419360.D+00 *U( 3) 4 +9676800.D+00 *U( 4) 5 -8467200.D+00 *U( 5) 6 +6773760.D+00 *U( 6) 7 -4233600.D+00 *U( 7) 8 +1935360.D+00 *U( 8) 9 -604800.D+00 *U( 9) A +115200.D+00 *U( 10) B -10080.D+00 *U( 11)) C... C... GRID POINT 4 UX( 4)=R10FDX* 1 (-10080.D+00 *U( 1) 2 +151200.D+00 *U( 2) 3 -1360800.D+00 *U( 3) 4 -2756160.D+00 *U( 4) 5 +6350400.D+00 *U( 5) 6 -3810240.D+00 *U( 6) 7 +2116800.D+00 *U( 7) 8 -907200.D+00 *U( 8) 9 +272160.D+00 *U( 9) A -50400.D+00 *U( 10) B +4320.D+00 *U( 11)) C... C... GRID POINT 5 UX( 5)=R10FDX* 1 (+4320.D+00 *U( 1) 2 -57600.D+00 *U( 2) 3 +388800.D+00 *U( 3) 4 -2073600.D+00 *U( 4) 5 -1330560.D+00 *U( 5) 6 +4354560.D+00 *U( 6) 7 -1814400.D+00 *U( 7) 8 +691200.D+00 *U( 8) 9 -194400.D+00 *U( 9) A +34560.D+00 *U( 10) B -2880.D+00 *U( 11)) C... C... GRID POINT I, I NE 1 TO 5, N-4 TO N DO 1 I=6,NM5 UX( I)=R10FDX* 1 (-2880.D+00 *U( I-5) 2 +36000.D+00 *U( I-4) 3 -216000.D+00 *U( I-3) 4 +864000.D+00 *U( I-2) 5 -3024000.D+00 *U( I-1) 6 +0.D+00 *U( I ) 7 +3024000.D+00 *U( I+1) 8 -864000.D+00 *U( I+2) 9 +216000.D+00 *U( I+3) A -36000.D+00 *U( I+4) B +2880.D+00 *U( I+5)) 1 CONTINUE C... C... GRID POINT N-4 UX(N-4)=R10FDX* 1 (+2880.D+00 *U(N-10) 2 -34560.D+00 *U( N-9) 3 +194400.D+00 *U( N-8) 4 -691200.D+00 *U( N-7) 5 +1814400.D+00 *U( N-6) 6 -4354560.D+00 *U( N-5) 7 +1330560.D+00 *U( N-4) 8 +2073600.D+00 *U( N-3) 9 -388800.D+00 *U( N-2) A +57600.D+00 *U( N-1) B -4320.D+00 *U( N )) C... C... GRID POINT N-3 UX(N-3)=R10FDX* 1 (-4320.D+00 *U(N-10) 2 +50400.D+00 *U( N-9) 3 -272160.D+00 *U( N-8) 4 +907200.D+00 *U( N-7) 5 -2116800.D+00 *U( N-6) 6 +3810240.D+00 *U( N-5) 7 -6350400.D+00 *U( N-4) 8 +2756160.D+00 *U( N-3) 9 +1360800.D+00 *U( N-2) A -151200.D+00 *U( N-1) B +10080.D+00 *U( N )) C... C... GRID POINT N-2 UX(N-2)=R10FDX* 1 (+10080.D+00 *U(N-10) 2 -115200.D+00 *U( N-9) 3 +604800.D+00 *U( N-8) 4 -1935360.D+00 *U( N-7) 5 +4233600.D+00 *U( N-6) 6 -6773760.D+00 *U( N-5) 7 +8467200.D+00 *U( N-4) 8 -9676800.D+00 *U( N-3) 9 +4419360.D+00 *U( N-2) A +806400.D+00 *U( N-1) B -40320.D+00 *U( N )) C... C... GRID POINT N-1 UX(N-1)=R10FDX* 1 (-40320.D+00 *U(N-10) 2 +453600.D+00 *U( N-9) 3 -2332800.D+00 *U( N-8) 4 +7257600.D+00 *U( N-7) 5 -15240960.D+00 *U( N-6) 6 +22861440.D+00 *U( N-5) 7 -25401600.D+00 *U( N-4) 8 +21772800.D+00 *U( N-3) 9 -16329600.D+00 *U( N-2) A +6636960.D+00 *U( N-1) B +362880.D+00 *U( N )) C... C... GRID POINT N UX(N )=R10FDX* 1 (+362880.D+00 *U(N-10) 2 -4032000.D+00 *U( N-9) 3 +20412000.D+00 *U( N-8) 4 -62208000.D+00 *U( N-7) 5 +127008000.D+00 *U( N-6) 6 -182891520.D+00 *U( N-5) 7 +190512000.D+00 *U( N-4) 8 -145152000.D+00 *U( N-3) 9 +81648000.D+00 *U( N-2) A -36288000.D+00 *U( N-1) B +10628640.D+00 *U( N )) RETURN END *DECK DSS012 SUBROUTINE DSS012(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS012 IS AN APPLICATION OF FIRST-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS WITH THE C... SIMPLEST FORM C... C... U + V*U = 0 (1) C... T X C... C... THE FIRST FIVE PARAMETERS, XL, XU, N, U AND UX, ARE THE SAME C... AS FOR SUBROUTINES DSS002 TO DSS010 AS DEFINED IN THOSE ROUTINES. C... THE SIXTH PARAMETER, V, MUST BE PROVIDED TO DSS012 SO THAT THE C... DIRECTION OF FLOW IN EQUATION (1) CAN BE USED TO SELECT THE C... APPROPRIATE FINITE DIFFERENCE APPROXIMATION FOR THE FIRST-ORDER C... SPATIAL DERIVATIVE IN EQUATION (1), U . THE CONVENTION FOR THE C... SIGN OF V IS X C... C... FLOW LEFT TO RIGHT V GT 0 C... (I.E., IN THE DIRECTION (I.E., THE SIXTH ARGUMENT IS C... OF INCREASING X) POSITIVE IN CALLING DSS012) C... C... FLOW RIGHT TO LEFT V LT 0 C... (I.E., IN THE DIRECTION (I.E., THE SIXTH ARGUMENT IS C... OF DECREASING X) NEGATIVE IN CALLING DSS012) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, U, UX, V, XL, XU, 1 DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFERENCE C... APPROXIMATION DEPENDING ON THE SIGN OF V IN EQUATION (1). THE C... ORIGIN OF THE FINITE DIFFERENCE APPROXIMATIONS USED BELOW IS GIVEN C... AT THE END OF SUBROUTINE DSS012. DX=(XU-XL)/DFLOAT(N-1) IF(V.LT.0.D+00)GO TO 10 C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V UX(1)=(U(2)-U(1))/DX DO 1 I=2,N UX(I)=(U(I)-U(I-1))/DX 1 CONTINUE RETURN C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V 10 NM1=N-1 DO 2 I=1,NM1 UX(I)=(U(I+1)-U(I))/DX 2 CONTINUE UX(N)=(U(N)-U(N-1))/DX RETURN C... C... THE BACKWARD DIFFERENCES IN SECTION (1) ABOVE ARE BASED ON THE C... TAYLOR SERIES C... C... 2 3 C... UI-1 = UI + UI (-DX) + UI (-DX) + UI (-DX) + ... C... X 1F 2X 2F 3X 3F C... C... 2 C... IF THIS SERIES IS TRUNCATED AFTER THE DX TERM AND THE RESULTING C... EQUATION SOLVED FOR U , WE OBTAIN IMMEDIATELY C... X C... C... UI = (UI - UI-1)/DX + O(DX) C... X C... C... WHICH IS THE FIRST-ORDER BACKWARD DIFFERENCE USED IN DO LOOP 1. C... THE DERIVATIVE U1 IS COMPUTED BY USING THE POINT TO THE RIGHT OF C... X C... U1, I.E., U2, SINCE THIS IS THE ONLY POINT AVAILABLE IF FICTITIOUS C... POINTS TO THE LEFT OF U1 ARE TO BE AVOIDED. C... C... THE FORWARD DIFFERENCES IN SECTION (2) ABOVE ARE BASED ON THE C... TAYLOR SERIES C... C... 2 3 C... UI+1 = UI + UI ( DX) + UI ( DX) + UI ( DX) + ... C... X 1F 2X 2F 3X 3F C... C... 2 C... IF THIS SERIES IS TRUNCATED AFTER THE DX TERM AND THE RESULTING C... EQUATION SOLVED FOR U , WE OBTAIN IMMEDIATELY C... X C... C... UI = (UI+1 - UI)/DX + O(DX) C... X C... C... WHICH IS THE FIRST-ORDER FORWARD DIFFERENCE USED IN DO LOOP 2. C... THE DERIVATIVE UN IS COMPUTED BY USING THE POINT TO THE LEFT OF C... X C... UN (UN-1), SINCE THIS IS THE ONLY POINT AVAILABLE IF FICTITIOUS C... POINTS TO THE RIGHT OF UN ARE TO BE AVOIDED. END *DECK DSS014 SUBROUTINE DSS014(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS014 IS AN APPLICATION OF SECOND-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS- C... CUSSED IN SUBROUTINE DSS012. THE COEFFICIENTS OF THE FINITE C... DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W. C... G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL C... GAZETTE, PP. 19-27, 1941, N = 2, M = 1, P = 0, 1, 2. C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R2FDX, U, UX, V, XL, 1 XU, DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION C... CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER- C... ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT). DX=(XU-XL)/DFLOAT(N-1) R2FDX=1.D+00/(2.D+00*DX) IF(V.LT.0.D+00)GO TO 10 C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V UX(1)=R2FDX* 1( -3.D+00 *U( 1) 2 +4.D+00 *U( 2) 3 -1.D+00 *U( 3)) UX(2)=R2FDX* 1( -1.D+00 *U( 1) 2 +0.D+00 *U( 2) 3 +1.D+00 *U( 3)) DO 1 I=3,N UX(I)=R2FDX* 1( +1.D+00 *U(I-2) 2 -4.D+00 *U(I-1) 3 +3.D+00 *U(I )) 1 CONTINUE RETURN C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V 10 NM2=N-2 DO 2 I=1,NM2 UX(I)=R2FDX* 1( -3.D+00 *U(I ) 2 +4.D+00 *U(I+1) 3 -1.D+00 *U(I+2)) 2 CONTINUE UX(N-1)=R2FDX* 1( -1.D+00 *U(N-2) 2 +0.D+00 *U(N-1) 3 +1.D+00 *U(N )) UX(N)=R2FDX* 1( +1.D+00 *U(N-2) 2 -4.D+00 *U(N-1) 3 +3.D+00 *U(N )) RETURN END *DECK DSS016 SUBROUTINE DSS016(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS016 IS AN APPLICATION OF THIRD-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS- C... CUSSED IN SUBROUTINE DSS012. THE COEFFICIENTS OF THE FINITE C... DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W. C... G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL C... GAZETTE, PP. 19-27, 1941, N = 3, M = 1, P = 0, 1, 2, 3. C... DOUBLE PRECISION DX, R3FDX, U, UX, V, XL, 1 XU, DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION C... CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER- C... ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT). DX=(XU-XL)/DFLOAT(N-1) R3FDX=1.D+00/(6.D+00*DX) IF(V.LT.0.D+00)GO TO 10 C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V UX( 1)=R3FDX* 1( -11.D+00 *U( 1) 2 +18.D+00 *U( 2) 3 -9.D+00 *U( 3) 4 +2.D+00 *U( 4)) UX( 2)=R3FDX* 1( -2.D+00 *U( 1) 2 -3.D+00 *U( 2) 3 +6.D+00 *U( 3) 4 -1.D+00 *U( 4)) UX( 3)=R3FDX* 1( +1.D+00 *U( 1) 2 -6.D+00 *U( 2) 3 +3.D+00 *U( 3) 4 +2.D+00 *U( 4)) DO 1 I=4,N UX( I)=R3FDX* 1( -2.D+00 *U(I-3) 2 +9.D+00 *U(I-2) 3 -18.D+00 *U(I-1) 4 +11.D+00 *U(I )) 1 CONTINUE RETURN C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V 10 NM3=N-3 DO 2 I=1,NM3 UX( I)=R3FDX* 1( -11.D+00 *U(I ) 2 +18.D+00 *U(I+1) 3 -9.D+00 *U(I+2) 4 +2.D+00 *U(I+3)) 2 CONTINUE UX(N-2)=R3FDX* 1( -2.D+00 *U(N-3) 2 -3.D+00 *U(N-2) 3 +6.D+00 *U(N-1) 4 -1.D+00 *U(N )) UX(N-1)=R3FDX* 1( +1.D+00 *U(N-3) 2 -6.D+00 *U(N-2) 3 +3.D+00 *U(N-1) 4 +2.D+00 *U(N )) UX( N)=R3FDX* 1( -2.D+00 *U(N-3) 2 +9.D+00 *U(N-2) 3 -18.D+00 *U(N-1) 4 +11.D+00 *U(N )) RETURN END *DECK DSS018 SUBROUTINE DSS018(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS018 IS AN APPLICATION OF THIRD-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS- C... CUSSED IN SUBROUTINE DSS012. THE COEFFICIENTS OF THE FINITE C... DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W. C... G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL C... GAZETTE, PP. 19-27, 1941, N = 3, M = 1, P = 0, 1, 2, 3. THE C... IMPLEMENTATION IS THE **FOUR-POINT BIASED UPWIND FORMULA** OF C... M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE C... ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69, C... AUGUST, 1978 C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R3FDX, U, UX, V, XL, 1 XU, DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION C... CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER- C... ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT). DX=(XU-XL)/DFLOAT(N-1) R3FDX=1.D+00/(6.D+00*DX) IF(V.LT.0.D+00)GO TO 10 C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V UX( 1)=R3FDX* 1( -11.D+00 *U( 1) 2 +18.D+00 *U( 2) 3 -9.D+00 *U( 3) 4 +2.D+00 *U( 4)) UX( 2)=R3FDX* 1( -2.D+00 *U( 1) 2 -3.D+00 *U( 2) 3 +6.D+00 *U( 3) 4 -1.D+00 *U( 4)) NM1=N-1 DO 1 I=3,NM1 UX( I)=R3FDX* 1( +1.D+00 *U(I-2) 2 -6.D+00 *U(I-1) 3 +3.D+00 *U(I ) 4 +2.D+00 *U(I+1)) 1 CONTINUE UX( N)=R3FDX* 1( -2.D+00 *U(N-3) 2 +9.D+00 *U(N-2) 3 -18.D+00 *U(N-1) 4 +11.D+00 *U(N )) RETURN C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V 10 UX( 1)=R3FDX* 1( -11.D+00 *U( 1) 2 +18.D+00 *U( 2) 3 -9.D+00 *U( 3) 4 +2.D+00 *U( 4)) NM2=N-2 DO 2 I=2,NM2 UX( I)=R3FDX* 1( -2.D+00 *U(I-1) 2 -3.D+00 *U(I ) 3 +6.D+00 *U(I+1) 4 -1.D+00 *U(I+2)) 2 CONTINUE UX(N-1)=R3FDX* 1( +1.D+00 *U(N-3) 2 -6.D+00 *U(N-2) 3 +3.D+00 *U(N-1) 4 +2.D+00 *U(N )) UX( N)=R3FDX* 1( -2.D+00 *U(N-3) 2 +9.D+00 *U(N-2) 3 -18.D+00 *U(N-1) 4 +11.D+00 *U(N )) RETURN END *DECK DSS020 SUBROUTINE DSS020(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS020 IS AN APPLICATION OF FOURTH-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS- C... CUSSED IN SUBROUTINE DSS012. THE COEFFICIENTS OF THE FINITE C... DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W. C... G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL C... GAZETTE, PP. 19-27, 1941, N = 4, M = 1, P = 0, 1, 2, 3, 4. THE C... IMPLEMENTATION IS THE **FIVE-POINT BIASED UPWIND FORMULA** OF C... M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE C... ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69, C... AUGUST, 1978 C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION DX, R4FDX, U, UX, V, XL, 1 XU, DFLOAT C... DIMENSION U(N),UX(N) C... C... COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION C... CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER- C... ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT). DX=(XU-XL)/DFLOAT(N-1) R4FDX=1.D+00/(12.D+00*DX) IF(V.LT.0.D+00)GO TO 10 C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V UX( 1)=R4FDX* 1( -25.D+00 *U( 1) 2 +48.D+00 *U( 2) 3 -36.D+00 *U( 3) 4 +16.D+00 *U( 4) 5 -3.D+00 *U( 5)) UX( 2)=R4FDX* 1( -3.D+00 *U( 1) 2 -10.D+00 *U( 2) 3 +18.D+00 *U( 3) 4 -6.D+00 *U( 4) 5 +1.D+00 *U( 5)) UX( 3)=R4FDX* 1( +1.D+00 *U( 1) 2 -8.D+00 *U( 2) 3 +0.D+00 *U( 3) 4 +8.D+00 *U( 4) 5 -1.D+00 *U( 5)) NM1=N-1 DO 1 I=4,NM1 UX( I)=R4FDX* 1( -1.D+00 *U(I-3) 2 +6.D+00 *U(I-2) 3 -18.D+00 *U(I-1) 4 +10.D+00 *U(I ) 5 +3.D+00 *U(I+1)) 1 CONTINUE UX( N)=R4FDX* 1( +3.D+00 *U(N-4) 2 -16.D+00 *U(N-3) 3 +36.D+00 *U(N-2) 4 -48.D+00 *U(N-1) 5 +25.D+00 *U(N )) RETURN C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V 10 UX( 1)=R4FDX* 1( -25.D+00 *U( 1) 2 +48.D+00 *U( 2) 3 -36.D+00 *U( 3) 4 +16.D+00 *U( 4) 5 -3.D+00 *U( 5)) NM3=N-3 DO 2 I=2,NM3 UX( I)=R4FDX* 1( -3.D+00 *U(I-1) 2 -10.D+00 *U(I ) 3 +18.D+00 *U(I+1) 4 -6.D+00 *U(I+2) 5 +1.D+00 *U(I+3)) 2 CONTINUE UX(N-2)=R4FDX* 1( +1.D+00 *U(N-4) 2 -8.D+00 *U(N-3) 3 +0.D+00 *U(N-2) 4 +8.D+00 *U(N-1) 5 -1.D+00 *U(N )) UX(N-1)=R4FDX* 1( -1.D+00 *U(N-4) 2 +6.D+00 *U(N-3) 3 -18.D+00 *U(N-2) 4 +10.D+00 *U(N-1) 5 +3.D+00 *U(N )) UX( N)=R4FDX* 1( +3.D+00 *U(N-4) 2 -16.D+00 *U(N-3) 3 +36.D+00 *U(N-2) 4 -48.D+00 *U(N-1) 5 +25.D+00 *U(N )) RETURN END *DECK DSS024 SUBROUTINE DSS024(XL,XU,N,U,UX,V) C... C... SUBROUTINE DSS024 IS AN APPLICATION OF FOURTH-ORDER DIRECTIONAL C... DIFFERENCING IN THE NUMERICAL METHOD OF LINES. IT IS INTENDED C... SPECIFICALLY FOR THE ANALYSIS OF CONVECTIVE SYSTEMS MODELLED BY C... FIRST-ORDER HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS AS DIS- C... CUSSED IN SUBROUTINE DSS012. THE COEFFICIENTS OF THE FINITE C... DIFFERENCE APPROXIMATIONS USED HEREIN ARE TAKEN FROM BICKLEY, W. C... G., FORMULAE FOR NUMERICAL DIFFERENTIATION, THE MATHEMATICAL C... GAZETTE, PP. 19-27, 1941, N = 4, M = 1, P = 0, 1, 2, 3, 4. THE C... IMPLEMENTATION IS THE **FIVE-POINT BIASED UPWIND FORMULA** OF C... M. B. CARVER AND H. W. HINDS, THE METHOD OF LINES AND THE C... ADVECTION EQUATION, SIMULATION, VOL. 31, NO. 2, PP. 59-69, C... AUGUST, 1978. C... C... DSS024 IS AN EXTENSION OF DSS020. THESE ROUTINES HAVE THE C... FOLLOWING APPROXIMATIONS: C... C... DSS020 - FIVE POINT APPROXIMATIONS OVER THE ENTIRE SPATIAL C... GRID. AT THE INTERIOR POINTS, THE FIVE POINT, BIASED C... UPWIND APPROXMTAION IS USED. C... C... DSS024 - A COMBINATION OF TWO POINT, UPWIND AND FIVE POINT, C... BIASED UPWIND APPROXIMATIONS IS USED. THE FIVE POINT C... APPROXIMATION IS USED AT THE INTERIOR POINTS TO C... PROVIDE RESOLUTION OF MOVING FRONTS (I.E., TO MINIMIZE C... NUMERICAL DIFFUSION). THE TWO POINT APPROXIMATION C... IS USED NEAR THE BOUNDARIES TO REDUCE NUMERICAL C... OSCILLATIONS FROM THE FIVE POINT APPROXIMATION. C... C... TYPE REAL VARIABLES AS DOUBLE PRECISION IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N),UX(N) C... C... NUMBER OF GRID POINTS AT THE LEFT AND RIGHT ENDS WITH TWO POINT C... UPWIND APPROXIMATIONS NL=5 NU=5 C... C... COMPUTE THE COMMON FACTOR FOR EACH FINITE DIFFERENCE APPROXIMATION C... CONTAINING THE SPATIAL INCREMENT, THEN SELECT THE FINITE DIFFER- C... ENCE APPROXIMATION DEPENDING ON THE SIGN OF V (SIXTH ARGUMENT). DX=(XU-XL)/DFLOAT(N-1) R4FDX=1.D+00/(12.D+00*DX) IF(V.GE.0.D+00)THEN C... C... (1) FINITE DIFFERENCE APPROXIMATION FOR POSITIVE V C... C... TWO POINT UPWIND APPROXIMATIONS AT LEFT END (FOR NL GRID POINTS) UX( 1)=(1.0D0/DX)* 1( -1.D+00 *U(1) 2 +1.D+00 *U(2)) DO 2 I=2,NL UX( I)=(1.0D0/DX)* 1( -1.D+00 *U(I-1) 2 +1.D+00 *U(I )) 2 CONTINUE C... C... FIVE POINT BIASED UPWIND APPROXIMATIONS FOR INTERIOR GRID POINTS C... NL+1 TO N-NU-1 DO 1 I=NL+1,N-NU-1 UX( I)=R4FDX* 1( -1.D+00 *U(I-3) 2 +6.D+00 *U(I-2) 3 -18.D+00 *U(I-1) 4 +10.D+00 *U(I ) 5 +3.D+00 *U(I+1)) 1 CONTINUE C... C... TWO POINT UPWIND APPROXIMATIONS AT RIGHT END (FOR NU GRID POINTS) DO 3 I=N-NU,N UX(I )=(1.0D0/DX)* 1( -1.D+00 *U(I-1) 2 +1.D+00 *U(I )) 3 CONTINUE C... C... (2) FINITE DIFFERENCE APPROXIMATION FOR NEGATIVE V C... ELSE IF(V.LT.0.D+00)THEN C... C... TWO POINT UPWIND APPROXIMATIONS AT LEFT END (FOR NL GRID POINTS) DO 12 I=1,NL UX( I)=(1.0D0/DX)* 1( -1.D+00 *U(I ) 2 +1.D+00 *U(I+1)) 12 CONTINUE C... C... FIVE POINT BIASED UPWIND APPROXIMATIONS FOR INTERIOR GRID POINTS C... NL+1 TO N-NU-1 DO 11 I=NL+1,N-NU-1 UX( I)=R4FDX* 1( -3.D+00 *U(I-1) 2 -10.D+00 *U(I ) 3 +18.D+00 *U(I+1) 4 -6.D+00 *U(I+2) 5 +1.D+00 *U(I+3)) 11 CONTINUE C... C... TWO POINT UPWIND APPROXIMATIONS AT RIGHT END (FOR NU GRID POINTS) DO 13 I=N-NU,N-1 UX( I)=(1.0D0/DX)* 1( -1.D+00 *U(I ) 2 +1.D+00 *U(I+1)) 13 CONTINUE UX( N)=(1.0D0/DX)* 1( -1.D+00 *U(N-1) 2 +1.D+00 *U(N )) END IF RETURN END *DECK DSS026 SUBROUTINE DSS026(XL,XU,N1,N2,ND,U2D,UX2D,V) C... C... SUBROUTINE DSS026 COMPUTES A PARTIAL DERIVATIVE OVER A TWO- C... DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE- C... POINT BIASED UPWIND APPROXIMATIONS. IT IS INTENDED PRIMARILY C... FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION C... OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN TWO DIMENSIONS. IT C... DIFFERS FROM DSS034 ONLY IN THAT DSS024 IS CALLED IN PLACE OF C... DSS020. IT THEREFORE SUPERSEDES DSS034 (IN THE SENSE THAT DSS024 C... SUPERSEDES DSS020). C... C... ARGUMENT LIST C... C... XL LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... XU UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... N1 NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT C... VARIABLE (INPUT) C... C... N2 NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT C... VARIABLE (INPUT) C... C... ND NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE C... PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... U2D TWO-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT VARI- C... ABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT TO C... INDEPENDENT VARIABLE ND (INPUT) C... C... UX2D TWO-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI- C... VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO C... INDEPENDENT VARIABLE ND (OUTPUT) C... C... V VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED C... OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE C... PARTIAL DERIVATIVE. V EQ 0 CALLS THE FIVE-POINT C... CENTERED APPROXIMATION. V NE 0 CALLS THE FIVE-POINT C... BIASED UPWIND APPROXIMATION (INPUT) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION UX1D, UX2D, U1D, U2D, V, XL, 1 XU C... C... THE FOLLOWING TWO-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U2D) AND ITS PARTIAL DERIVATIVE (UX2D) DIMENSION U2D(N1,N2), UX2D(N1,N2) C... C... THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D). IN EACH C... CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE C... OTHER INDEPENDENT VARIABLE VARIES OVER ITS TOTAL INTERVAL. C... THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE C... ONE-DIMENSIONAL ROUTINES DSS004 AND DSS024. C... C... NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE- C... FORE HAVE TO BE INCREASED IN SIZE. HOWEVER, WITH A SIZE C... OF 101, THE TWO-DIMENSIONAL PROBLEM COULD HAVE A GRID OF C... 101 X 101 POINTS, THEREBY GENERATING AN APPROXIMATING ODE C... SYSTEM WITH A MULTIPLE OF 101 X 101 EQUATIONS, DEPENDING ON C... THE NUMBER OF SIMULTANEOUS PDES. THIS IS A VERY LARGE ODE C... PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING C... IS CONSIDERED ADEQUATE FOR MOST PROBLEMS. PARAMETER (NSIZE=101) DIMENSION U1D(NSIZE), UX1D(NSIZE) C... C... GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED C... WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE IF(ND.EQ.2)GO TO 2 C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N1 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 GRID POINTS VIA NESTED DO LOOPS 10, 11 AND 12 DO 10 J=1,N2 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS024 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 11 I=1,N1 U1D(I)=U2D(I,J) 11 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS024(XL,XU,N1,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE TWO-DIMENSIONAL ARRAY UX2D DO 12 I=1,N1 UX2D(I,J)=UX1D(I) 12 CONTINUE C... C... THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE SECOND INDE- C... PENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE CALCULATION FOR C... THE NEXT VALUE OF THE SECOND INDEPENDENT VARIABLE 10 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL C... DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D RETURN C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N2 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 GRID POINTS VIA NESTED DO LOOPS 20 AND 21 2 DO 20 I=1,N1 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS024 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 21 J=1,N2 U1D(J)=U2D(I,J) 21 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS024(XL,XU,N2,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE TWO-DIMENSIONAL ARRAY UX2D DO 22 J=1,N2 UX2D(I,J)=UX1D(J) 22 CONTINUE C... C... THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE FIRST INDE- C... PENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE CALCULATION FOR C... THE NEXT VALUE OF THE FIRST INDEPENDENT VARIABLE 20 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL C... DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D RETURN END *DECK DSS032 SUBROUTINE DSS032(NCALL,NPTS,IGRID,N,X,U,UX,C) C... C... SUBROUTINE DSS032 COMPUTES A NUMERICAL DERIVATIVE BY FIVE-POINT C... CENTERED AND NONCENTERED APPROXIMATIONS BASED ON THE FOURTH-ORDER C... LAGRANGE INTERPOLATON POLYNOMIAL. THE MATHEMATICAL DETAILS OF C... THE APPROXIMATIONS ARE GIVEN IN THE SUBSEQUENT COMMENTS, AND IN C... THE COMMENTS IN SUBROUTINE DSS32A, THE CORE DIFFERENTIATOR CALLED C... BY DSSO32. C... C... ARGUMENT LIST C... C... NCALL INTEGER INDEX FOR SPECIFYING THE CALCULATION OF THE C... WEIGHTING COEFFICIENTS IN THE DIFFERENTIATION FORMULAS C... DURING THE FIRST CALL TO SUBROUTINE DSS032. SPECIFI- C... CALLY, C... C... NCALL = 1 - THE FIRST CALL TO SUBROUTINE DSS032 C... FOR WHICH THE WEIGHTING COEFFICIENTS C... ARE CALCULATED C... C... NCALL NE 1 - A CALL TO SUBROUTINE DSS032 AFTER C... THE FIRST CALL SO THAT THE WEIGHTING C... COEFFICIENTS ARE NOT RECALCULATED C... C... (INPUT) C... C... NPTS NUMBER OF POINTS USED IN THE NUMERICAL DIFFERENTIATION C... FORMULAS. PRESENTLY, ONLY FIVE-POINT FORMULAS ARE USED C... IN SUBROUTINE DSS032 SO NPTS = 5 (INPUT) C... C... IGRID ONE-DIMENSIONAL ARRAY CONTAINING THE IDENTIFIERS FOR C... THE FIVE BASIC GRIDS DESCRIBED ABOVE (INPUT) C... C... N NUMBER OF GRID POINTS OVER WHICH THE INDEPENDENT C... VARIABLE, X, AND DEPENDENT VARIABLE, U, ARE DEFINED, C... INCLUDING THE END POINTS (INPUT) C... C... X ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE C... INDEPENDENT VARIABLE (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE C... DEPENDENT VARIABLE (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE N VALUES OF THE C... NUMERICAL DERIVATIVE OF U WITH RESPECT TO X (OUTPUT) C... C... C TWO-DIMENSIONAL ARRAY CONTAINING THE WEIGHTING COEFFI- C... CIENTS IN THE DIFFERENTIATION FORMULAS (OUTPUT) C... C... SUBROUTINE DSS032 IS A NUMERICAL DIFFERENTIATOR BASED ON THE FIVE- C... POINT LAGRANGIAN INTERPOLATION POLYNOMIAL. A NUMERICAL DERIVATIVE C... UX IS COMPUTED BASED ON ONE OF THE FOLLOWING FIVE-POINT GRIDS C... C... (1) UX C... I I+1 I+2 I+3 I+4 C... C... (2) UX C... I-1 I I+1 I+2 I+3 C... C... (3) UX C... I-2 I-1 I I+1 I+2 C... C... (4) UX C... I-3 I-2 I-1 I I+1 C... C... (5) UX C... I-4 I-3 I-2 I-1 I C... C... FOR GRID (1), FOUR POINTS ARE USED TO THE RIGHT OF THE POINT WHERE C... THE NUMERICAL DERIVATIVE, UX, IS COMPUTED. FOR GRID (2), ONE C... POINT TO THE LEFT AND THREE POINTS TO THE RIGHT OF THE POINT WHERE C... THE NUMERICAL DERIVATIVE IS COMPUTED ARE USED, ETC. THUS, GRID C... (3) IS A FIVE-POINT CENTERED APPROXIMATION, GRIDS (2) AND (4) ARE C... FIVE-POINT BIASED UPWIND APPROXIMATIONS, AND GRIDS (1) AND (5) ARE C... FIVE-POINT UPWIND APPROXIMATIONS. C... C... THE USER MAY SELECT ANY COMBINATION OF THESE GRIDS OVER A DOMAIN C... XL LE X LE XU CONSISTING OF A TOTAL OF N GRID POINTS. NORMALLY, C... GRIDS (1) AND (5) ARE USED AT THE LEFT AND RIGHT BOUNDARIES, X = C... XL AND X = XU, RESPECTIVELY. GRIDS (2) AND (4) ARE USED AT ONE C... GRID POINT TO THE RIGHT AND LEFT OF THE LEFT AND RIGHT BOUNDARIES, C... RESPECTIVELY. C... C... GRID (3) CAN BE USED IN THE INTERIOR OF THE X DOMAIN IN PROBLEMS C... FOR WHICH A CENTERED APPROXIMATION IS APPROPRIATE, E.G., A PARA- C... BOLIC PROBLEM IN PARTIAL DIFFERENTIAL EQUATIONS (PDES). C... C... GRID (4) CAN BE USED THROUGHOUT THE INTERIOR OF THE X DOMAIN FOR C... FIRST-ORDER HYPERBOLIC PDES WITH A POSITIVE VELOCITY, V, IN THE C... DIFFERENTIAL GROUP UT + V*UX, I.E., BIASED UPWIND APPROXIMATION OF C... THE SPATIAL DERIVATIVE, UX, WITH GREATER WEIGHTING TO THE LEFT C... (UPWIND) OF THE POINT WHERE UX IS CALCULATED. CONVERSELY, GRID C... (2) CAN BE USED THROUGHOUT THE INTERIOR OF THE X DOMAIN FOR FIRST- C... ORDER HYPERBOLIC PDES WITH A NEGATIVE VELOCITY, I.E., BIASED UP- C... WIND APPROXIMATION WITH GREATER WEIGHTING TO THE RIGHT (UPWIND) OF C... THE POINT WHERE UX IS CALCULATED. C... C... THE SELECTION OF THE FIVE BASIC GRIDS IS MADE BY THE USER THROUGH C... GRID IDENTIFIER ARRAY IGRID, THE THIRD ARGUMENT OF SUBROUTINE C... DSS032. FOR EXAMPLE, IF IGRID(1) = 1, IGRID(2) = 2, IGRID(I) = 3, C... IGRID(N-1) = 4, IGRID(N) = 5, THE FIRST AND SECOND GRID POINTS C... WILL RETURN THE DERIVATIVE UX ACCORDING TO GRIDS (1) AND (2) C... ABOVE. THE (N-1)ST AND NTH GRID POINTS WILL RETURN UX ACCORDING C... TO GRIDS (4) AND (5) ABOVE. THE INTERIOR GRID POINTS, WITH SUB- C... SCRIPT I, 3 LE I LE N-2, WILL RETURN UX ACCORDING TO GRID (3) C... ABOVE. C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION C, C5, U, UX, X, X5 C... C... VARIABLE DIMENSION THE ARRAYS FOR THE ARGUMENTS LISTED ABOVE DIMENSION IGRID(N), X(N), U(N), UX(N), C(NPTS,N) C... C... ABSOLUTE DIMENSION THE ARRAYS USED BY SUBROUTINE DSS32A WHICH IS C... CALLED BY DSS032 DIMENSION X5(5), C5(5) C... C... IF THIS IS THE FIRST CALL TO SUBROUTINE DSS032, GO TO STATEMENT C... 17 TO CALCULATE THE COEFFICIENTS IN THE DIFFERENTIATION FORMULAS C... VIA SUBROUTINE DSS32A IF(NCALL.EQ.1)GO TO 17 C... C... THE NUMERICAL DERIVATIVES ARE CALCULATED ON THE N GRID POINTS C... USING THE WEIGHTING COEFFICIENTS CALCULATED PREVIOUSLY BY SUB- C... ROUTINE DSS32A 7 DO 6 I=1,N ITYPE=IGRID(I) GO TO(1,2,3,4,5),ITYPE C... C... GRID TYPE 1 1 UX(I)=C(1,I)*U(I) 1 +C(2,I)*U(I+1) 2 +C(3,I)*U(I+2) 3 +C(4,I)*U(I+3) 4 +C(5,I)*U(I+4) GO TO 6 C... C... GRID TYPE 2 2 UX(I)=C(1,I)*U(I-1) 1 +C(2,I)*U(I) 2 +C(3,I)*U(I+1) 3 +C(4,I)*U(I+2) 4 +C(5,I)*U(I+3) GO TO 6 C... C... GRID TYPE 3 3 UX(I)=C(1,I)*U(I-2) 1 +C(2,I)*U(I-1) 2 +C(3,I)*U(I) 3 +C(4,I)*U(I+1) 4 +C(5,I)*U(I+2) GO TO 6 C... C... GRID TYPE 4 4 UX(I)=C(1,I)*U(I-3) 1 +C(2,I)*U(I-2) 2 +C(3,I)*U(I-1) 3 +C(4,I)*U(I) 4 +C(5,I)*U(I+1) GO TO 6 C... C... GRID TYPE 5 5 UX(I)=C(1,I)*U(I-4) 1 +C(2,I)*U(I-3) 2 +C(3,I)*U(I-2) 3 +C(4,I)*U(I-1) 4 +C(5,I)*U(I) 6 CONTINUE RETURN C... C... DURING THE FIRST CALL TO SUBROUTINE DSS032, THE COEFFICIENTS IN C... THE DIFFERENTIATION FORMULAS, STORED IN ARRAY C, ARE CALCULATED C... BY A CALL TO SUBROUTINE DSS32A 17 DO 16 I=1,N C... C... DEFINE THE VALUES OF X USED BY SUBROUTINE DSS32A TO CALCULATE C... THE WEIGHTING COEFFICIENTS IN THE DIFFERENTIATION FORMULAS ITYPE=IGRID(I) GO TO(11,12,13,14,15),ITYPE C... C... GRID TYPE 1 11 X5(1)=X(I) X5(2)=X(I+1) X5(3)=X(I+2) X5(4)=X(I+3) X5(5)=X(I+4) GO TO 19 C... C... GRID TYPE 2 12 X5(1)=X(I-1) X5(2)=X(I) X5(3)=X(I+1) X5(4)=X(I+2) X5(5)=X(I+3) GO TO 19 C... C... GRID TYPE 3 13 X5(1)=X(I-2) X5(2)=X(I-1) X5(3)=X(I) X5(4)=X(I+1) X5(5)=X(I+2) GO TO 19 C... C... GRID TYPE 4 14 X5(1)=X(I-3) X5(2)=X(I-2) X5(3)=X(I-1) X5(4)=X(I) X5(5)=X(I+1) GO TO 19 C... C... GRID TYPE 5 15 X5(1)=X(I-4) X5(2)=X(I-3) X5(3)=X(I-2) X5(4)=X(I-1) X5(5)=X(I) C... C... CALL SUBROUTINE DSS32A TO CALCULATE THE WEIGHTING COEFFICIENTS IN C... DIFFERENTIATION FORMULA FOR GRID POINT I 19 CALL DSS32A(ITYPE,X5,C5) C... C... STORE THE FIVE WEIGHTING COEFFICIENTS IN ARRAY C5 IN THE TWO- C... DIMENSIONAL ARRAY C DO 18 IC=1,5 C(IC,I)=C5(IC) 18 CONTINUE 16 CONTINUE C... C... ALL OF THE WEIGHTING COEFFICIENTS FOR THE N GRID POINTS HAVE BEEN C... CALCULATED AND STORED IN ARRAY C SO RETURN TO DO LOOP 6 TO CALCU- C... LATE THE NUMERICAL DERIVATIVES GO TO 7 END SUBROUTINE DSS32A(I,X,C) C... C... SUBROUTINE DSS32A IS THE CORE DIFFERENTIATOR CALLED BY SUBROUTINE C... DSS032. IT BASICALLY COMPUTES THE WEIGHTING COEFFCIENTS IN THE C... FIVE-POINT CENTERED AND NONCENTERED DIFFERENTIATION FORMULAS FOR C... USE IN SUBROUTINE DSS032. C... C... ARGUMENT LIST C... C... I TYPE OF GRID AS DISCUSSED AT THE BEGINNING OF SUBROU- C... TINE DSS032 WITH POSSIBLE VALUES 1, 2, 3, 4, 5 (INPUT) C... C... X ONE-DIMENSIONAL ARRAY CONTAINING THE FIVE VALUES OF THE C... INDEPENDENT VARIABLE, X, USED IN THE CALCULATION OF THE C... NUMERICAL DERIVATIVE (INPUT) C... C... C ONE-DIMENSIONAL ARRAY CONTAINING THE FIVE COEFFICIENTS C... OF THE DIFFERENTIATION FORMULA FOR GRID TYPE I (OUTPUT) C... C... SUBROUTINE DSS32A COMPUTES THE WEIGHTING COEFFICIENTS IN THE C... DIFFERENTIATION FORMULAS USED IN SUBROUTINE DSS032. DSS032 C... CALLS DSS32A JUST ONCE TO CALCULATE THE WEIGHTING COEFFICIENTS C... AFTER THE VARIABLE GRID IS DEFINED SINCE THE WEIGHTING CO- C... EFFICIENTS ARE CONSTANTS FOR A GIVEN VARIABLE GRID. C... C... THE GENERAL FIVE-POINT LAGRANGE INTERPOLATION POLYNOMIAL SERVES C... AS THE MATHEMATICAL BASIS FOR THIS SUBROUTINE C... C... C... (X - X2)(X - X3)(X - X4)(X - X5) C... U(X) = .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... (X - X1)(X - X3)(X - X4)(X - X5) C... + .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... (X - X1)(X - X2)(X - X4)(X - X5) C... + .................................... U(X3) (1) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... (X - X1)(X - X2)(X - X3)(X - X5) C... + .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... (X - X1)(X - X2)(X - X3)(X - X4) C... + .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... DIFFERENTIATION OF EQUATION (1) GIVES THE FIVE-POINT DIFFERENTIA- C... TION FORMULAS USED IN SUBROUTINE DSS032 FOR THE CALCULATION OF C... NUMERICAL DERIVATIVES OVER A VARIABLE GRID C... C... C... UX(X) = ((X - X2)(X - X3)(X - X4) + C... C... (X - X2)(X - X3)(X - X5) + C... C... (X - X2)(X - X4)(X - X5) + C... C... (X - X3)(X - X4)(X - X5)) C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... + ((X - X1)(X - X3)(X - X4) + C... C... (X - X1)(X - X3)(X - X5) + C... C... (X - X1)(X - X4)(X - X5) + C... C... (X - X3)(X - X4)(X - X5)) C... .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... + ((X - X1)(X - X2)(X - X4) + C... C... (X - X1)(X - X2)(X - X5) + C... C... (X - X1)(X - X4)(X - X5) + C... C... (X - X2)(X - X4)(X - X5)) C... .................................... U(X3) (2) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... + ((X - X1)(X - X2)(X - X3) + C... C... (X - X1)(X - X2)(X - X5) + C... C... (X - X1)(X - X3)(X - X5) + C... C... (X - X2)(X - X3)(X - X5)) C... .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... + ((X - X1)(X - X2)(X - X3) + C... C... (X - X1)(X - X2)(X - X4) + C... C... (X - X1)(X - X3)(X - X4) + C... C... (X - X2)(X - X3)(X - X4)) C... .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... IF X = X1 IN EQUATION (2), THE FORMULA FOR THE FIRST DERIVATIVE C... AT X = X1, DENOTED AS UX(X1), RESULTS C... C... C... UX(X1) = ((X1 - X2)(X1 - X3)(X1 - X4) + C... C... (X1 - X2)(X1 - X3)(X1 - X5) + C... C... (X1 - X2)(X1 - X4)(X1 - X5) + C... C... (X1 - X3)(X1 - X4)(X1 - X5)) C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... ((X1 - X3)(X1 - X4)(X1 - X5)) C... + .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... ((X1 - X2)(X1 - X4)(X1 - X5)) C... + .................................... U(X3) (3) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... ((X1 - X2)(X1 - X3)(X1 - X5)) C... + .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... ((X1 - X2)(X1 - X3)(X1 - X4)) C... + .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... SIMILARLY, THE SUBSTITUTIONS X = X2, X = X3, X = X4, X = X5 IN C... EQUATION (2) GIVES RESPECTIVELY THE DIFFERENTIATION FORMULAS FOR C... UX(X2), UX(X3), UX(X4), UX(X5) C... C... C... UX(X2) = ((X2 - X3)(X2 - X4)(X2 - X5)) C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... + ((X2 - X1)(X2 - X3)(X2 - X4) + C... C... (X2 - X1)(X2 - X3)(X2 - X5) + C... C... (X2 - X1)(X2 - X4)(X2 - X5) + C... C... (X2 - X3)(X2 - X4)(X2 - X5)) C... .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... ((X2 - X1)(X2 - X4)(X2 - X5)) C... + .................................... U(X3) (4) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... ((X2 - X1)(X2 - X3)(X2 - X5)) C... + .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... ((X2 - X1)(X2 - X3)(X2 - X4)) C... + .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... C... UX(X3) = ((X3 - X2)(X3 - X4)(X3 - X5)) C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... ((X3 - X1)(X3 - X4)(X3 - X5)) C... + .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... + ((X3 - X1)(X3 - X2)(X3 - X4) + C... C... (X3 - X1)(X3 - X2)(X3 - X5) + C... C... (X3 - X1)(X3 - X4)(X3 - X5) + C... C... (X3 - X2)(X3 - X4)(X3 - X5)) C... .................................... U(X3) (5) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... ((X3 - X1)(X3 - X2)(X3 - X5)) C... + .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... ((X3 - X1)(X3 - X2)(X3 - X4)) C... + .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... C... UX(X4) = ((X4 - X2)(X4 - X3)(X4 - X5)) C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... ((X4 - X1)(X4 - X3)(X4 - X5)) C... + .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... ((X4 - X1)(X4 - X2)(X4 - X5)) C... + .................................... U(X3) (6) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... + ((X4 - X1)(X4 - X2)(X4 - X3) + C... C... (X4 - X1)(X4 - X2)(X4 - X5) + C... C... (X4 - X1)(X4 - X3)(X4 - X5) + C... C... (X4 - X2)(X4 - X3)(X4 - X5)) C... .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... ((X4 - X1)(X4 - X2)(X4 - X3)) C... + .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... C... UX(X5) = ((X5 - X2)(X5 - X3)(X5 - X4) + C... .................................... U(X1) C... (X1 - X2)(X1 - X3)(X1 - X4)(X1 - X5) C... C... C... ((X5 - X1)(X5 - X3)(X5 - X4)) C... + .................................... U(X2) C... (X2 - X1)(X2 - X3)(X2 - X4)(X2 - X5) C... C... C... ((X5 - X1)(X5 - X2)(X5 - X4)) C... + .................................... U(X3) (7) C... (X3 - X1)(X3 - X2)(X3 - X4)(X3 - X5) C... C... C... ((X5 - X1)(X5 - X2)(X5 - X3)) C... + .................................... U(X4) C... (X4 - X1)(X4 - X2)(X4 - X3)(X4 - X5) C... C... C... + ((X5 - X1)(X5 - X2)(X5 - X3) + C... C... (X5 - X1)(X5 - X2)(X5 - X4) + C... C... (X5 - X1)(X5 - X3)(X5 - X4) + C... C... (X5 - X2)(X5 - X3)(X5 - X4)) C... .................................... U(X5) C... (X5 - X1)(X5 - X2)(X5 - X3)(X5 - X4) C... C... C... EQUATIONS (3) TO (7) ARE THE FIVE-POINT DIFFERENTIATION FORMULAS C... FOR A VARIABLE GRID. C... C... NOTE THAT THE FACTORS WHICH MULTIPLY U(X1), U(X2), U(X3), U(X4) C... AND U(X5) IN EQUATIONS (3) TO (7) ARE JUST CONSTANTS ONCE THE C... VALUES OF X1, X2, X3, X4 AND X5 ARE DEFINED. THUS, EQUATIONS C... (3) TO (7) CAN BE WRITTEN IN THE GENERAL FORM C... C... UX(XI) = C1*U(X1) + C2*U(X2) + C3*U(X3) + C4*U(X4) + C5*U(X5) C... C... WHERE XI = X1, X2, X3, X4, X5 FOR EQUATIONS (3) TO (7) RESPECT- C... IVELY, AND C1, C2, C3, C4 AND C5 ARE THE WEIGHTING COEFFICIENTS C... CALCULATED IN DSS32A AND PASSED TO DSS032 AS THE THIRD ARGUMENT C... (ARRAY C). C... C... IF THE FOLLOWING SUBSTUTIONS FOR A UNIFORM GRID ARE MADE IN EQUA- C... TIONS (3) TO (7) C... C... C... DX = (X2 - X1) = (X3 - X2) = (X4 - X3) = (X5 - X4) C... C... 2DX = (X3 - X1) = (X4 - X2) = (X5 - X3) C... C... 3DX = (X4 - X1) = (X5 - X2) C... C... 4DX = (X5 - X1) C... C... -DX = (X1 - X2) = (X2 - X3) = (X3 - X4) = (X4 - X5) C... C... -2DX = (X1 - X3) = (X2 - X4) = (X3 - X5) C... C... -3DX = (X1 - X4) = (X2 - X5) C... C... -4DX = (X1 - X5) C... C... C... C... EQUATIONS (3) TO (7) REDUCE TO THE FOLLOWING FIVE-POINT DIFFER- C... TIATION FORMULAE FOR A UNIFORM GRID C... C... C... UX(X1) = ((-1)(-2)(-3) + C... C... (-1)(-2)(-4) + C... C... (-1)(-3)(-4) + C... C... (-2)(-3)(-4)) C... ................. U(X1)/DX C... (-1)(-2)(-3)(-4) C... C... C... ((-2)(-3)(-4)) C... + ................. U(X2)/DX C... ( 1)(-1)(-2)(-3) C... C... C... ((-1)(-3)(-4)) C... + ................. U(X3)/DX C... ( 2)( 1)(-1)(-2) C... C... C... ((-1)(-2)(-4)) C... + ................. U(X4)/DX C... ( 3)( 2)( 1)(-1) C... C... C... ((-1)(-2)(-3)) C... + ................. U(X5)/DX C... ( 4)( 3)( 2)( 1) C... C... C... OR NUMERICALLY EVALUATING THE COEFFICIENTS OF THE U(X1) TO U(X5) C... TERMS (WITH 4F = 1*2*3*4 = 4 FACTORIAL) C... C... C... UX(X1) = (-50/4F)U(X1)/DX C... C... + ( 96/4F)U(X2)/DX C... C... + (-72/4F)U(X3)/DX (8) C... C... + ( 32/4F)U(X4)/DX C... C... + ( -6/4F)U(X5)/DX C... C... C... SIMILARLY, FROM EQUATIONS (4) TO (7), THE FOLLOWING DIFFERENTIA- C... TION FORMULAS FOR A UNIFORM GRID RESULT C... C... C... UX(X2) = ((-1)(-2)(-3)) C... ................. U(X1)/DX C... (-1)(-2)(-3)(-4) C... C... C... + (( 1)(-1)(-2) + C... C... ( 1)(-1)(-3) + C... C... ( 1)(-2)(-3) + C... C... (-1)(-2)(-3)) C... ................. U(X2)/DX C... ( 1)(-1)(-2)(-3) C... C... C... (( 1)(-2)(-3)) C... + ................. U(X3)/DX C... ( 2)( 1)(-1)(-2) C... C... C... (( 1)(-1)(-3)) C... + ................. U(X4)/DX C... ( 3)( 2)( 1)(-1) C... C... C... (( 1)(-1)(-2)) C... + ................. U(X5)/DX C... ( 4)( 3)( 2)( 1) C... C... C... UX(X2) = ( -6/4F)U(X1)/DX C... C... + (-20/4F)U(X2)/DX C... C... + ( 36/4F)U(X3)/DX (9) C... C... + (-12/4F)U(X4)/DX C... C... + ( 2/4F)U(X5)/DX C... C... C... C... UX(X3) = (( 1)(-1)(-2)) C... + ................. U(X1)/DX C... (-1)(-2)(-3)(-4) C... C... C... (( 2)(-1)(-2)) C... + ................. U(X2)/DX C... ( 1)(-1)(-2)(-3) C... C... C... + (( 2)( 1)(-1) + C... C... ( 2)( 1)(-2) + C... C... ( 2)(-1)(-2) + C... C... ( 1)(-1)(-2)) C... ................. U(X3)/DX C... ( 2)( 1)(-1)(-2) C... C... C... (( 2)( 1)(-2)) C... + ................. U(X4)/DX C... ( 3)( 2)( 1)(-1) C... C... C... (( 2)( 1)(-1)) C... + ................. U(X5)/DX C... ( 4)( 3)( 2)( 1) C... C... C... UX(X3) = ( 2/4F)U(X1)/DX C... C... + (-16/4F)U(X2)/DX C... C... + ( 0/4F)U(X3)/DX (10) C... C... + ( 16/4F)U(X4)/DX C... C... + ( -2/4F)U(X5)/DX C... C... C... C... UX(X4) = (( 2)( 1)(-1)) C... ................. U(X1)/DX C... (-1)(-2)(-3)(-4) C... C... C... (( 3)( 1)(-1)) C... + ................. U(X2)/DX C... ( 1)(-1)(-2)(-3) C... C... C... (( 3)( 2)(-1)) C... + ................. U(X3)/DX C... ( 2)( 1)(-1)(-2) C... C... C... + (( 3)( 2)( 1) + C... C... ( 3)( 2)(-1) + C... C... ( 3)( 1)(-1) + C... C... ( 2)( 1)(-1)) C... ................. U(X4)/DX C... ( 3)( 2)( 1)(-1) C... C... C... (( 3)( 2)( 1)) C... + ................. U(X5)/DX C... ( 4)( 3)( 2)( 1) C... C... C... UX(X4) = ( -2/4F)U(X1)/DX C... C... + ( 12/4F)U(X2)/DX C... C... + (-36/4F)U(X3)/DX (11) C... C... + ( 20/4F)U(X4)/DX C... C... + ( 6/4F)U(X5)/DX C... C... C... C... UX(X5) = (( 3)( 2)( 1)) C... ................. U(X1)/DX C... (-1)(-2)(-3)(-4) C... C... C... (( 4)( 2)( 1)) C... + ................. U(X2)/DX C... ( 1)(-1)(-2)(-3) C... C... C... (( 4)( 3)( 1)) C... + ................. U(X3)/DX C... ( 2)( 1)(-1)(-2) C... C... C... (( 4)( 3)( 2)) C... + ................. U(X4)/DX C... ( 3)( 2)( 1)(-1) C... C... C... + (( 4)( 3)( 2) + C... C... ( 4)( 3)( 1) + C... C... ( 4)( 2)( 1) + C... C... ( 3)( 2)( 1)) C... ................. U(X5)/DX C... ( 4)( 3)( 2)( 1) C... C... C... UX(X5) = ( 6/4F)U(X1)/DX C... C... + (-32/4F)U(X2)/DX C... C... + ( 72/4F)U(X3)/DX (12) C... C... + (-96/4F)U(X4)/DX C... C... + ( 50/4F)U(X5)/DX C... C... C... EQUATIONS (8) TO (12) CAN BE SUMMARIZED IN MATRIX FORM AS C... C... C... .. .. .. .. .. .. C... .UX(X1). . -50 96 -72 32 -6 . .U(X1). C... . . . . . . C... .UX(X2). . -6 -20 36 -12 2 . .U(X2). C... . . . . . . C... .UX(X3). = (1/4FDX). 2 -16 0 16 -2 .*.U(X3). (13) C... . . . . . . C... .UX(X4). . -2 12 -36 20 6 . .U(X4). C... . . . . . . C... .UX(X5). . 6 -32 72 -96 50 . .U(X5). C... .. .. .. .. .. .. C... C... C... THE 5 X 5 COEFFICIENT MATRIX, WITH THE FACTOR (1/4F), IS A BICKLEY C... MATRIX FOR N = 4, M = 1, P = 0, 1, 2, 3, 4 (BICKLEY, W. G., FORMU- C... LAE FOR NUMERICAL DIFFERENTIATION, MATH. GAZ., VOL. 25, 1941) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION C, SUM1, SUM2, SUM3, SUM4, SUM5, 1 TERM, TERM1, TERM2, TERM3, TERM4, TERM5, 2 X C... C... ABSOLUTE DIMENSION THE ARRAYS USED BY SUBROUTINE DSS32A DIMENSION X(5), C(5) C... C... SELECT ONE OF THE DIFFERENTATION FORMULAS, EQUATIONS (3) TO (7) GO TO(1,2,3,4,5),I C... C... BEGIN THE CALCULATION OF THE WEIGHTING COEFFICIENTS FOR FIVE- C... POINT NUMERICAL DIFFERENTIATION OVER A VARIABLE GRID C... C... ****************************************************************** C... C... EQUATION (3) C... C... COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3), C... U(X4), U(X5) TERMS IN EQUATION (3) BY DO LOOP 18, STARTING WITH C... THE DENOMINATOR OF EACH COEFFICIENT 1 DO 18 I=1,5 TERM=1.D+00 DO 11 J=1,5 IF(J.EQ.I)GO TO 11 TERM=TERM/(X(I)-X(J)) 11 CONTINUE IF(I.EQ.1)GO TO 13 C... C... COMPLETE THE COEFFICIENTS OF THE U(X2), U(X3), U(X4), U(X5) TERMS C... IN EQUATION (3) BY DO LOOP 12, ENDING WITH THE NUMERATOR OF EACH C... COEFFICIENT DO 12 J=1,5 IF(J.EQ.1)GO TO 12 IF(J.EQ.I)GO TO 12 TERM=(X(1)-X(J))*TERM 12 CONTINUE C(I)=TERM GO TO 18 C... C... COMPLETE THE COEFFICIENT OF THE U(X1) TERM IN EQUATION (3) BY DO C... LOOP 15 13 SUM1=0.D+00 DO 15 JSUM=1,5 IF(JSUM.EQ.1)GO TO 15 TERM1=1.D+00 DO 14 J=1,5 IF(J.EQ.1)GO TO 14 IF(J.EQ.JSUM)GO TO 14 TERM1=TERM1*(X(1)-X(J)) 14 CONTINUE SUM1=SUM1+TERM1 15 CONTINUE C(I)=SUM1*TERM 18 CONTINUE RETURN C... C... ****************************************************************** C... C... EQUATION (4) C... C... COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3), C... U(X4), U(X5) TERMS IN EQUATION (4) BY DO LOOP 28, STARTING WITH C... THE DENOMINATOR OF EACH COEFFICIENT 2 DO 28 I=1,5 TERM=1.D+00 DO 21 J=1,5 IF(J.EQ.I)GO TO 21 TERM=TERM/(X(I)-X(J)) 21 CONTINUE IF(I.EQ.2)GO TO 23 C... C... COMPLETE THE COEFFICIENTS OF THE U(X1), U(X3), U(X4), U(X5) TERMS C... IN EQUATION (4) BY DO LOOP 22, ENDING WITH THE NUMERATOR OF EACH C... COEFFICIENT DO 22 J=1,5 IF(J.EQ.2)GO TO 22 IF(J.EQ.I)GO TO 22 TERM=(X(2)-X(J))*TERM 22 CONTINUE C(I)=TERM GO TO 28 C... C... COMPLETE THE COEFFICIENT OF THE U(X2) TERM IN EQUATION (4) BY DO C... LOOP 25 23 SUM2=0.D+00 DO 25 JSUM=1,5 IF(JSUM.EQ.2)GO TO 25 TERM2=1.D+00 DO 24 J=1,5 IF(J.EQ.2)GO TO 24 IF(J.EQ.JSUM)GO TO 24 TERM2=TERM2*(X(2)-X(J)) 24 CONTINUE SUM2=SUM2+TERM2 25 CONTINUE C(I)=SUM2*TERM 28 CONTINUE RETURN C... C... ****************************************************************** C... C... EQUATION (5) C... C... COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3), C... U(X4), U(X5) TERMS IN EQUATION (5) BY DO LOOP 38, STARTING WITH C... THE DENOMINATOR OF EACH COEFFICIENT 3 DO 38 I=1,5 TERM=1.D+00 DO 31 J=1,5 IF(J.EQ.I)GO TO 31 TERM=TERM/(X(I)-X(J)) 31 CONTINUE IF(I.EQ.3)GO TO 33 C... C... COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X4), U(X5) TERMS C... IN EQUATION (5) BY DO LOOP 32, ENDING WITH THE NUMERATOR OF EACH C... COEFFICIENT DO 32 J=1,5 IF(J.EQ.3)GO TO 32 IF(J.EQ.I)GO TO 32 TERM=(X(3)-X(J))*TERM 32 CONTINUE C(I)=TERM GO TO 38 C... C... COMPLETE THE COEFFICIENT OF THE U(X3) TERM IN EQUATION (5) BY DO C... LOOP 35 33 SUM3=0.D+00 DO 35 JSUM=1,5 IF(JSUM.EQ.3)GO TO 35 TERM3=1.D+00 DO 34 J=1,5 IF(J.EQ.3)GO TO 34 IF(J.EQ.JSUM)GO TO 34 TERM3=TERM3*(X(3)-X(J)) 34 CONTINUE SUM3=SUM3+TERM3 35 CONTINUE C(I)=SUM3*TERM 38 CONTINUE RETURN C... C... ****************************************************************** C... C... EQUATION (6) C... C... COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3), C... U(X4), U(X5) TERMS IN EQUATION (6) BY DO LOOP 48, STARTING WITH C... THE DENOMINATOR OF EACH COEFFICIENT 4 DO 48 I=1,5 TERM=1.D+00 DO 41 J=1,5 IF(J.EQ.I)GO TO 41 TERM=TERM/(X(I)-X(J)) 41 CONTINUE IF(I.EQ.4)GO TO 43 C... C... COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X3), U(X5) TERMS C... IN EQUATION (6) BY DO LOOP 42, ENDING WITH THE NUMERATOR OF EACH C... COEFFICIENT DO 42 J=1,5 IF(J.EQ.4)GO TO 42 IF(J.EQ.I)GO TO 42 TERM=(X(4)-X(J))*TERM 42 CONTINUE C(I)=TERM GO TO 48 C... C... COMPLETE THE COEFFICIENT OF THE U(X4) TERM IN EQUATION (6) BY DO C... LOOP 45 43 SUM4=0.D+00 DO 45 JSUM=1,5 IF(JSUM.EQ.4)GO TO 45 TERM4=1.D+00 DO 44 J=1,5 IF(J.EQ.4)GO TO 44 IF(J.EQ.JSUM)GO TO 44 TERM4=TERM4*(X(4)-X(J)) 44 CONTINUE SUM4=SUM4+TERM4 45 CONTINUE C(I)=SUM4*TERM 48 CONTINUE RETURN C... C... ****************************************************************** C... C... EQUATION (7) C... C... COMPUTE THE WEIGHTING COEFFICIENTS OF THE U(X1), U(X2), U(X3), C... U(X4), U(X5) TERMS IN EQUATION (7) BY DO LOOP 58, STARTING WITH C... THE DENOMINATOR OF EACH COEFFICIENT 5 DO 58 I=1,5 TERM=1.D+00 DO 51 J=1,5 IF(J.EQ.I)GO TO 51 TERM=TERM/(X(I)-X(J)) 51 CONTINUE IF(I.EQ.5)GO TO 53 C... C... COMPLETE THE COEFFICIENTS OF THE U(X1), U(X2), U(X3), U(X4) TERMS C... IN EQUATION (7) BY DO LOOP 52, ENDING ITH THE NUMERATOR OF EACH C... COEFFICIENT DO 52 J=1,5 IF(J.EQ.5)GO TO 52 IF(J.EQ.I)GO TO 52 TERM=(X(5)-X(J))*TERM 52 CONTINUE C(I)=TERM GO TO 58 C... C... COMPLETE THE COEFFICIENT OF THE U(X5) TERM IN EQUATION (7) BY DO C... LOOP 55 53 SUM5=0.D+00 DO 55 JSUM=1,5 IF(JSUM.EQ.5)GO TO 55 TERM5=1.D+00 DO 54 J=1,5 IF(J.EQ.5)GO TO 54 IF(J.EQ.JSUM)GO TO 54 TERM5=TERM5*(X(5)-X(J)) 54 CONTINUE SUM5=SUM5+TERM5 55 CONTINUE C(I)=SUM5*TERM 58 CONTINUE RETURN C... C... ****************************************************************** C... END *DECK DSS034 SUBROUTINE DSS034(XL,XU,N1,N2,ND,U2D,UX2D,V) C... C... SUBROUTINE DSS034 COMPUTES A PARTIAL DERIVATIVE OVER A TWO- C... DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE- C... POINT BIASED UPWIND APPROXIMATIONS. IT IS INTENDED PRIMARILY C... FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION C... OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN TWO DIMENSIONS. C... C... SUBROUTINE DSS034 IS USED IN ESSENTIALLY THE SAME WAY AS C... SUBROUTINE DSS024. IT IS MORE GENERAL THAN DSS024 SINCE IT C... USES CENTERED OR BIASED UPWIND APPROXIMATIONS AS SELECTED BY C... THE USER. ALSO, THE CODING IN DSS034 IS MUCH SIMPLER THAN IN C... DSS024 BECAUSE CALLS TO THE ONE-DIMENSIONAL ROUTINES DSS004 C... (FIVE-POINT CENTERED) AND DSS020 (FIVE-POINT BIASED UPWIND) C... ARE USED IN PLACE OF MUCH OF THE COMPLEX CODING OF DSS024. C... BECAUSE OF THIS RELATIVE SIMPLICITY, DSS034 SHOULD BE AT LEAST C... AS EFFICIENT AS DSS024. THEREFORE, DSS034 IS RECOMMENDED FOR C... USE IN PLACE OF DSS024, PARTICULARLY FOR CONVECTIVE-DIFFUSION C... PROBLEMS WHICH REQUIRE COMBINATIONS OF BIASED UPWIND AND CENTERED C... APPROXIMATIONS FOR FIRST AND SECOND-ORDER SPATIAL DERIVATIVES IN C... PDES. C... C... ARGUMENT LIST C... C... XL LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... XU UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... N1 NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT C... VARIABLE (INPUT) C... C... N2 NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT C... VARIABLE (INPUT) C... C... ND NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE C... PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... U2D TWO-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT VARI- C... ABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT TO C... INDEPENDENT VARIABLE ND (INPUT) C... C... UX2D TWO-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI- C... VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO C... INDEPENDENT VARIABLE ND (OUTPUT) C... C... V VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED C... OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE C... PARTIAL DERIVATIVE. V EQ 0 CALLS THE FIVE-POINT C... CENTERED APPROXIMATION. V NE 0 CALLS THE FIVE-POINT C... BIASED UPWIND APPROXIMATION (INPUT) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION UX1D, UX2D, U1D, U2D, V, XL, 1 XU C... C... THE FOLLOWING TWO-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U2D) AND ITS PARTIAL DERIVATIVE (UX2D) DIMENSION U2D(N1,N2), UX2D(N1,N2) C... C... THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D). IN EACH C... CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE C... OTHER INDEPENDENT VARIABLE VARIES OVER ITS TOTAL INTERVAL. C... THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE C... ONE-DIMENSIONAL ROUTINES DSS004 AND DSS020. C... C... NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE- C... FORE HAVE TO BE INCREASED IN SIZE. HOWEVER, WITH A SIZE C... OF 51, THE TWO-DIMENSIONAL PROBLEM COULD HAVE A GRID OF C... 51 X 51 POINTS, THEREBY GENERATING AN APPROXIMATING ODE C... SYSTEM WITH A MULTIPLE OF 51 X 51 EQUATIONS, DEPENDING ON C... THE NUMBER OF SIMULTANEOUS PDES. THIS IS A VERY LARGE ODE C... PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING C... IS CONSIDERED ADEQUATE FOR MOST PROBLEMS. DIMENSION U1D(51), UX1D(51) C... C... GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED C... WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE IF(ND.EQ.2)GO TO 2 C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N1 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 GRID POINTS VIA NESTED DO LOOPS 10, 11 AND 12 DO 10 J=1,N2 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 11 I=1,N1 U1D(I)=U2D(I,J) 11 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS020(XL,XU,N1,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE TWO-DIMENSIONAL ARRAY UX2D DO 12 I=1,N1 UX2D(I,J)=UX1D(I) 12 CONTINUE C... C... THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE SECOND INDE- C... PENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE CALCULATION FOR C... THE NEXT VALUE OF THE SECOND INDEPENDENT VARIABLE 10 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL C... DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D RETURN C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N2 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 GRID POINTS VIA NESTED DO LOOPS 20 AND 21 2 DO 20 I=1,N1 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE TWO-DIMENSIONAL ARRAY U2D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 21 J=1,N2 U1D(J)=U2D(I,J) 21 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS020(XL,XU,N2,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE TWO-DIMENSIONAL ARRAY UX2D DO 22 J=1,N2 UX2D(I,J)=UX1D(J) 22 CONTINUE C... C... THE PARTIAL DERIVATIVE AT A PARTICULAR VALUE OF THE FIRST INDE- C... PENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE CALCULATION FOR C... THE NEXT VALUE OF THE FIRST INDEPENDENT VARIABLE 20 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE PARTIAL C... DERIVATIVE IN THE TWO-DIMENSIONAL ARRAY UX2D RETURN END *DECK DSS036 SUBROUTINE DSS036(XL,XU,N1,N2,N3,ND,U3D,UX3D,V) C... C... SUBROUTINE DSS036 COMPUTES A PARTIAL DERIVATIVE OVER A THREE- C... DIMENSIONAL DOMAIN USING EITHER FIVE-POINT CENTERED OR FIVE- C... POINT BIASED UPWIND APPROXIMATIONS. IT IS INTENDED PRIMARILY C... FOR THE NUMERICAL METHOD OF LINES (NMOL) NUMERICAL INTEGRATION C... OF PARTIAL DIFFERENTIAL EQUATIONS (PDES) IN THREE DIMENSIONS. C... C... SUBROUTINE DSS036 IS CALLED IN ESSENTIALLY THE SAME WAY AS C... SUBROUTINE DSS034. THE ONLY DIFFERENCE IS AN ADDITIONAL ARGU- C... MENT, N3, TO DEFINE THE NUMBER OF GRID POINTS IN THE THIRD C... DIMENSION. THE COMMENTS IN DSS034 SHOULD THEREFORE BE USEFUL C... IN UNDERSTANDING THE OPERATION OF DSS036. IN PARTICULAR, C... DSS036 CALLS SUBROUTINES DSS004 AND DSS020 TO IMPLEMENT THE C... FIVE-POINT CENTERED APPROXIMATION AND FIVE-POINT BIASED UPWIND C... APPROXIMATION OF THE PARTIAL DERIAVTIVE, RESPECTIVELY. C... C... ARGUMENT LIST C... C... XL LOWER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... XU UPPER VALUE OF THE INDEPENDENT VARIABLE FOR WHICH C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... N1 NUMBER OF GRID POINTS FOR THE FIRST INDEPENDENT C... VARIABLE (INPUT) C... C... N2 NUMBER OF GRID POINTS FOR THE SECOND INDEPENDENT C... VARIABLE (INPUT) C... C... N3 NUMBER OF GRID POINTS FOR THE THIRD INDEPENDENT C... VARIABLE (INPUT) C... C... ND NUMBER OF THE INDEPENDENT VARIABLE FOR WHICH THE C... PARTIAL DERIVATIVE IS TO BE COMPUTED (INPUT) C... C... U3D THREE-DIMENSIONAL ARRAY CONTAINING THE DEPENDENT C... VARIABLE WHICH IS TO BE DIFFERENTIATED WITH RESPECT C... TO INDEPENDENT VARIABLE ND (INPUT) C... C... UX3D THREE-DIMENSIONAL ARRAY CONTAINING THE PARTIAL DERI- C... VATIVE OF THE DEPENDENT VARIABLE WITH RESPECT TO C... INDEPENDENT VARIABLE ND (OUTPUT) C... C... V VARIABLE TO SELECT EITHER THE FIVE-POINT CENTERED C... OR FIVE-POINT BIASED UPWIND APPROXIMATION FOR THE C... PARTIAL DERIVATIVE. V EQ 0 CALLS THE FIVE-POINT C... CENTERED APPROXIMATION. V NE 0 CALLS THE FIVE-POINT C... BIASED UPWIND APPROXIMATION (INPUT) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION UX1D, UX3D, U1D, U3D, V, XL, 1 XU C... C... THE FOLLOWING THREE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U3D) AND ITS PARTIAL DERIVATIVE (UX3D) DIMENSION U3D(N1,N2,N3), UX3D(N1,N2,N3) C... C... THE FOLLOWING ONE-DIMENSIONAL ARRAYS CONTAIN THE DEPENDENT C... VARIABLE (U1D) AND ITS PARTIAL DERIVATIVE (UX1D). IN EACH C... CASE, ONE OF THE INDEPENDENT VARIABLES IS CONSTANT AND THE C... OTHER TWO INDEPENDENT VARIABLES VARY OVER THEIR TOTAL INTERVALS. C... THESE ARRAYS ARE USED FOR TEMPORARY STORAGE IN CALLING THE C... ONE-DIMENSIONAL ROUTINES DSS004 AND DSS020. C... C... NOTE THAT THE ARRAYS HAVE ABSOLUTE DIMENSIONS AND MAY THERE- C... FORE HAVE TO BE INCREASED IN SIZE. HOWEVER, WITH A SIZE C... OF 51, THE THREE-DIMENSIONAL PROBLEM COULD HAVE A GRID OF C... 51 X 51 X 51 POINTS, THEREBY GENERATING AN APPROXIMATING ODE C... SYSTEM WITH A MULTIPLE OF 51 X 51 X 51 EQUATIONS, DEPENDING ON C... THE NUMBER OF SIMULTANEOUS PDES. THIS IS A VERY LARGE ODE C... PROBLEM, AND THEREFORE THE FOLLOWING ABSOLUTE DIMENSIONING C... IS CONSIDERED ADEQUATE FOR MOST PROBLEMS. DIMENSION U1D(51), UX1D(51) C... C... GO TO STATEMENT 2 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED C... WITH RESPECT TO THE SECOND INDEPENDENT VARIABLE IF(ND.EQ.2)GO TO 2 C... C... GO TO STATEMENT 3 IF THE PARTIAL DERIVATIVE IS TO BE COMPUTED C... WITH RESPECT TO THE THIRD INDEPENDENT VARIABLE IF(ND.EQ.3)GO TO 3 C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... FIRST INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N1 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 X N3 GRID POINTS VIA NESTED DO LOOPS 10, 11, 12 AND 13 DO 10 J=1,N2 DO 11 K=1,N3 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 12 I=1,N1 U1D(I)=U3D(I,J,K) 12 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N1,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS020(XL,XU,N1,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE THREE-DIMENSIONAL ARRAY UX3D DO 13 I=1,N1 UX3D(I,J,K)=UX1D(I) 13 CONTINUE C... C... THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE SECOND AND C... THIRD INDEPENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE C... CALCULATION FOR THE OTHER VALUES OF THE SECOND AND THIRD C... INDEPENDENT VARIABLES 11 CONTINUE 10 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 X N3 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE C... PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D RETURN C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... SECOND INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N2 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 X N3 GRID POINTS VIA NESTED DO LOOPS 20, 21, 22 AND 23 2 DO 20 I=1,N1 DO 21 K=1,N3 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 22 J=1,N2 U1D(J)=U3D(I,J,K) 22 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N2,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS020(XL,XU,N2,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE THREE-DIMENSIONAL ARRAY UX3D DO 23 J=1,N2 UX3D(I,J,K)=UX1D(J) 23 CONTINUE C... C... THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE FIRST AND C... THIRD INDEPENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE C... CALCULATION FOR THE OTHER VALUES OF THE FIRST AND THIRD C... INDEPENDENT VARIABLES 21 CONTINUE 20 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 X N3 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE C... PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D RETURN C... C... ****************************************************************** C... C... THE PARTIAL DERIVATIVE IS TO BE COMPUTED WITH RESPECT TO THE C... THIRD INDEPENDENT VARIABLE DEFINED OVER AN INTERVAL CONSISTING C... OF N3 GRID POINTS. COMPUTE THE PARTIAL DERIVATIVE AT THE N1 X C... N2 X N3 GRID POINTS VIA NESTED DO LOOPS 30, 31, 32 AND 33 3 DO 30 I=1,N1 DO 31 J=1,N2 C... C... TRANSFER THE DEPENDENT VARIABLE IN THE THREE-DIMENSIONAL ARRAY U3D C... TO THE ONE-DIMENSIONAL ARRAY U1D SO THAT SUBROUTINES DSS004 AND C... DSS020 CAN BE USED TO CALCULATE THE PARTIAL DERIVATIVE DO 32 K=1,N3 U1D(K)=U3D(I,J,K) 32 CONTINUE C... C... IF V EQ 0, A FIVE-POINT CENTERED APPROXIMATION IS USED FOR THE C... PARTIAL DERIVATIVE IF(V.EQ.0.D+00)CALL DSS004(XL,XU,N3,U1D,UX1D) C... C... IF V NE 0, A FIVE-POINT BIASED UPWIND APPROXIMATION IS USED FOR C... THE PARTIAL DERIVATIVE IF(V.NE.0.D+00)CALL DSS020(XL,XU,N3,U1D,UX1D,V) C... C... RETURN THE PARTIAL DERIVATIVE IN THE ONE-DIMENSIONAL ARRAY UX1D C... TO THE THREE-DIMENSIONAL ARRAY UX3D DO 33 K=1,N3 UX3D(I,J,K)=UX1D(K) 33 CONTINUE C... C... THE PARTIAL DERIVATIVE AT PARTICULAR VALUES OF THE FIRST AND C... SECOND INDEPENDENT VARIABLE HAS BEEN CALCULATED. REPEAT THE C... CALCULATION FOR THE OTHER VALUES OF THE FIRST AND SECOND C... INDEPENDENT VARIABLES 31 CONTINUE 30 CONTINUE C... C... THE PARTIAL DERIVATIVE HAS BEEN CALCULATED OVER THE ENTIRE N1 X C... N2 X N3 GRID. THEREFORE RETURN TO THE CALLING PROGRAM WITH THE C... PARTIAL DERIVATIVE IN THE THREE-DIMENSIONAL ARRAY UX3D RETURN END *DECK DSS038 SUBROUTINE DSS038(N,X,U,UX,B,C,D) C... C... SUBROUTINE DSS038 IS A NUMERICAL DIFFERENTAITION ROUTINE BASED ON C... THE CUBIC SPLINE OF FORSYTHE, G. E., M. A. MALCOLM, AND C. B. C... MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE- C... HALL, ENGLEWOOD CLIFFS, N.J., 1977, PP. 76-79. C... C... ARGUMENT LIST C... C... N NUMBER OF SPATIAL GRID POINTS (INPUT) C... C... X ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE C... SPATIAL (BOUNDARY-VALUE) INDEPENDENT VARIABLE (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE C... DEPENDENT VARIABLE (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE DERIVATIVE OF U C... WITH RESPECT TO X (OUTPUT) C... C... B,C,D ONE-DIMENSIONAL ARRAYS CONTAINING THE CUBIC SPLINE CO- C... EFFICIENTS (OUTPUT) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION B, C, D, U, UX, X C... C... VARIABLE DIMENSION THE ARRAYS PASSED AS ARGUMENTS DIMENSION X(N), U(N), UX(N), B(N), C(N), D(N) C... C... CALL SUBROUTINE SPLINE TO COMPUTE THE SPLINE COEFFICIENTS CALL SPLINE(N,X,U,B,C,D) C... C... THE NUMERICAL DERIVATIVE IS OBTAINED BY DIFFERENTIATION OF THE C... CUBIC SPLINE C... C... S(X) = U(I) + B(I)*(X - X(I)) + C(I)*(X - X(I))**2 C... C... + D(I)*(X - X(I))**3 C... THUS C... C... UX(I) = DS(X(I))/DX = B(I) C... DO 1 I=1,N UX(I)=B(I) 1 CONTINUE RETURN END *DECK DSS040 SUBROUTINE DSS040(N,X,U,UX,UXX,B,C,D) C... C... SUBROUTINE DSS040 IS A NUMERICAL DIFFERENTAITION ROUTINE BASED ON C... THE CUBIC SPLINE OF FORSYTHE, G. E., M. A. MALCOLM, AND C. B. C... MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE- C... HALL, ENGLEWOOD CLIFFS, N.J., 1977, PP. 76-79. C... C... ARGUMENT LIST C... C... N NUMBER OF SPATIAL GRID POINTS (INPUT) C... C... X ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE C... SPATIAL (BOUNDARY-VALUE) INDEPENDENT VARIABLE (INPUT) C... C... U ONE-DIMENSIONAL ARRAY CONTAINING THE VALUES OF THE C... DEPENDENT VARIABLE (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY CONTAINING THE FIRST DERIVATIVE C... OF U WITH RESPECT TO X (OUTPUT) C... C... UXX ONE-DIMENSIONAL ARRAY CONTAINING THE SECOND DERIVATIVE C... OF U WITH RESPECT TO X (OUTPUT) C... C... B,C,D ONE-DIMENSIONAL ARRAYS CONTAINING THE CUBIC SPLINE CO- C... EFFICIENTS (OUTPUT) C... C... TYPE SELECTED REAL VARIABLES AS DOUBLE PRECISION DOUBLE PRECISION B, C, D, U, UX, UXX, 1 X C... C... VARIABLE DIMENSION THE ARRAYS PASSED AS ARGUMENTS DIMENSION X(N), U(N), UX(N),UXX(N), B(N), C(N), D(N) C... C... CALL SUBROUTINE SPLINE TO COMPUTE THE SPLINE COEFFICIENTS CALL SPLINE(N,X,U,B,C,D) C... C... THE NUMERICAL DERIVATIVES, UX AND UXX, ARE COMPUTED BY DIFFERENT- C... IATION OF THE CUBIC SPLINE C... C... S(X) = U(I) + B(I)*(X - X(I)) + C(I)*(X - X(I))**2 C... C... + D(I)*(X - X(I))**3 C... THUS C... C... UX(I) = DS(X(I))/DX = B(I) C... C... UXX(I) = D(DS(X(I))/DX)/DX = 2*C(I) C... DO 1 I=1,N UX(I)=B(I) UXX(I)=2.D+00*C(I) 1 CONTINUE RETURN END *DECK SPLINE SUBROUTINE SPLINE (N, X, Y, B, C, D) C... C... THIS SUBROUTINE APPEARS IN C... C... FORSYTHE, G. E., M. A. MALCOLM AND C. B. MOLER, COMPUTER C... METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-HALL, INC., C... ENGLEWOOD CLIFFS, N.J., 1977, 76-78 C... INTEGER N DOUBLE PRECISION X(N), Y(N), B(N), C(N), D(N) C C THE COEFFICIENTS B(I), C(I), AND D(I), I=1,2,...,N ARE COMPUTED C FOR A CUBIC INTERPOLATING SPLINE C C S(X) = Y(I) + B(I)*(X-X(I)) + C(I)*(X-X(I))**2 + D(I)*(X-X(I))**3 C C FOR X(I) .LE. X .LE. X(I+1) C C INPUT.. C C N = THE NUMBER OF DATA POINTS OR KNOTS (N.GE.2) C X = THE ABSCISSAS OF THE KNOTS IN STRICTLY INCREASING ORDER C Y = THE ORDINATES OF THE KNOTS C C OUTPUT.. C C B, C, D = ARRAYS OF SPLINE COEFFICIENTS AS DEFINED ABOVE. C C USING P TO DENOTE DIFFERENTIATION, C C Y(I) = S(X(I)) C B(I) = SP(X(I)) C C(I) = SPP(X(I))/2 C D(I) = SPPP(X(I))/6 (DERIVATIVE FROM THE RIGHT) C C THE ACCOMPANYING FUNCTION SUBPROGRAM SEVAL CAN BE USED C TO EVALUATE THE SPLINE. C C INTEGER NM1, IB, I DOUBLE PRECISION T C NM1 = N-1 IF ( N .LT. 2 ) RETURN IF ( N .LT. 3 ) GO TO 50 C C SET UP TRIDIAGONAL SYSTEM C C B = DIAGONAL, D = OFFDIAGONAL, C = RIGHT HAND SIDE. C D(1) = X(2) - X(1) C(2) = (Y(2) - Y(1))/D(1) DO 10 I = 2, NM1 D(I) = X(I+1) - X(I) B(I) = 2.D+00*(D(I-1) + D(I)) C(I+1) = (Y(I+1) - Y(I))/D(I) C(I) = C(I+1) - C(I) 10 CONTINUE C C END CONDITIONS. THIRD DERIVATIVES AT X(1) AND X(N) C OBTAINED FROM DIVIDED DIFFERENCES C B(1) = -D(1) B(N) = -D(N-1) C(1) = 0.D+00 C(N) = 0.D+00 IF ( N .EQ. 3 ) GO TO 15 C(1) = C(3)/(X(4)-X(2)) - C(2)/(X(3)-X(1)) C(N) = C(N-1)/(X(N)-X(N-2)) - C(N-2)/(X(N-1)-X(N-3)) C(1) = C(1)*D(1)**2/(X(4)-X(1)) C(N) = -C(N)*D(N-1)**2/(X(N)-X(N-3)) C C FORWARD ELIMINATION C 15 DO 20 I = 2, N T = D(I-1)/B(I-1) B(I) = B(I) - T*D(I-1) C(I) = C(I) - T*C(I-1) 20 CONTINUE C C BACK SUBSTITUTION C C(N) = C(N)/B(N) DO 30 IB = 1, NM1 I = N-IB C(I) = (C(I) - D(I)*C(I+1))/B(I) 30 CONTINUE C C C(I) IS NOW THE SIGMA(I) OF THE TEXT C C COMPUTE POLYNOMIAL COEFFICIENTS C B(N) = (Y(N) - Y(NM1))/D(NM1) + D(NM1)*(C(NM1) + 2.D+00*C(N)) DO 40 I = 1, NM1 B(I) = (Y(I+1) - Y(I))/D(I) - D(I)*(C(I+1) + 2.D+00*C(I)) D(I) = (C(I+1) - C(I))/D(I) C(I) = 3.D+00*C(I) 40 CONTINUE C(N) = 3.D+00*C(N) D(N) = D(N-1) RETURN C 50 B(1) = (Y(2)-Y(1))/(X(2)-X(1)) C(1) = 0.D+00 D(1) = 0.D+00 B(2) = B(1) C(2) = 0.D+00 D(2) = 0.D+00 RETURN END *DECK DSS042 SUBROUTINE DSS042(XL,XU,N,U,UX,UXX,NL,NU) C... C... SUBROUTINE DSS042 COMPUTES A SECOND-ORDER APPROXIMATION OF A C... SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE C... AT THE BOUNDARY. C... C... DOUBLE PRECISION C... C... ARGUMENT LIST C... C... XL LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT) C... C... XU RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT) C... C... N NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END C... POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE C... DIFFERENTIATED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U. C... THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN C... NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU, C... DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE- C... SCRIPTION OF NL AND NU BELOW) C... C... UXX ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U C... (OUTPUT) C... C... NL INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT C... X = XL (INPUT). THE ALLOWABLE VALUES ARE C... C... 1 - DIRICHLET BOUNDARY CONDITION AT X = XL C... (UX(1) IS NOT USED) C... C... 2 - NEUMANN BOUNDARY CONDITION AT X = XL C... (UX(1) IS USED) C... C... NU INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT C... X = XU (INPUT). THE ALLOWABLE VALUES ARE C... C... 1 - DIRICHLET BOUNDARY CONDITION AT X = XU C... (UX(N) IS NOT USED) C... C... 2 - NEUMANN BOUNDARY CONDITION AT X = XU C... (UX(N) IS USED) C... C... TYPE REAL VARIABLES AS DOUBLE PRECISION IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... DIMENSION U(N), UX(N), UXX(N) C... C... THE FOLLOWING DERIVATION IS FOR A SET OF SECOND-ORDER, FOUR-POINT C... APPROXIMATIONS FOR A SECOND DERIVATIVE THAT CAN BE USED AT THE C... BOUNDARIES OF A SPATIAL DOMAIN. THESE APPROXIMATIONS HAVE THE C... FEATURES C... C... (1) ONLY INTERIOR AND BOUNDARY POINTS ARE USED (I.E., NO C... FICTITIOUS POINTS ARE USED) C... C... (2) THE NORMAL DERIVATIVE AT THE BOUNDARY IS INCLUDED AS PART C... OF THE APPROXIMATION FOR THE SECOND DERIVATIVE C... C... (3) APPROXIMATIONS FOR THE BOUNDARY CONDITIONS ARE NOT USED. C... C... THE DERIVATION IS BY PROFESSOR GILBERT A. STENGLE, DEPARTMENT OF C... MATHEMATICS, LEHIGH UNIVERSITY, BETHLEHEM, PA 18015, AND WAS DONE C... ON DECEMBER 7, 1985. C... C... FOR AN APPROXIMATION AT THE LEFT BOUNDARY, INVOLVING THE POINTS C... I, I+1, I+2 AND I+3, CONSIDER THE FOLLOWING TAYLOR SERIES EXPAN- C... SIONS C... C... UX(I)( DX) UXX(I)( DX)**2 UXXX(I)( DX)**3 C... U(I+1) = U(I) + ---------- + -------------- + --------------- +... C... 1 2 6 C... C... C... UX(I)(2DX) UXX(I)(2DX)**2 UXXX(I)(2DX)**3 C... U(I+2) = U(I) + ---------- + -------------- + --------------- +... C... 1 2 6 C... C... IF WE NOW FORM THE FOLLOWING LINEAR COMBINATION, INVOLVING CON- C... STANTS A, B, C AND D TO BE DETERMINED, AND USE THE PRECEDING TWO C... TAYLOR SERIES, C... C... A*U(I) + B*UX(I) + C*U(I+1) + D*U(I+2) C... C... WE HAVE C... C... A*U(I) + B*UX(I) + C*U(I+1) + D*U(I+2) = C... C... (A + B + C + D)*U(I) + C... C... (B + DX*C + 2*DX*D)*UX(I) + C... C... (C*(DX**2)/2 + D*((2*DX)**2)/2)*UXX(I) + C... C... (C*(DX**3)/6 + D*((2*DX)**3)/6)*UXXX(I) + O(DX**4) C... C... THE THIRD DERIVATIVE, UXXX(I), CAN BE DROPPED BY TAKING C... C... C = -8*D C... C... THE SECOND DERIVATIVE, UXX(I), CAN BE RETAINED BY TAKING C... C... (DX**2)(C/2 + 2*D) = 1 C... C... WHICH, WHEN COMBINED WITH THE PRECEDING RESULT GIVES C... C... D = -1/(2*(DX**2)) C... C... C = 4/(DX**2) C... C... THE FIRST DERIVATIVE, UX(I), CAN BE DROPPED BY TAKING C... C... B + DX*C + 2*DX*D = 0 C... C... OR C... C... B = -DX*C - 2*DX*D = -4/DX - 2*DX*(-1/(2*(DX**2))) = -3/DX C... C... FINALLY, U(I), CAN BE DROPPED BY TAKING C... C... A = - C - D = 8*D - D = -7*D = -7/(2*(DX**2)) C... C... IF WE NOW SOLVE FOR THE DERIVATIVE OF INTEREST, UXX(I), C... C... UXX(I) = -7/(2(DX**2))*U(I) - 3/DX*UX(I) C... C... + 8/(DX**2)*U(I+1) - 1/(2*(DX**2))U(I+2) + O(DX**2) C... C... = (1/(2*(DX**2)))*(-U(I+2) + 8*U(I+1) - 7*U(I) - 6*DX*UX(I)) C... C... + O(DX**2) C... C... WHICH IS THE FOUR-POINT, SECOND-ORDER APPROXIMATION FOR THE SECOND C... DERIVATIVE, UXX(I), INCLUDING THE FIRST DERIVATIVE, UX(I). C... C... FOUR CHECKS OF THIS APPROXIMATION CAN EASILY BE MADE FOR U(I) = C... 1, U(I) = X, U(I) = X**2 AND U(I) = X**3 C... C... UXX(I) = (1/(2*(DX**2)))*(-1 + 8*1 - 7*1 - 6*DX*0) = 0 C... C... UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX) + 8*(X + DX) C... C... -7*X - 6*DX*1) = 0 C... C... UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX)**2 + 8*(X + DX)**2 C... C... - 7*(X**2) - 6*DX*(2*X)) C... C... = (- X**2 - 4*X*DX - 4*DX**2 C... C... + 8*X**2 + 16*X*DX + 8*DX**2 C... C... - 7*X**2 - 12*X*DX)/(2*(DX**2)) = 2 C... C... UXX(I) = (1/(2*(DX**2)))*(-(X + 2*DX)**3 + 8*(X + DX)**3 C... C... - 7*(X**3) - 6*DX*(3*X**2)) C... C... = (1/(2*(DX**2)))*(- X**3 - 6*DX*X**2 - 12*X*DX**2 C... C... - 8*DX**3 + 8*X**3 + 24*DX*X**2 + 24*X*DX**2 + 8*DX**3 C... C... - 7*X**3 - 18*DX*X**2) C... C... = (1/(2*(DX**2)))*(12*X*DX**2) = 6*X C... C... THE PRECEDING APPROXIMATION FOR UXX(I) CAN BE APPLIED AT THE C... LEFT BOUNDARY VALUE OF X BY TAKING I = 1. AN APPROXIMATION AT C... THE RIGHT BOUNDARY IS OBTAINED BY TAKING DX = -DX AND REVERSING C... THE SUBSCRIPTS IN THE PRECEDING APPROXIMATION, WITH I = N C... C... UXX(I) C... C... = (1/(2*(DX**2)))*(-U(I-2) + 8*U(I-1) - 7*U(I) + 6*DX*UX(I)) C... C... + O(DX**2) C... C... TO OBTAIN APPROXIMATIONS OF THE SECOND DERVIAVTIVE WHICH DO NOT C... INVOLVE THE FIRST DERIVATIVE, WE TAKE AS THE LINEAR COMBINATION C... C... A*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) = C... C... WE HAVE C... C... A*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) = C... C... (A + B + C + D)*U(I)+ C... C... (DX*B + 2*DX*C + 4*DX*D)*UX(I)+ C... C... (B*(DX**2)/2 + C*((2*DX)**2)/2 + D*((3*DX)**2)/2)*UXX(I) + C... C... (B*(DX**3)/6 + C*((2*DX)**3)/6 + D*((3*DX)**3)/6)*UXX(I) + C... C... O(DX**4) C... C... THE THIRD DERIVATIVE, UXXX(I), CAN BE DROPPED BY TAKING C... C... B + 8*C + 27*D = 0 C... C... THE SECOND DERIVATIVE, UXX(I), CAN BE RETAINED BY TAKING C... C... (DX**2)*(B/2 + 2*C + (9/2)*D) = 1 C... C... THE FIRST DERIVATIVE CAN BE DROPPED BY TAKING C... C... B + 2*C + 3*D = 0 C... C... SOLUTION OF THE PRECEDING EQUATIONS FOR C AND D BY ELIMINATION OF C... B GIVES C... C... 6*C + 24*D = 0 C... C... 4*C + 18*D = -2/(DX**2) C... C... THEN, ELIMINATING C, GIVES C... C... (18 - 16)*D = -2/(DX**2) C... C... OR C... C... D = -1/(DX**2) C... C... C = (24/6)/(DX**2) = 4/(DX**2) C... C... B = -8/(DX**2) + 3/(DX**2) = -5/(DX**2) C... C... U(I) CAN BE DROPPED BY TAKING C... C... A + B + C + D = 0 C... C... OR C... C... A = (5 - 4 + 1)/(DX**2) = 2/(DX**2) C... C... IF WE NOW SOLVE FOR THE DERIVATIVE OF INTEREST, UXX(I), C... C... UXX(I) = (1/DX**2)*(2*U(I) - 5*U(I+1) + 4*U(I+2) - 1*U(I+3)) C... C... + O(DX**2) C... C... WHICH IS THE FOUR-POINT, SECOND-ORDER APPROXIMATION FOR THE SECOND C... DERIVATIVE, UXX(I), WITHOUT THE FIRST DERIVATIVE, UX(I). C... C... FOUR CHECKS OF THIS APPROXIMATION CAN EASILY BE MADE FOR U(I) = C... 1, U(I) = X, U(I) = X**2 AND U(I) = X**3 C... C... UXX(I) = (1/DX**2)*(2 - 5 + 4 - 1) = 0 C... C... UXX(I) = (1/DX**2)*(2*X - 5*(X + DX) + 4*(X + 2*DX) C... C... - 1*(X + 3*DX)) = 0 C... C... UXX(I) = (1/DX**2)*(2*X**2 - 5*(X + DX)**2 + 4*(X + 2*DX)**2 C... C... - 1*(X + 3*DX)**2) = 2 C... C... UXX(I) = (1/DX**2)*(2*X**3 - 5*(X + DX)**3 + 4*(X + 2*DX)**3 C... C... - 1*(X + 3*DX)**3) C... C... = (1/DX**2)*(2*X**3 - 5*X**3 - 15*X*DX**2 C... C... - 15*DX*X**2 - 5*DX**3 + 4*X**3 + 24*DX*X**2 C... C... + 48*X*DX**2 + 32*DX**3 - X**3 - 9*DX*X**2 C... C... - 27*X*DX**2 - 27DX**3) C... C... = (1/DX**2)*(6*X*DX**2) = 6*X C... C... THE PRECEDING APPROXIMATION FOR UXX(I) CAN BE APPLIED AT THE C... LEFT BOUNDARY VALUE OF X BY TAKING I = 1. AN APPROXIMATION AT C... THE RIGHT BOUNDARY IS OBTAINED BY TAKING DX = -DX AND REVERSING C... THE SUBSCRIPTS IN THE PRECEDING APPROXIMATION, WITH I = N C... C... UXX(I) = (1/DX**2)*(2*U(I) - 5*U(I-1) + 4*U(I-2) - 1*U(I-3)) C... C... + O(DX**2) C... C... GRID SPACING DX=(XU-XL)/DFLOAT(N-1) C... C... CALCULATE UXX AT THE LEFT BOUNDARY, WITHOUT UX IF(NL.EQ.1)THEN UXX(1)=(( 2.D0)*U( 1) 1 +( -5.D0)*U( 2) 2 +( 4.D0)*U( 3) 3 +( -1.D0)*U( 4))/(DX**2) C... C... CALCULATE UXX AT THE LEFT BOUNDARY, INCLUDING UX ELSE IF(NL.EQ.2)THEN UXX(1)=(( -7.D0)*U( 1) 1 +( 8.D0)*U( 2) 2 +( -1.D0)*U( 3))/(2.D0*DX**2) 3 +( -6.D0)*UX( 1) /(2.D0*DX) END IF C... C... CALCULATE UXX AT THE RIGHT BOUNDARY, WITHOUT UX IF(NU.EQ.1)THEN UXX(N)=(( 2.D0)*U(N ) 1 +( -5.D0)*U(N-1) 2 +( 4.D0)*U(N-2) 3 +( -1.D0)*U(N-3))/(DX**2) C... C... CALCULATE UXX AT THE RIGHT BOUNDARY, INCLUDING UX ELSE IF(NU.EQ.2)THEN UXX(N)=(( -7.D0)*U(N ) 1 +( 8.D0)*U(N-1) 2 +( -1.D0)*U(N-2))/(2.D0*DX**2) 3 +( 6.D0)*UX(N ) /(2.D0*DX) END IF C... C... CALCULATE UXX AT THE INTERIOR GRID POINTS DO 1 I=2,N-1 UXX(I)=(U(I+1)-2.D0*U(I)+U(I-1))/DX**2 1 CONTINUE RETURN END *DECK DSS044 SUBROUTINE DSS044(XL,XU,N,U,UX,UXX,NL,NU) C... C... SUBROUTINE DSS044 COMPUTES A FOURTH-ORDER APPROXIMATION OF A C... SECOND-ORDER DERIVATIVE, WITH OR WITHOUT THE NORMAL DERIVATIVE C... AT THE BOUNDARY. C... C... SINGLE PRECISION C... C... ARGUMENT LIST C... C... XL LEFT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT) C... C... XU RIGHT VALUE OF THE SPATIAL INDEPENDENT VARIABLE (INPUT) C... C... N NUMBER OF SPATIAL GRID POINTS, INCLUDING THE END C... POINTS (INPUT) C... C... U ONE-DIMENSIONAL ARRAY OF THE DEPENDENT VARIABLE TO BE C... DIFFERENTIATED (INPUT) C... C... UX ONE-DIMENSIONAL ARRAY OF THE FIRST DERIVATIVE OF U. C... THE END VALUES OF UX, UX(1) AND UX(N), ARE USED IN C... NEUMANN BOUNDARY CONDITIONS AT X = XL AND X = XU, C... DEPENDING ON THE ARGUMENTS NL AND NU (SEE THE DE- C... SCRIPTION OF NL AND NU BELOW) C... C... UXX ONE-DIMENSIONAL ARRAY OF THE SECOND DERIVATIVE OF U C... (OUTPUT) C... C... NL INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT C... X = XL (INPUT). THE ALLOWABLE VALUES ARE C... C... 1 - DIRICHLET BOUNDARY CONDITION AT X = XL C... (UX(1) IS NOT USED) C... C... 2 - NEUMANN BOUNDARY CONDITION AT X = XL C... (UX(1) IS USED) C... C... NU INTEGER INDEX FOR THE TYPE OF BOUNDARY CONDITION AT C... X = XU (INPUT). THE ALLOWABLE VALUES ARE C... C... 1 - DIRICHLET BOUNDARY CONDITION AT X = XU C... (UX(N) IS NOT USED) C... C... 2 - NEUMANN BOUNDARY CONDITION AT X = XU C... (UX(N) IS USED) C... C... TYPE REAL VARIABLES AS DOUBLE PRECISION IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... DIMENSION U(N), UX(N), UXX(N) C... C... THE FOLLOWING DERIVATION WAS COMPLETED BY W. E. SCHIESSER, DEPTS C... OF CHE AND CSEE, LEHIGH UNIVERSITY, BETHLEHEM, PA 18015, USA, ON C... DECEMBER 15, 1986. ADDITIONAL DETAILS ARE GIVEN IN SUBROUTINE C... DSS042. C... C... ****************************************************************** C... C... (1) UXX AT THE INTERIOR POINTS 3, 4,..., N-2 C... C... TO DEVELOP A SET OF FOURTH-ORDER CORRECT DIFFERENTIATION FORMULAS C... FOR THE SECOND DERIVATIVE UXX, WE CONSIDER FIRST THE INTERIOR C... GRID POINTS AT WHICH A SYMMETRIC FORMULA CAN BE USED. C... C... IF WE CONSIDER A FORMULA OF THE FORM C... C... A*U(I-2) + B*U(I-1) + E*U(I) + C*U(I+1) + D*U(I+2) C... C... TAYLOR SERIES EXPANSIONS OF U(I-2), U(I-1), U(I+1) AND U(I+2) C... CAN BE SUBSTITUTED INTO THIS FORMULA. WE THEN CONSIDER THE C... LINEAR ALBEGRAIC EQUATIONS RELATING A, B, C AND D WHICH WILL C... RETAIN CERTAIN TERMS, I.E., UXX, AND DROP OTHERS, E.G., UXXX, C... UXXXX AND UXXXXX. C... C... THUS, FOR GRID POINTS 3, 4,..., N-2 C... C... TO RETAIN UXX C... C... 4*A + B + C + 4*D = 2 (1) C... C... TO DROP UXXX C... C... -8*A - B + C + 8*D = 0 (2) C... C... TO DROP UXXXX C... C... 16*A + B + C + 16*D = 0 (3) C... C... TO DROP UXXXXX C... C... -32*A - B + C + 32*D = 0 (4) C... C... EQUATIONS (1) TO (4) CAN BE SOLVED FOR A, B, C AND D. IF EQUA- C... TION (1) IS ADDED TO EQUATION (2) C... C... -4*A + 2*C + 12*D = 2 (5) C... C... IF EQUATION (1) IS SUBTRACTED FROM EQUATION (3) C... C... 12*A + 12*D = -2 (6) C... C... IF EQUATION (1) IS ADDED TO EQUATION (4) C... C... -28*A + 2*C + 36*D = 2 (7) C... C... EQUATIONS (5) TO (7) CAN BE SOLVED FOR A, C AND D. IF EQUATION C... (5) IS SUBTRACTED FROM EQUATION (7), AND THE RESULT COMBINED C... WITH EQUATION (6) C... C... 12*A + 12*D = -2 (6) C... C... -24*A + 24*D = 0 (8) C... C... EQUATIONS (6) AND (8) CAN BE SOLVED FOR A AND D. FROM (8), A C... = D. FROM EQUATION (6), A = -1/12 AND D = -1/12. THEN, FROM C... EQUATION (5), C = 4/3, AND FROM EQUATION (1), B = 4/3. C... C... THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS C... C... (-1/12)*U(I-2) + (4/3)*U(I-1) + C... C... (4/3)*U(I+1) + (-1/12)*U(I+2) C... C... (-1/12 + 4/3 - 1/12 + 4/3)*U(I) + UXX(I)*(DX**2) + O(DX**6) C... C... OR C... C... UXX(I) = (1/(12*DX**2))*(-1*U(I-2) + 16*U(I-1) C... C... - 30*U(I) + 16*U(I+1) - 1*U(I+2) (9) C... C... + O(DX**4) C... C... NOTE THAT THE UX TERM DROPS OUT, I.E., THE BASIC EQUATION IS C... C... -2*A - B + C + 2*D = C... C... -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0 C... C... EQUATION (9) WAS OBTAINED BY DROPPING ALL TERMS IN THE UNDERLYING C... TAYLOR SERIES UP TO AND INCLUDING THE FIFTH DERIVATIVE, UXXXXX. C... THUS, EQUATION (9) IS EXACT FOR POLYNOMIALS UP TO AND INCLUDING C... FIFTH ORDER. THIS CAN BE CHECKED BY SUBSTITUTING THE FUNCTIONS C... 1, X, X**2, X**3, X**4 AND X**5 IN EQUATION (9) AND COMPUTING THE C... CORRESPONDING DERIVATIVES FOR COMPARISON WITH THE KNOWN SECOND C... DERIVATIVES. THIS IS DONE FOR 1 MERELY BY SUMMING THE WEIGHTING C... COEFFICIENTS IN EQUATION (9), WHICH SHOULD SUM TO ZERO, I.E., C... -1 + 16 - 30 + 16 -1 = 0. C... C... FOR THE REMAINING FUNCTIONS, THE ALGEBRA IS RATHER INVOLVED, BUT C... THESE FUNCTIONS CAN BE CHECKED NUMERICALLY, I.E., NUMERICAL VALUES C... OF X**2, X**3, X**4 AND X**5 CAN BE SUBSTITUTED IN EQUATION (9) C... AND THE COMPUTED DERIVATIVES CAN BE COMPARED WITH THE KNOW NUMERI- C... CAL SECOND DERIVATIVES. THIS IS NOT A PROOF OF CORRECTNESS OF C... EQUATION (9), BUT WOULD LIKELY DETECT ANY ERRORS IN EQUATION (9). C... C... CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (1) TO (4) C... CALL CHECK(1) C... C... ****************************************************************** C... C... (2) UXX AT THE INTERIOR POINTS I = 2 AND N-1 C... C... FOR GRID POINT 2, WE CONSIDER A FORMULA OF THE FORM C... C... A*U(I-1) + F*U(I) + B*U(I+1) + C*U(I+2) + D*U(I+3) + E*U(I+4) C... C... TAYLOR SERIES EXPANSIONS OF U(I-1), U(I+1), U(I+2), U(I+3) AND C... U(I+4) WHEN SUBSTITUTED INTO THIS FORMULA GIVE LINEAR ALGEBRAIC C... EQUATIONS RELATING A, B, C, D AND E. C... C... TO DROP UX C... C... -A + B + 2*C + 3*D + 4*E = 0 (10) C... C... TO RETAIN UXX C... C... A + B + 4*C + 9*D + 16*E = 2 (11) C... C... TO DROP UXXX C... C... -A + B + 8*C + 27*D + 64*E = 0 (12) C... C... TO DROP UXXXX C... C... A + B + 16*C + 81*D + 256*E = 0 (13) C... C... TO DROP UXXXXX C... C... -A + B + 32*C + 243*D + 1024*E = 0 (14) C... C... EQUATIONS (11), (12), (13) AND (14) CAN BE SOLVED FOR A, B, C, C... D AND E. IF EQUATION (10) IS ADDED TO EQUATION (11) C... C... 2*B + 6*C + 12*D +20*E = 2 (15) C... C... IF EQUATION (10) IS SUBTRACTED FROM EQUATION (12) C... C... 6*C + 24*D + 60*E = 0 (16) C... C... IF EQUATION (10) IS ADDED TO EQUATION (13) C... C... 2*B + 18*C + 84*D + 260*E = 0 (17) C... C... IF EQUATION (10) IS SUBTRACTED FROM EQUATION (14) C... C... 30*C + 240*D + 1020*E = 0 (18) C... C... EQUATIONS (15), (16), (17) AND (18) CAN BE SOLVED FOR B, C, D C... AND E. C... C... 6*C + 24*D + 60*E = 0 (16) C... C... IF EQUATION (15) IS SUBTRACTED FROM EQUATION (17) C... C... 12*C + 72*D + 240*E = -2 (19) C... C... 30*C + 240*D + 1020*E = 0 (18) C... C... EQUATIONS (16), (18) AND (19) CAN BE SOLVED FOR C, D AND E. IF C... TWO TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (19), C... C... 24*D + 120*E = -2 (20) C... C... IF FIVE TIMES EQUATION (16) IS SUBTRACTED FROM EQUATION (18), C... C... 120*D + 720*E = 0 (21) C... C... EQUATIONS (20) AND (21) CAN BE SOLVED FOR D AND E. FROM (21), C... E = (-1/6)*D. SUBSTITUTION IN EQUATION (20) GIVES D = -1/2. C... THUS, E = 1/12. FROM EQUATION (16), C = 7/6. FROM EQUATION C... (15), B = -1/3. FROM EQUATION (10), A = 5/6. C... C... THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS C... C... (5/6)*U(I-1) + (-1/3)*U(I+1) + (7/6)*U(I+2) + (-1/2)*U(I+3) C... C... + (1/12)*U(I+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*U(I) C... C... + UXX*(DX**2) + O(DX**6) C... C... OR C... C... UXX(I) = (1/12*DX**2)*(10*U(I-1) - 15*U(I) - 4*U(I+1) C... (22) C... + 14*U(I+2) - 6*U(I+3) + 1*U(I+4)) + O(DX**4) C... C... EQUATION (22) WILL BE APPLIED AT I = 2 AND N-1. THUS C... C... UXX(2) = (1/12*DX**2)*(10*U(1) - 15*U(2) - 4*U(3) C... (23) C... + 14*U(4) - 6*U(5) + 1*U(6)) + O(DX**4) C... C... UXX(N-1) = (1/12*DX**2)*(10*U(N) - 15*U(N-1) - 4*U(N-2) C... (24) C... + 14*U(N-3) - 6*U(N-4) + 1*U(N-5)) + O(DX**4) C... C... CHECK THE VALUES OF A, B, C, D AND E IN EQUATIONS (10) TO (14) C... CALL CHECK(2) C... C... ****************************************************************** C... C... (3) UXX AT THE BOUNDARY POINTS 1 AND N C... C... FINALLY, FOR GRID POINT 1, AN APPROXIMATION WITH A NEUMANN BOUND- C... ARY CONDITON OF THE FORM C... C... A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*UX(I) + F*U(I) C... C... WILL BE USED. THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... TO DROP UX C... C... A + 2*B + 3*C + 4*D + E = 0 (25) C... C... TO RETAIN UXX C... C... A + 4*B + 9*C + 16*D = 2 (26) C... C... TO DROP UXXX C... C... A + 8*B + 27*C + 64*D = 0 (27) C... C... TO DROP UXXXX C... C... A + 16*B + 81*C + 256*D = 0 (28) C... C... TO DROP UXXXXX C... C... A + 32*B + 243*C + 1024*D = 0 (29) C... C... EQUATIONS (25) TO (29) CAN BE SOLVED FOR A, B, C, D AND E. IF C... C... EQUATION (26) IS SUBTRACTED FROM EQUATIONS (27), (28) AND (29), C... C... 4*B + 18*C + 48*D = -2 (30) C... C... 12*B + 72*C + 240*D = -2 (31) C... C... 28*B + 234*C + 1008*D = -2 (32) C... C... EQUATIONS (30), (31) AND (32) CAN BE SOLVED FOR B, C AND D C... C... 18*C + 96*D = 4 (33) C... C... 108*C + 672*D = 12 (34) C... C... EQUATIONS (3) AND (34) CAN BE SOLVED FOR C AND D, C = 8/9, D = C... -1/8. C... C... FROM EQUATION (30), B = -3. FROM EQUATION (26), A = 8. FROM C... EQUATION (25), E = -25/6. C... C... THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS C... C... 8*U(I+1) - 3*U(I+2) + (8/9)*U(I+3) - (1/8)*U(I+4) C... C... - (25/6)*UX(I)*DX C... C... = (8 - 3 + (8/9) - (1/8))*U(I) + UXX*(DX**2) + O(DX**6) C... C... OR C... C... UXX(I) = (1/12*DX**2)*((-415/6)*U(I) + 96*U(I+1) - 36*U(I+2) C... (35) C... + (32/3)*U(I+3) - (3/2)*U(I+4) - 50*UX(I)*DX) + O(DX**4) C... C... EQUATION (35) WILL BE APPLIED AT I = 1 AND I = N C... C... UXX(1) = (1/12*DX**2)*((-415/6)*U(1) + 96*U(2) - 36*U(3) C... (36) C... + (32/3)*U(4) - (3/2)*U(5) - 50*UX(1)*DX) + O(DX**4) C... C... UXX(N) = (1/12*DX**2)*((-415/6)*U(N) + 96*U(N-1) - 36*U(N-2) C... (37) C... + (32/3)*U(N-3) - (3/2)*U(N-4) + 50*UX(N)*DX) + O(DX**4) C... C... CHECK THE VALUES OF A, B, C AND D IN EQUATIONS (25) TO (29) C... CALL CHECK(3) C... C... ALTERNATIVELY, FOR GRID POINT 1, AN APPROXIMATION WITH A DIRICHLET C... BOUNDARY CONDITION OF THE FORM C... C... A*U(I+1) + B*U(I+2) + C*U(I+3) + D*U(I+4) + E*U(I+5) + F*U(I) C... C... CAN BE USED. THE CORRESPONDING ALGEBRAIC EQUATIONS ARE C... C... TO DROP UX C... C... A + 2*B + 3*C + 4*D + 5*E = 0 (38) C... C... TO RETAIN UXX C... C... A + 4*B + 9*C + 16*D + 25*E = 2 (39) C... C... TO DROP UXXX C... C... A + 8*B + 27*C + 64*D + 125*E = 0 (40) C... C... TO DROP UXXXX C... C... A + 16*B + 81*C + 256*D + 625*E = 0 (41) C... C... TO DROP UXXXXX C... C... A + 32*B + 243*C + 1024*D + 3125*E = 0 (42) C... C... EQUATIONS (38), (39), (40), (41) AMD (42) CAN BE SOLVED FOR A, C... B, C, D AND E. C... C... 2*B + 6*C + 12*D + 20*E = 2 (43) C... C... 6*B + 24*C + 60*D + 120*E = 0 (44) C... C... 14*B + 78*C + 252*D + 620*E = 0 (45) C... C... 30*B + 240*C + 1020*D + 3120*E = 0 (46) C... C... EQUATIONS (43), (44), (45) AND (46) CAN BE SOLVED FOR B, C, D C... AND E C... C... 6*C + 24*D + 60*E = -6 (47) C... C... 36*C + 168*D + 480*E = -14 (48) C... C... 150*C + 840*D + 2820*E = -30 (49) C... C... EQUATIONS (47), (48) AND (49) CAN BE SOLVED FOR C, D AND E C... C... 24*D + 120*E = 22 (50) C... C... 240*D + 1320*E = 120 (51) C... C... FROM EQUATIONS (50) AND (51), D = 61/12, E = -5/6. FROM EQUATION C... (47), C = -13. FROM EQUATION (43), B = 107/6. FROM EQUATION C... (38), A = -77/6. C... C... THE FINAL DIFFERENTIATION FORMULA IS THEN OBTAINED AS C... C... (-77/6)*U(I+1) + (107/6)*U(I+2) - 13*U(I+3) + (61/12)*U(I+4) C... C... - (5/6)*U(I+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*U(I) + C... C... UXX(I)*(DX**2) + O(DX**6) C... C... OR C... C... UXX(I) = (1/12*DX**2)*(45*U(I) - 154*U(I+1) + 214*U(I+2) C... (52) C... - 156*U(I+3) + 61*U(I+4) - 10*U(I+5)) + O(DX**4) C... C... EQUATION (52) WILL BE APPLIED AT I = 1 AND I = N C... C... UXX(1) = (1/12*DX**2)*(45*U(1) - 154*U(2) + 214*U(3) C... (53) C... - 156*U(4) + 61*U(5) - 10*U(6)) + O(DX**4) C... C... UXX(N) = (1/12*DX**2)*(45*U(N) - 154*U(N-1) + 214*U(N-2) C... (54) C... -156*U(N-3) + 61*U(N-4) - 10*U(N-5)) + O(DX**4) C... C... CHECK THE VALUES OF A, B, C, D, AND E IN EQUATIONS (38) TO (42) C... CALL CHECK(4) C... C... ****************************************************************** C... C... GRID SPACING DX=(XU-XL)/DFLOAT(N-1) C... C... 1/(12*DX**2) FOR SUBSEQUENT USE R12DXS=1./(12.0D0*DX**2) C... C... UXX AT THE LEFT BOUNDARY C... C... WITHOUT UX (EQUATION (53)) IF(NL.EQ.1)THEN UXX(1)=R12DXS* 1 ( 45.0D0*U(1) 2 -154.0D0*U(2) 3 +214.0D0*U(3) 4 -156.0D0*U(4) 5 +61.0D0*U(5) 6 -10.0D0*U(6)) C... C... WITH UX (EQUATION (36)) ELSE IF(NL.EQ.2)THEN UXX(1)=R12DXS* 1 (-415.0D0/6.0D0*U(1) 2 +96.0D0*U(2) 3 -36.0D0*U(3) 4 +32.0D0/3.0D0*U(4) 5 -3.0D0/2.0D0*U(5) 6 -50.0D0*UX(1)*DX) END IF C... C... UXX AT THE RIGHT BOUNDARY C... C... WITHOUT UX (EQUATION (54)) IF(NU.EQ.1)THEN UXX(N)=R12DXS* 1 ( 45.0D0*U(N ) 2 -154.0D0*U(N-1) 3 +214.0D0*U(N-2) 4 -156.0D0*U(N-3) 5 +61.0D0*U(N-4) 6 -10.0D0*U(N-5)) C... C... WITH UX (EQUATION (37)) ELSE IF(NU.EQ.2)THEN UXX(N)=R12DXS* 1 (-415.0D0/6.0D0*U(N ) 2 +96.0D0*U(N-1) 3 -36.0D0*U(N-2) 4 +32.0D0/3.0D0*U(N-3) 5 -3.0D0/2.0D0*U(N-4) 6 +50.0D0*UX(N )*DX) END IF C... C... UXX AT THE INTERIOR GRID POINTS C... C... I = 2 (EQUATION (23)) UXX(2)=R12DXS* 1 ( 10.0D0*U(1) 2 -15.0D0*U(2) 3 -4.0D0*U(3) 4 +14.0D0*U(4) 5 -6.0D0*U(5) 6 +1.0D0*U(6)) C... C... I = N-1 (EQUATION (24)) UXX(N-1)=R12DXS* 1 ( 10.0D0*U(N ) 2 -15.0D0*U(N-1) 3 -4.0D0*U(N-2) 4 +14.0D0*U(N-3) 5 -6.0D0*U(N-4) 6 +1.0D0*U(N-5)) C... C... I = 3, 4,..., N-2 (EQUATION (9)) DO 1 I=3,N-2 UXX(I)=R12DXS* 1 ( -1.0D0*U(I-2) 2 +16.0D0*U(I-1) 3 -30.0D0*U(I ) 4 +16.0D0*U(I+1) 5 -1.0D0*U(I+2)) 1 CONTINUE RETURN END