DIFF算法语句 fortran

怪咖先生

C SOLVES 1D TRANSTENT HEAT CONDUCTION EQUATION USING
C FTCS SCHEME
C
DIMENSION TN(41),DUM(41),TD(41),X(41),TE(41)
C
C INPUT DATA;
C
C JMAX = THE NUMBER OF POINTS ALONG THE ROD
C MAXEX = THE NUMBER OF TERMS IN THE EXACT SOLUTION
C NMAX = THE MAXIMUM NUMBER OF TIME STEPS
C ALPH = THE THERMAL DIFFUSIVITY
C S = ALPH*DELT/DELX/DELX
C TMAX = THE MAXIMUM TIME
C
OPEN(1, FILE='DIFF.DAT')
OPEN(6, FILE='DIFF.OUT')
READ(1,1) JMAX, MAXEX ,NMAX,ALPH,S,TMAX
1 FORMAT(3I5,E10,3,F5,2,F5,0)
PI = 3.1415927
C
C TD = DIMENSIONAL TEMPERATURE
C TN = NONDIMENSIONAL TEMPERTURE
C
JMAP = JMAX-1
AJM = JMAP
DELX = 1./AJM
DELT = DELX*DELX*S/ALPH
WRITE(6,2)JMAX,MAXEX,NMAX,TMAX
2 FORMAT(' JMAX= ', I5, ' MAXEX= ', I5, 'NMAX= ', I5, ' TMAX= ', F8.2)
WAITE ( 6,3) S, ALPH, DELT, DELX
3 FORMA(' S= ', F5.3, ' ALPH= ', E10.3, ' DELT= ' , E10.3,
1 ' DELX= ', E10.3，//)
WEITE(6,4)S
4 FORMAT(' FTCS(EXPLICIT) SCHEME , ' , 5X, ' S= ', F5.3, //)
C
C SET INITIAL CONDITIONS
C
DO 5 J = 1 , JMAP
5 TM(J) = 0.
N = 0
T = 0.
SJ = 1.0-2.0*S
C
C SET BOUNDARY CONDITIONS
C EACH TIME STEP STARTS AT STATEMENT 6
C
6 TN(1) = 1.
TN(JMAX) = 1.
IF(T .LT. 0.01)TN(1) =0.5
IF(T .LT. 0.01)TN(JMAX) =0.5
TD(1) = 100.*TN(1)
TD(JMAX) = 100.*TN(JMAX)
C
C COMPUTE F.D SOLUTION
C
DO 7 J = 2, JMAP
DUM(J) = SJ*TN(J) + S*(TN(J-1) + TN(J+1))
7 CONTINUE
DO 8 J = 2, JMAP
8 TN(J) = DUM(J)
C
DO 9 J = 2, JMAP
9 TD(J) =100.*TN(J)
T = T+ DELT
WRITE(6,10)T,(TD(J), J=1, JMAX)
10 FORMAT(' T= ', F5.0, ' TD= ', 11F6.2)
C
C IF MAXIMUM TIME OR MAXIMUM NUMBER OF TIME-STEPS EXCEEDED
C EXIT FROM LOOP
C
IF(N .GE. NMAX)GOTO 11
IF(T .LT. TMAX)GOTO 6
C
C OBTAIN EXACT SOLUTION AND COMPARE
C
11 SUM = 0.
DO 13 J = 1,JMAX
AJ =J - 1
X(J) = DELX*AJ
TE(J) = 100.
DO 12 M = 1,MAXEX
AM = M
DAM = (2.*AM - 1.)
DXM = DAM*PI*X(J)
DTM = -ALPH*DAM*DAM*PI*PI*T
C
C LIMIT THE ARGUMENT SIZE OF EXP(DTM)
C
IF(DTM .LT. -87.)DTM = -87.0
12 TE(J) = TE(J) - 400./DAM/PI*SIN(DXM)*EXP(DTM)
SUM = SUM + (TE(J) -TD(J))*2
13 CONTINUE
WRITE(6,14)T, (TE(J), J=1, JMAX)
14 FORMAT(/,' T= ', F5.0,' TE= ', 11F6.2,//)
C
C RMS IS THE RMS ERROR
C
AVS = SUM/(1. + AJM)
RMS = SQRT(AVS)
WRITE(6,15)RMS
15 FORMAT(' RMS DIF = ', E11.4,//)
STOP
END
运行结果：




如有错误 请大家给指出来

转自：http://blog.sina.com.cn/s/blog_a319f5ff0101iho9.html