SUBROUTINE INITAL C... C... ONE-DIMENSIONAL CONVECTION AND DIFFUSION C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... SPATIAL DIFFERENTIATOR C... C... NCASE = 1 - DSS020 C... C... NCASE = 2 - DSS024 C... NCASE=1 C... C... SYSTEM LENGTH ZL=1.0D0 C... C... SYSTEM PARAMETERS PE=10.0D0 C... C... INITIAL CONDITIONS DO 1 I=1,NZ C(I)=0.0D0 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV C... C... INITIALIZE COUNTERS FOR PRINTING AND PLOTTING IP=0 NP=81 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... CZ IF(NCASE.EQ.1)CALL DSS020(0.0D0,ZL,NZ,C,CZ,1.0D0) IF(NCASE.EQ.2)CALL DSS024(0.0D0,ZL,NZ,C,CZ,1.0D0) C... C... BOUNDARY CONDITION, Z = 0 NL=2 CZ(1)=PE*(C(1)-FB(T)) C... C... BOUNDARY CONDITION, Z = 1 NU=2 CZ(NZ)=0.0D0 C... C... CZZ CALL DSS044(0.0D0,ZL,NZ,C,CZ,CZZ,NL,NU) C... C... PDE DO 1 I=1,NZ CT(I)=CZZ(I)-PE*CZ(I) 1 CONTINUE RETURN END DOUBLE PRECISION FUNCTION FB(T) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... UNIT STEP FUNCTION IF(T.LT.0.0D0)THEN FB=0.0D0 ELSE +IF(T.GE.0.0D0)THEN FB=1.0D0 END IF RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... MONITOR THE OUTPUT WRITE(*,*)IP C... C... WRITE NUMERICAL AND SERIES SOLUTIONS EVERY 40TH CALL TO PRINT IP=IP+1 IF (((IP-1)/40*40).EQ.(IP-1))THEN WRITE(NO,1)T,(C(I),I=1,NZ) 1 FORMAT(/,' t = ',F5.2,/,F10.5,/,(5F10.5)) END IF C... C... MATLAB PLOTTING OF THE SOLUTION CALL PLOTM C... C... SERIES SOLUTION AT END OF NUMERICAL SOLUTION IF(IP.EQ.NP)CALL SERIES(NI,NO) RETURN END SUBROUTINE PLOTM IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... DEFINE A FILE FOR MATLAB PLOTTING OF THE NUMERICAL SOLUTION IF(IP.EQ.1)OPEN(1,FILE='ma2b.mat') C... C... WRITE SOLUTION AT Z = 1 FOR MATLAB PLOTTING (EVERY FIFTH POINT IS C... USED SO THAT THE INDIVIDUAL PLOTTED POINTS CAN BE DISTINGUISHED) IF(((IP-1)/5*5).EQ.(IP-1))WRITE(1,1)T,C(NZ) 1 FORMAT(2F10.4) RETURN END SUBROUTINE SERIES(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... DEFINE A FILE FOR MATLAB PLOTTING OF THE ANALYTICAL SOLUTION OPEN(2,FILE='ma2a.out') C... C... NUMBER OF EIGENVALUES NS=10 C... C... NUMBER OF PANELS IN THE NUMERICAL QUADRATURE (INTEGRATION) NQ=1000 C... C... LIMITS OF INTEGRATION ZL=0.0D0 ZU=1.0D0 C... C... SYSTEM LENGTH ZLU=1.0D0 C... C... SYSTEM PARAMETERS PE=10.0D0 C... C... COMPUTE AND PRINT NS EIGENVALUES CALL EIGVAL WRITE(NO,1) WRITE( *,1) 1 FORMAT(//,' Eigenvalues',/) WRITE(NO,2)(I,ZMU(I),I=1,NS) 2 FORMAT(I5,F10.5) C... C... TEST ORTHOGONALITY OF EIGENFUNCTIONS BY COMPUTING THE INTEGRAL C... OF THE PRODUCT OF TWO EIGENFUNCTIONS FOR I AND J = 1, 2,..., NS WRITE(NO,13) WRITE( *,13) 13 FORMAT(//,' Product integrals',/) DO 10 I=1,NS DO 11 J=1,NS CALL ORTHO(I,J,PIJ) WRITE(NO,12)I,J,PIJ 12 FORMAT(' i = ',I3,' j = ',I3,' Pij = ',F10.5) IF(J.EQ.NS)WRITE(NO,14) 14 FORMAT(/) 11 CONTINUE 10 CONTINUE C... C... COMPUTE AND PRINT NS FOURIER COEFFICIENTS NS=200 CALL EIGVAL WRITE(NO,22) WRITE( *,22) 22 FORMAT(//,' Fourier coefficients',/) DO 20 I=1,NS CALL COEFF(I) WRITE(NO,21)I,CF(I) 21 FORMAT(I5,F10.7) 20 CONTINUE C... C... FOURIER SERIES FOR UNITY WITH NS TERMS AT A SERIES OF VALUES OF Z NS=200 WRITE(NO,30) WRITE( *,30) 30 FORMAT(//,' Fourier series for unity',/) C... C... STEP THROUGH NZV VALUES OF Z NZV=11 DO 32 J=1,NZV C... C... SET Z Z=DFLOAT(J-1)/DFLOAT(NZV-1) C... C... COMPUTE AND PRINT SERIES FOR UNITY ONE=UNITY(Z) WRITE(NO,31)Z,ONE 31 FORMAT(F10.2,F10.5) 32 CONTINUE C... C... EVALUATE AND PRINT THE SERIES SOLUTION AT Z = 1 FOR A SERIES C... OF VALUES OF T USING NS TERMS IN THE SERIES WRITE(NO,42) WRITE( *,42) 42 FORMAT(//,' Series solution at z = 1, phi1(1,t)',/) Z=1.0D0 DT=0.0025D0 NTM=81 DO 40 I=1,NTM C... C... T AND SERIES SOLUTION T=DT*DFLOAT(I-1) UZT=PHI1(Z,T) C... C... PRINT SOLUTION WRITE(NO,41)Z,T,UZT 41 FORMAT(F10.3,F10.4,F10.5) C... C... WRITE SOLUTION FOR MATLAB PLOTTING WRITE(2,43)T,UZT 43 FORMAT(2F10.4) 40 CONTINUE RETURN END SUBROUTINE EIGVAL IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... FUNCTION FOR THE TRANSCENDENTAL EQUATION THAT DEFINES THE C... EIGENVALUES EXTERNAL F1 C... C... PI PI=4.0*DATAN(1.0D0) C... C... ERROR TOLERANCE EPS=1.0D-06 C... C... COMPUTE NS EIGENVALUES DO 1 I=1,NS C... C... BOUND ROOT (EIGENVALUE) ZP=DFLOAT(I-1)*PI ZN=DFLOAT(I )*PI C... C... COMPUTE ROOT AND STORE AS AN EIGENVALUE ZMU(I)=ROOT(F1,ZN,ZP,EPS) C... C... NEXT EIGENVALUE 1 CONTINUE RETURN END SUBROUTINE ORTHO(I,J,PIJ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... EXTERNAL THE FUNCTION FOR THE INTEGRAND OF THE PRODUCT INTEGRAL EXTERNAL F2 C... C... SET THE INDICES OF THE EIGENFUNCTIONS IN=I JN=J C... C... PERFORM THE INTEGRATION OF THE PRODUCT OF TWO EIGENFUNCTIONS C... FROM ZL = 0 TO ZU = 1 PIJ=SIMP(F2,ZL,ZU,NQ) RETURN END SUBROUTINE COEFF(I) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... EXTERNAL THE FUNCTION FOR THE INTEGRAND OF THE NUMERATOR AND C... DENOMINATOR INTEGRALS OF THE FOURIER COEFFICIENT EXTERNAL F3, F4 C... C... SET THE INDEX OF THE FOURIER COEFFICIENT IN=I C... C... PERFORM THE INTEGRATION FOR THE NUMERATOR AND DENOMINATOR C... INTEGRALS FROM ZL = 0 TO ZU = 1 ZNUM=SIMP(F3,ZL,ZU,NQ) ZDEN=SIMP(F4,ZL,ZU,NQ) C... C... FOURIER COEFFICIENT CF(I)=ZNUM/ZDEN RETURN END DOUBLE PRECISION FUNCTION UNITY(Z) C... C... FUNCTION UNITY COMPUTES THE SERIES EXPANSION OF UNITY (THE C... NUMBER ONE) AT Z C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... INITIALIZE SERIES SUM=0.0D0 C... C... PRECOMPUTE THE COMMON EXPONENTIAL FACTOR EXPZ=DEXP((PE/2.0D0)*Z) C... C... COMPUTE SERIES WITH NS TERMS DO 1 I=1,NS C... C... EIGENFUNCTION PSI(I) PSII=PE/(2.0D0*ZMU(I))*DSIN(ZMU(I)*Z)+DCOS(ZMU(I)*Z) C... C... RUNNING SUM SUM=SUM+CF(I)*PSII C... C... CONTINUE SUMMATION 1 CONTINUE C... C... INCLUDE EXPONENTIAL FACTOR UNITY=EXPZ*SUM C... C... CALCULATION IS COMPLETE RETURN END DOUBLE PRECISION FUNCTION PHI1(Z,T) C... C... FUNCTION PHI1 COMPUTES THE SERIES SOLUTION AT Z AND T C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ TM, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... INITIALIZE SERIES SUM=0.0D0 C... C... PRECOMPUTE THE COMMON EXPONENTIAL FACTOR EXPZ=DEXP((PE/2.0D0)*Z) C... C... COMPUTE SERIES WITH NS TERMS DO 1 I=1,NS C... C... EXPONENTIAL IN T TLAM=ZMU(I)**2+PE**2/4.0D0 C... C... TEST MAGNITUDE OF TLAM*T TO AVOID UNDERFLOW OF EXP FUNCTION TLAM=-TLAM*T IF(TLAM.LT.-1.0D+20)THEN C... C... TERMS ARE NEGLIGIBLY SMALL SO TERMINATE SUMMING SERIES GO TO 2 ELSE C... C... CONTINUE TO SUM SERIES EXPT=DEXP(TLAM) END IF C... C... EIGENFUNCTION PSI(I) PSII=PE/(2.0D0*ZMU(I))*DSIN(ZMU(I)*Z)+DCOS(ZMU(I)*Z) C... C... RUNNING SUM SUM=SUM+CF(I)*(1.0D0-EXPT)*PSII C... C... CONTINUE SUMMATION 1 CONTINUE C... C... INCLUDE EXPONENTIAL FACTOR 2 PHI1=EXPZ*SUM C... C... CALCULATION IS COMPLETE RETURN END DOUBLE PRECISION FUNCTION F1(Z) C... C... F1 COMPUTES THE TRANSCENDENTAL FUNCTION THAT DEFINES THE C... EIGENVALUES. F1 IS CALL BY FUNCTION ROOT VIA SUBROUTINE C... EIGVAL TO FIND THE ZEROS OF THE TRANSCENDENTAL FUNCTION. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... COMPUTE FUNCTION TO BE ZEROED F1=1.0D0/DTAN(Z)-Z/PE+PE/(4.0D0*Z) RETURN END DOUBLE PRECISION FUNCTION F2(Z) C... C... FUNCTION F2 COMPUTES THE PRODUCT OF TWO EIGENFUNCTIONS. F2 C... IS CALLED BY FUNCTION SIMP VIA SUBROUTINE ORTHO. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... EIGENFUNCTION PSI(I) PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z) C... C... EIGENFUNCTION PSI(J) PSIJ=PE/(2.0D0*ZMU(JN))*DSIN(ZMU(JN)*Z)+DCOS(ZMU(JN)*Z) C... C... PRODUCT OF THE EIGENFUNCTIONS F2=PSII*PSIJ RETURN END DOUBLE PRECISION FUNCTION F3(Z) C... C... FUNCTION F3 COMPUTES THE NUMERATOR INTEGRAND OF THE FOURIER C... COEFFICIENT. F3 IS CALLED BY FUNCTION SIMP VIA SUBROUTINE C... COEFF. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... EIGENFUNCTION PSI(I) PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z) C... C... PRODUCT OF EXPONENTIAL AND EIGENFUNCTION F3=DEXP((-PE/2.0D0)*Z)*PSII RETURN END DOUBLE PRECISION FUNCTION F4(Z) C... C... FUNCTION F4 COMPUTES THE DENOMINATOR INTEGRAND OF THE FOURIER C... COEFFICIENT. F4 IS CALLED BY FUNCTION SIMP VIA SUBROUTINE C... COEFF. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NZ=81,NA=200) COMMON/T/ T, NSTOP, NORUN + /Y/ C(NZ) + /F/ CT(NZ) + /S/ CZ(NZ), CZZ(NZ) + /C/ ZL, ZU, ZLU, PE, + ZMU(NA), CF(NA) + /I/ NS, NQ, NP, IP, + NCASE, IN, JN C... C... EIGENFUNCTION PSI(I) PSII=PE/(2.0D0*ZMU(IN))*DSIN(ZMU(IN)*Z)+DCOS(ZMU(IN)*Z) C... C... EIGENFUNCTION SQUARED F4=PSII**2 RETURN END DOUBLE PRECISION FUNCTION ROOT(F,ZN,ZP,EPS) C... C... FUNCTION ROOT COMPUTES THE ROOT ON AN EQUATION (ZERO OF A C... FUNCTION) BY THE SECANT METHOD. THE CODE WAS ORIGINALLY C... PROVIDED BY YOUNG-SOON CHOE, AND MODIFIED SLIGHTLY BY W. E. C... SCHIESSER. C... C... ARGUMENTS C... C... F EXTERNAL TO COMPUTE THE FUNCTION TO BE ZEROED (INPUT) C... C... ZN LOWER BOUND ON THE ROOT (INPUT) C... C... ZP UPPER BOUND ON THE ROOT (INPUT) C... C... EPS ABSOLUTE ERROR TOLERANCE FOR THE ROOT (INPUT) C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... EXTERNAL FOR THE FUNCTION TO BE ZEROED EXTERNAL F C... C... BEGIN SEARCH FOR THE ROOT BY THE SECANT METHOD IP=0 IN=0 C... C... INITIAL ESTIMATE OF THE ROOT 2 Z=(ZP+ZN)/2.0D0 C... C... COMPUTE FUNCTION TO BE ZEROED FZ=F(Z) C... C... CONTINUE SEARCH IF(FZ.LE.0.0D0)THEN ZN=Z FM=FZ IN=1 ELSE ZP=Z FP=FZ IP=1 END IF IF((IP*IN).EQ.0) GO TO 2 IF((FP-FM).GT.0.1D0)GO TO 2 C... C... UPDATE Z AND COMPUTE FUNCTION 3 SL=(FP-FM)/(ZP-ZN) Z=ZP-FP/SL FZ=F(Z) C... C... TEST FOR CONVERGENCE IF((DABS(FZ).LT.EPS).OR.(DABS(ZP-ZN).LT.EPS))THEN C... C... CONVERGENCE ACHIEVED SO STORE ROOT AND RETURN TO CALLING C... PROGRAM ROOT=Z C... C... CONTINUE SEARCH ELSE IF(FZ.GE.0.0D0)THEN FM=FZ ZN=Z GO TO 3 ELSE FP=FZ ZP=Z GO TO 3 END IF END IF RETURN END DOUBLE PRECISION FUNCTION SIMP(F,A,B,NQ) C... C... FUNCTION SIMP PERFORMS A NUMERICAL QUADRATURE (INTEGRATION) BY C... SIMPSON'S RULE C... C... ARGUMENTS C... C... F EXTERNAL TO COMPUTE THE INTEGRAND (INPUT) C... C... A LOWER LIMIT OF THE INTEGRAL (INPUT) C... C... B UPPER LIMIT OF THE INTEGRAL (INPUT) C... C... NQ NUMBER OF PANELS USED IN SIMPSON'S RULE (MUST BE C... EVEN) (INPUT) C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... EXTERNAL FUNCTION FOR INTEGRAND EXTERNAL F C... C... CHECK FOR AN EVEN NUMBER OF PANELS IF((NQ/2*2).NE.NQ)THEN WRITE(*,2)NQ 2 FORMAT(//,' NQ = ',I3,' (NOT EVEN)') SIMP=0.0D0 RETURN END IF C... C... INTEGRATION INTERVAL DZ=(B-A)/DFLOAT(NQ) C... C... INITIALIZE SUM SIMP=(F(A)+4.0D0*F(A+DZ)+F(B)) C... C... CONTINUE SUM IF(NQ.GT.2)THEN DO 1 I=2,NQ-1,2 SIMP=SIMP + +2.0D0*F(A+DFLOAT(I )*DZ) + +4.0D0*F(A+DFLOAT(I+1)*DZ) 1 CONTINUE END IF C... C... INTEGRAL SIMP=(DZ/3.0D0)*SIMP RETURN END