*DECK MATH42 SUBROUTINE INITAL C... C... AN INTRODUCTION TO FINITE INTEGRAL TRANSFORMS C... C... THE FOLLOWING DEVELOPMENT ILLUSTRATES THE DERIVATION OF A FINITE C... INTEGRAL TRANSFORM IN TERMS OF AN APPLICATION. THE INTENT IS TO C... DEMONSTRATE HOW THESE TRANSFORMS CAN BE DEVELOPED AND USED IN C... THE SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS. C... C... CONSIDER THE LINEAR, PARABOLIC PARTIAL DIFFERENTIAL EQUATION (PDE) C... I.E., FOURIER'S SECOND LAW C... C... U = U (1) C... T XX C... C... THE BOUNDARY CONDITIONS FOR EQUATION (1) ARE TAKEN AS C... C... U(0,T) = 0, U (XL,T) = 0 (2)(3) C... X C... C... AND THE INITIAL CONDITION AS C... C... U(X,0) = 1 (4) C... C... NOTE THAT BOUNDARY CONDITION (2) AND INITIAL CONDITION (4) ARE C... INCONSISTENT, I.E., U(0,0) IS AMBIGUOUS. FORTUNATELY, DISCON- C... TINUITIES OF THIS TYPE (FROM U(X,0) = 1 TO U(0,T) = 0 AT X = T C... = 0) ARE GENERALLY DAMPED BY PARABOLIC PDES AND HAVE LITTE EFFECT C... AFTER T = 0. THIS IS NOT TRUE FOR HYPERBOLIC PDES, HOWEVER, C... WHICH TEND TO PROPAGATE DISCONTINUITIES. C... C... NUMERICAL AND ANALYTICAL SOLUTIONS TO EQUATIONS (1) TO (4) ARE C... COMPUTED AND COMPARED IN THE FOLLOWING CODING. IN THIS CASE, C... THE NUMERICAL SOLUTION MAY BE MORE ACCURATE BECAUSE OF THE SLOW C... CONVERGENCE OF THE INFINITE SERIES IN THE ANALYTICAL SOLUTION C... (THE FOURIER COEFFICIENTS ARE PROPORTIONAL TO 1/N WHERE N IS THE C... SUMMATION INDEX BECAUSE OF THE DISCONTINUITY) AND THE GIBBS C... PHENOMENON NEAR X = 0. C... C... THE ANALYTICAL SOLUTION TO EQUATIONS (1) TO (4) CAN BE DERIVED C... BY CLASSICAL FOURIER SERIES ANALYSIS. A SOLUTION WHICH SATISFIES C... THE PDE, EQUATION (1), AND THE TWO BOUNDARY CONDITIONS, EQUATIONS C... (2) AND (3), IS OBTAINED BY STRAIGHTFORWARD SEPARATION OF VARI- C... ABLE AS C... C... INF C... U(X,T) = SUM CN*EXP(-(EN**2)*T)*SIN(EN*X) (5) C... 1 C... C... WHERE THE EIGENVALUES ARE EN = (N - 1/2)*PI/XL, N = 1, 2,... C... APPLICATION OF THE INITIAL CONDITION, EQUATION (4), TO EQUATION C... (5) REQUIRES FOR THE FOURIER EXPANSION OF 1 C... C... INF C... 1 = SUM CN*SIN(EN*X) (6) C... 1 C... C... MULTIPLICATION OF BOTH SIDES OF EQUATION (6) BY SIN(EM*X), AND C... INTEGRATION FROM X = 0 TO X = XL, KEEPING IN MIND THE ORTHOGON- C... ALITY OF THE SINE FUNCTIONS (SO THAT ONLY THE CASE N = M NEED BE C... CONSIDERED) C... C... XL XL 2 C... INT 1*SIN(EN*X)*DX = CN*INT SIN (EN*X)*DX (7) C... 0 0 C... C... THE LHS INTEGRAL OF EQUATION (7) IS EASILY EVALUATED AS 1/EN. C... THE RHS INTEGRAL OF EQUATION (7) IS EASILY EVALUATED BY THE C... WELL KNOW TRIGONOMETRIC IDENTITY C... C... 2 C... SIN (EN*X) = (1/2)*(1 - COS (2*EN*X)) C... C... THUS C... C... XL 2 XL C... INT SIN (EN*X)*DX = (1/2) INT (1 - COS(2*EN*X))*DX C... 1 1 C... C... XL C... . C... = (1/2)*(X - 1/(2*EN)*SIN(2*EN*X)). C... . C... 0 C... C... = (1/2)*(XL - 1/(2*EN)*SIN(2*(N - 1/2)*PI)) = (1/2)*XL C... C... C... THEREFORE, THE FOURIER COEFFICIENTS ARE C... C... CN=2/(EN*XL) (8) C... C... EQUATIONS (5) AND (8) ARE THE ANALYTICAL SOLUTION. C... C... AN ALTERNATIVE APPROACH TO THIS CONVENTIONAL SEPARATION OF C... SERIES SOLUTION IS TO DEFINE A FINITE INTEGRAL TRANSFORM, G(EN), C... OF A FUNCTION F(X) C... C... XL C... G(EN) = INT F(X)*SIN(EN*X)*DX (9) C... 0 C... C... WITH THE INVERSE TRANSFORM C... C... INF C... F(X) = (2/XL) SUM G(EN)*SIN(EN*X) (10) C... N=1 C... C... WHERE, AGAIN, EN = PI*(N - 1/2)/XL, N = 1, 2,... C... C... USING THE DEFINITION OF THE TRANSFORM, EQUATION (9), THE C... FOLLOWING TRANSFORM PAIRS CAN BE DERIVED C... C... F(X) G(EN) C... C... A*F(X) A*G(EN) C... C... A*F1(X) + B*F2(X) A*G1(EN) + B*G2(EN) C... C... 1 1/EN C... C... F''(X) (-1**(N+1))*F'(XL) C... C... + EN*F(0) C... C... - (EN**2)*G(EN) C... C... THE TRANSFORM OF F''(X), WHICH IS CENTRAL TO THE SOLUTION OF C... SECOND-ORDER PDES, IS OBTAINED IN THE FOLLOWING WAY C... C... XL C... G(EN) = INT F''(X)*SIN(EN*X)*DX C... 0 C... C... XL C... . XL C... = F'(X)*SIN(EN*X). - EN INT F'(X)*COS(EN*X)*DX C... . 0 C... 0 C... C... = (-1**(N+1))*F'(XL) C... C... XL C... . XL C... - EN*(F(X)*COS(EN*X). + (EN**2) INT F(X)*SIN(EN*X)DX C... . 0 C... 0 C... C... = (-1**(N+1))*F'(XL) + EN*F(0) - (EN**2)*G(EN) C... C... THIS RESULT SUGGESTS THAT TRANSFORM PAIRS (9) AND (10) CAN BE C... USED IN THE SOLUTION OF SECOND-ORDER PDES WITH A DIRICHLET C... BOUNDARY CONDITION AT X = 0 (DUE TO THE TERM EN*F(0)) AND A C... NEUMANN BOUNDARY CONDITION AT X = XL (DUE TO THE TERM (-1**(N+1) C... *F'(XL)). C... C... SPECIFICALLY, FOR THE PROBLEM DEFINED BY EQUATIONS (1) TO (4), C... IF A FINITE INTEGRAL TRANSFORM OF U(X,T) WITH RESPECT TO X IS C... DEFINED AS V(EN,T), EQUATIONS (1) TO (4) TRANSFORM TO C... C... V = (EN**2)*V (11) C... T C... C... V(EN,0) = 1/EN (12) C... C... THE SOLUTION TO EQUATIONS (11) AND (12) IS EASILY OBTAINED AS C... C... V(EN,T) = (1/EN)*EXP(-(EN**2)*T) C... C... WHICH, WHEN SUBSTITUTED IN THE INVERSION FORMULAS, EQUATION (10), C... GIVES THE PREVIOUS SEPARATION OF VARIABLES SOLUTION, EQUATIONS C... (5) AND (8) C... C... INF C... U(X,T) = (2/XL) SUM (1/EN)*EXP(-(EN**2)*T)*SIN(EN*X) (13) C... 1 C... C... THE PRECEDING DISCUSSION OF THE FINITE INTEGRAL TRANSFORM C... SUGGESTS A GENERALIZATION LEADING TO FOUR TRANSFORMS, TWO FOR C... THE SINE KERNEL FUNCTION SIN(EN*X), WITH EN = PI*N/XL OR C... EN = PI*(N - 1/2)/XL, AND TWO FOR THE COSINE KERNEL FUNCTION C... COS(EN*X), WITH EN = PI*N/XL OR EN = PI*(N - 1/2)/XL. C... C... FOR THE SINE FUNCTION, THE TRANSFORM OF THE SECOND DERIVATIVE C... OF A FUNCTION F(X) IS C... C... INT F''(X)*SIN(EN*X)*DX = C... C... XL C... . XL C... F'(X)*SIN(EN*X). - EN INT F'(X)*COS(EN*X)*DX C... . 0 C... 0 C... C... = F'(XL)*SIN(EN*XL) C... C... XL C... . XL C... - EN*(F(X)*COS(EN*X). + EN INT F(X)*SIN(EN*X)*DX) C... . 0 C... 0 C... C... = F'(XL)*SIN(EN*XL) - EN*F(XL)*COS(EN*XL) + EN*F(0) C... C... - (EN**2)*G(EN) C... C... FOR EN = PI*N/XL, THIS RESULT REDUCES TO C... C... -EN*F(XL)*(-1**N) + EN*F(0) - (EN**2)*G(EN) (14) C... C... WHILE FOR EN = PI*(N - 1/2)/XL, WE OBTAIN THE PREVIOUS RESULT C... C... F'(XL)*(-1**(N+1)) + EN*F(0) - (EN**2)*G(EN) (15) C... C... SIMILARLY, FOR THE COSINE FUNCTION C... C... XL C... INT F''(X)*COS(EN*X)*DX = C... 0 C... C... XL C... . XL C... F'(X)*COS(EN*X). + EN INT F'(X)*SIN(EN*X)*DX C... . 0 C... 0 C... C... = F'(XL)*COS(EN*XL) - F'(0) C... C... XL C... . XL C... + EN*(F(X)*SIN(EN*X). - EN INT F(X)*COS(EN*X)*DX) C... . 0 C... 0 C... C... = F'(XL)*COS(EN*XL) - F'(0) C... C... + EN*F(XL)*SIN(EN*XL) - (EN**2)*G(EN) C... C... FOR EN = PI*N/XL, THIS RESULT REDUCES TO C... C... F'(XL)*(-1**N) - F'(0) - (EN**2)*G(EN) (16) C... C... WHILE FOR EN = PI*(N - 1/2)/XL, C... C... - F'(0) + EN*F(XL)*(-1**N) - (EN**2)*G(EN) (17) C... C... THE FOUR TRANSFORMS CAN BE SUMMARIZED AS C... C... C... F(X) G(EN) C... C... INT F(X)*SIN(EN*X)*DX (2/XL) SUM G(EN)*SIN(EN*X) C... C... EN = PI*N/XL C... C... F''(X) - EN*F(XL)*(-1**N) + EN*F(0) C... C... - (EN**2)*G(EN) C... C... USE WITH DIRICHLET BC AT X = 0 AND DIRICHLET BC AT X = XL C... C... C... EN = PI*(N - 1/2)/XL C... C... F''(X) F'(XL)*(-1**(N+1)) + EN*F(0) C... C... - (EN**2)*G(EN) C... C... USE WITH DIRICHLET BC AT X = 0 AND NEUMANN BC AT X = XL C... C... C... XL C... INT F(X)*COS(EN*X)*DX (2/XL) SUM G(EN)*COS(EN*X) C... 0 C... C... EN = PI*N/XL C... C... F''(X) F'(XL)*(-1**N) - F'(0) C... C... - (EN**2)*G(EN) C... C... USE WITH NEUMANN BC AT X = 0 AND NEUMANN BC AT X = XL C... C... C... EN = PI*(N - 1/2)/XL C... C... F''(X) - F'(0) + EN*F(XL)*(-1**N) C... C... - (EN**2)*G(EN) C... C... USE WITH NEUMANN BC AT X = 0 AND DIRICHLET BC AT X = XL C... C... C... OF COURSE, MANY OTHER FINITE INTEGRAL TRANSFORMS ARE POSSIBLE C... THROUGH THE USE OF VARIOUS ORTHOGONAL FUNCTIONS, E.G., THE C... SO-CALLED 'HANKEL TRANSFORMS' BASED ON BESSEL FUNCTIONS. THE C... PRINCIPAL ADVANTAGES OF FINITE TRANSFORMS ARE C... C... (1) THEY SYSTEMATIZE THE SEPARATION OF VARIABLES SOLUTION C... OF LINEAR PARTIAL DIFFERENTIAL EQUATIONS. C... C... (2) THE INVERSION OF THE TRANSFORM IS PARTICULAR EASY (THE C... DESCRIPTION 'FINITE (INTEGRAL) TRANSFORM' REFERS TO C... INVERSION WHICH REQUIRES ONLY A SUBSTITUTION INTO AN C... INFINITE SERIES, RATHER THAN THE EVALUATION OF AN C... INVERSION INTEGRAL, E.G., A FOURIER OR LAPLACE TRANS- C... FORM). C... C... GENERALLY THOUGH, FINITE INTEGRAL TRANSFORMS ARE LIMITED TO C... PROBLEMS FOR WHICH SEPARATION OF VARIABLES CAN BE APPLIED, C... I.E., THEY CANNOT GO BEYOND WHAT CAN BE ACCOMPLISHED WITH THE C... SEPARATION OF VARIABLES. C... C... THE FOLLOWING CODING IMPLEMENTS A NUMERICAL METHOD OF LINES C... SOLUTION TO EQUATIONS (1) TO (4), WHICH IS COMPARED WITH THE C... ANALYTICAL SOLUTION, EQUATIONS (5) AND (8) (OR EQUATION (13)). C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX),UXX(NX) 4 /C/ XL, IP C... C... TOTAL LENGTH XL=1.0D0 C... C... INITIAL CONDITION DO 1 I=1,NX U(I)=1.0D0 1 CONTINUE RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ T, NSTOP, NORUN C... C... STAGEWISE DIFFERENTIATION IF(NORUN.EQ.1)CALL DERV1 C... C... DIRECT DIFFERENTIATION IF(NORUN.EQ.2)CALL DERV2 RETURN END SUBROUTINE DERV1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX),UXX(NX) 4 /C/ XL, IP C... C... BOUNDARY CONDITION AT X = 0 U(1)=0.0D0 C... C... UX CALL DSS004(0.0D0,XL,NX,U,UX) C... C... BOUNDARY CONDITION AT X = XL UX(NX)=0.0D0 C... C... UXX CALL DSS004(0.0D0,XL,NX,UX,UXX) C... C... PDE DO 1 I=1,NX UT(I)=UXX(I) 1 CONTINUE RETURN END SUBROUTINE DERV2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX),UXX(NX) 4 /C/ XL, IP C... C... BOUNDARY CONDITION AT X = 0 NL=1 U(1)=0.0D0 C... C... BOUNDARY CONDITION AT X = XL NU=2 UX(NX)=0.0D0 C... C... UXX CALL DSS044(0.0D0,XL,NX,U,UX,UXX,NL,NU) C... C... PDE DO 1 I=1,NX UT(I)=UXX(I) 1 CONTINUE RETURN END DOUBLE PRECISION FUNCTION UANAL(X,T,XL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... SUM 200 TERMS IN THE SERIES SOLUTION. THIS LARGE NUMBER WAS C... SELECTED BECAUSE OF THE SLOW CONVERGENCE OF THE SERIES, I.E., C... THE FOURIER COEFFICIENTS DECREASE AS 1/I DUE TO THE DISCONTINUITY C... AT X = 0 SUM=0.0D0 DO 1 I=1,200 C... C... EIGENVALUES EN=(DFLOAT(I)-0.5D0)*3.1415927D0/XL C... C... FOURIER COEFFICIENTS CN=2.0D0/(EN*XL) C... C... ARGUMENT OF THE EXPONENTIAL ARG=-(EN**2)*T C... C... BRANCH TO AVOID AN UNDERFLOW OF THE EXPONENTIAL FUNCTION IF(ARG.LT.-200.0D0)THEN EX=0. ELSE EX=DEXP(ARG) END IF C... C... NEXT TERM IN THE SERIES TERM=CN*EX*DSIN(EN*X) C... C... SUM THE SERIES SUM=SUM+TERM 1 CONTINUE C... C... 200 TERMS HAVE BEEN SUMMED AS THE ANALYTICAL SOLUTION UANAL=SUM RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ T, NSTOP, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX),UXX(NX) 4 /C/ XL, IP DIMENSION XA(6), UA(6) C... C... INITIALIZE A COUNTER FOR THE PRINTING AND PLOTTING DATA IP/0/ IP=IP+1 IF(IP.EQ.1)THEN C... C... MAP THE ODE JACOBIAN MATRIX AT THE BEGINNING OF THE SOLUTION CALL MAP END IF C... C... MONITOR THE CALCULATIONS WRITE(*,*)NORUN,T C... C... COMPUTE THE ANALYTICAL SOLUTION AT X = 0, 0.2*XL, 0.4*XL, 0.6*XL, C... 0.8*XL, XL XA(1)=0.D0 XA(2)=0.2D0*XL XA(3)=0.4D0*XL XA(4)=0.6D0*XL XA(5)=0.8D0*XL XA(6)=XL DO 1 I=1,6 UA(I)=UANAL(XA(I),T,XL) 1 CONTINUE C... C... PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS AT X = 0, 0.2*XL, C... 0.4*XL, 0.6*XL, 0.8*XL, XL IF(T.EQ.0.)WRITE(NO,4) 4 FORMAT(10X,'T',8X,'X=0',3X,'X=0.2*XL',3X,'X=0.4*XL', 1 3X,'X=0.6*XL',3X,'X=0.8*XL', 2 7X,'X=XL') WRITE(NO,2)T,(U(I),I=1,NX,4) 2 FORMAT(F11.2,/,' NUMERICAL',6F11.4) WRITE(NO,3)(UA(I),I=1,6) 3 FORMAT(' ANALYTICAL',6F11.4,//) C... C... MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION IF(IP.LT.11)RETURN CALL MAP C... C... RESET THE INTEGER COUNTER FOR THE NEXT RUN IP=0 RETURN END SUBROUTINE MAP PARAMETER (N=21) COMMON/T/ T 1 /Y/ Y(N) 2 /F/ F(N) C... C... SUBROUTINE MAP CALLS SUBROUTINE JMAP TO MAP THE JACOBIAN MATRIX C... OF AN NTH-ORDER ODE SYSTEM. C... C... DECLARE SELECTED VARIABLES DOUBLE PRECISION DOUBLE PRECISION T, Y, F, YOLD(N), FOLD(N) C... C... DEFINE SINGLE PRECISION ARRAY REQUIRED BY SUBROUTINE JMAP (A) REAL A(N,N) C... C... MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE C... DERV CALL JMAP(N,A,Y,YOLD,F,FOLD) RETURN END