*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
