PROGRAM H1 C... C... H1 - ERF SOLUTION OF FOURIER'S SECOND LAW C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... EXTERNAL FUNCTIONS TO CALCULATION INTEGRANDS FOR VARIOUS CASES C... C... F3 THIRD ORDER POLYNOMIAL (CUBIC FUNCTION) C... C... F4 FOURTH ORDER POLYNOMIAL (QUARTIC FUNCTION) C... C... FE INTEGRAND OF ERF(X) C... EXTERNAL F3, F4, FE C... C... OPEN OUTPUT FILE NO=7 OPEN(NO,FILE='output') C... C... INTEGRATION LIMITS A=0.0D0 B=1.0D0 C... C... STEP THROUGH A SERIES OF RUNS FOR POLYNOMIALS WITH DIFFERING C... NUMBERS OF PANELS DO 2 NORUN=1,4 IF(NORUN.EQ.1)NP=2 IF(NORUN.EQ.2)NP=4 IF(NORUN.EQ.3)NP=10 IF(NORUN.EQ.4)NP=100 C... C... INTEGRATION OF THIRD ORDER POLYNOMIAL CI3=SIMP(F3,A,B,NP) C... C... INTEGRATION OF FOURTH ORDER POLYNOMIAL CI4=SIMP(F4,A,B,NP) C... C... PRINT NUMERICAL INTEGRALS WRITE (NO,1)NP, + A,B,CI3, + A,B,CI4 1 FORMAT(' NP = ',I3,//, + ' A = ',F6.2,' B = ',F6.2,' I3 = ',F10.7,/, + ' A = ',F6.2,' B = ',F6.2,' I4 = ',F10.7,/) 2 CONTINUE C... C... HEAT CONDUCTION PROBLEM C... C... THREE RUNS FOR T = 50, 100, 200 S (SECONDS) DO 3 NORUN=1,3 IF(NORUN.EQ.1)T= 50.0D0 IF(NORUN.EQ.2)T=100.0D0 IF(NORUN.EQ.3)T=200.0D0 C... C... THERMAL DIFFUSIVITY, INITIAL, SURFACE TEMPERATURES ALPHA=0.0001D0 T0= 25.0D0 TS=250.0D0 C... C... PI PI=4.0D0*DATAN(1.0D0) C... C... X COORDINATES X1=0.25D0 X2=1.00D0 C... C... LIMITS A=0.0D0 B1=X1/DSQRT(4.0D0*ALPHA*T) B2=X2/DSQRT(4.0D0*ALPHA*T) C... C... EVALUATION OF ERF(X) NP=100 ERFB1=SIMP(FE,A,B1,NP) ERFB2=SIMP(FE,A,B2,NP) C... C... INCLUDE FACTOR 2/SQRT(PI) ERFB1=2.0D0/DSQRT(PI)*ERFB1 ERFB2=2.0D0/DSQRT(PI)*ERFB2 C... C... TEMPERATURES AT X1, X2 T1=TS+(T0-TS)*ERFB1 T2=TS+(T0-TS)*ERFB2 C... C... PRINT NUMERICAL RESULTS WRITE(NO,4)NP,T,A,B1,ERFB1,T1,A,B2,ERFB2,T2 4 FORMAT(//,' NP = ',I3,' T = ',F6.1,//, + ' A = ',F6.2,' B1 = ',F6.2, + ' ERF(B1) = ',F10.7,' T1 = ',F10.7,/, + ' A = ',F6.2,' B2 = ',F6.2, + ' ERF(B2) = ',F10.7,' T2 = ',F10.7,/) C... C... ERF(X) FROM HASTINGS APPROXIMATION ERFB1=ERF(B1) ERFB2=ERF(B2) WRITE(NO,5)B1,ERFB1,B2,ERFB2 5 FORMAT(//,' ERF FROM HASTINGS APPROXIMATION',/, + ' B1 = ',F6.2,' ERF(B1) = ',F10.7,/, + ' B2 = ',F6.2,' ERF(B2) = ',F10.7,/) C... C... NEXT TIME 3 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