PROGRAM H1 C... C... H1 - ERF SOLUTION OF FOURIER'S SECOND LAW; MATLAB PLOTTING C... OF THE SOLUTION C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... SPATIAL GRID PARAMETER(NX=51) DIMENSION X(NX) C... C... TIME GRID PARAMETER(NT=5) DIMENSION T(NT) C... C... EXTERNAL FUNCTION TO CALCULATE INTEGRAND C... C... FE INTEGRAND OF ERF(X) C... EXTERNAL FE C... C... OPEN OUTPUT FILE NO=7 OPEN(NO,FILE='output') C... C... OPEN A FILE FOR MATLAB PLOTTING OPEN(1,FILE='h1.out') C... C... HEAT CONDUCTION PROBLEM C... C... DEFINE SPATIAL GRID XL=1.0D0 DO 1 I=1,NX X(I)=XL*DFLOAT(I-1)/DFLOAT(NX-1) 1 CONTINUE C... C... DEFINE TIME GRID TF=500.0D0 DO 2 I=1,NT T(I)=TF*DFLOAT(I)/DFLOAT(NT) 2 CONTINUE C... C... THERMAL DIFFUSIVITY, INITIAL, SURFACE TEMPERATURES ALPHA=0.0001D0 T0= 25.0D0 TS=250.0D0 C... C... PI C... C... STEP THROUGH TIMES DO 10 I=1,NT PI=4.0D0*DATAN(1.0D0) C... C... LABEL FOR TABULATED OUTPUT WRITE(NO,6) 6 FORMAT(4X,'I',4X,'J',6X,'T(I)',6X,'X(J)', + 5X,'ERFAB',5X,'ERFHA',8X,'T1') C... C... STEP THROUGH SPACE DO 11 J=1,NX C... C... LIMITS A=0.0D0 B=X(J)/DSQRT(4.0D0*ALPHA*T(I)) C... C... EVALUATION OF ERF(X) NP=100 ERFAB=SIMP(FE,A,B,NP) C... C... INCLUDE FACTOR 2/SQRT(PI) ERFAB=2.0D0/DSQRT(PI)*ERFAB C... C... TEMPERATURE T1=TS+(T0-TS)*ERFAB C... C... ERF(X) FROM HASTINGS APPROXIMATION ERFHA=ERF(B) C... C... MONITOR CALCULATION WRITE(*,*)I,J C... C... WRITE NUMERICAL RESULTS WRITE(NO,4)I,J,T(I),X(J),ERFAB,ERFHA,T1 4 FORMAT(2I5,F10.1,F10.3,2F10.5,F10.2) C... C... WRITE THE OUTPUT FOR MATLAB PLOTTING WRITE(1,5)X(J),T1 5 FORMAT(2F10.5) C... C... NEXT SPACE POINT 11 CONTINUE C... C... NEXT TIME 10 CONTINUE C... C... CALCULATIONS COMPLETE STOP END DOUBLE PRECISION FUNCTION SIMP(F,A,B,NP) C... C... FUNCTION SIMP COMPUTES A SIMPSON'S RULE QUADRATURE C... C... ARGUMENT LIST C... C... F EXTERNAL FUNCTION TO DEFINE THE INTEGRAND (INPUT) C... C... A LOWER LIMIT OF THE INTEGRAL (INPUT) C... C... B UPPER LIMIT OF THE INTEGRAL (INPUT) C... C... NP NUMBER OF PANELS (MUST BE EVEN) (INPUT) C... C... DOUBLE PRECISION CODING IS USED 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((NP/2*2).NE.NP)THEN WRITE(*,2)NP 2 FORMAT(//,' NP = ',I3,' (NOT EVEN)') SIMP=0.0D0 RETURN END IF C... C... INTEGRATION INTERVAL DX=(B-A)/DFLOAT(NP) C... C... INITIALIZE SUM SIMP=(F(A)+4.0D0*F(A+DX)+F(B)) C... C... CONTINUE SUM IF(NP.GT.2)THEN DO 1 I=2,NP-1,2 SIMP=SIMP + +2.0D0*F(A+DFLOAT(I )*DX) + +4.0D0*F(A+DFLOAT(I+1)*DX) 1 CONTINUE END IF C... C... INTEGRAL SIMP=(DX/3.0D0)*SIMP RETURN END DOUBLE PRECISION FUNCTION F3(X) IMPLICIT DOUBLE PRECISION (A-H,O-Z) F3=X**3 RETURN END DOUBLE PRECISION FUNCTION F4(X) IMPLICIT DOUBLE PRECISION (A-H,O-Z) F4=X**4 RETURN END DOUBLE PRECISION FUNCTION FE(X) IMPLICIT DOUBLE PRECISION (A-H,O-Z) FE=DEXP(-X**2) RETURN END DOUBLE PRECISION FUNCTION DFLOAT(I) C... C... FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO A C... DOUBLE PRECISION FLOATING POINT. THIS FUNCTION IS PROVIDED C... IN CASE THE USER'S FORTRAN COMPILER DOES NOT INCLUDE DFLOAT. C... DFLOAT=DBLE(FLOAT(I)) RETURN END DOUBLE PRECISION FUNCTION ERF(X) C... C... FUNCTION ERF COMPUTES THE ERROR FUNCTION ACCORDING TO THE C... HASTINGS APPROXIMATION: C... C... GREENBERG, M. D. 1978. FOUNDATIONS OF APPLIED MATHEMATICS. C... ENGLEWOOD CLIFFS, NJ: PRENTICE-HALL. C... C... HASTINGS, C., JR. 1955. APPROXIMATIONS FOR DIGITAL COMPUTERS. C... PRINCETON, NJ: PRINCETON UNIVERSITY PRESS. C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... THE FOLLOWING CODING IS A STRAIGHTFORWARD IMPLEMENTATION OF C... THE HASTINGS APPROXIMATION CITED ABOVE T= 1.0D0/(1.0D0+0.3275911D0*X) A1= 0.254829592D0 A2=-0.284496736D0 A3= 1.421413741D0 A4=-1.453152027D0 A5= 1.061405429D0 ERF=1.0D0-(A1*T+A2*T**2+A3*T**3+A4*T**4+A5*T**5)*DEXP(-X**2) RETURN END