2 с - /

3 С DIFIM.FOR SOLVES THE DIFFUSION EQUATION USINil;

4 С VARIOUS IMPLICIT SCHEMES /

6 REAL*8 SUM,AVS,RMS,DSQRT,DMP

7 DIMENSION DUM(65),X(41),TE(41),A(5,65),D(6,5)

8 DIMENSION ELX(3),EMX(3),T0L(41),TN(41),TD(41)

10 OPEN(1,FILE- DIFIM.DAT)

11 OPEN(6,FILE- DIFIM.OUT)

12 READ(1,1)ME,IPR,JMAX, MAXEX,NMAX,ALPH,S,TMAX,TST,GAM

13 1 FORMAT(5I5,E10.3,F5.3,3F5.2)

14 О

15 PI - 3.1415927

16 JMAP - JMAX - 1

17 AJM - JMAP

18 DELX - l./AJM

19 DELT - DELX*DELX*S/ALPH

20 С

21 IF(ME .EQ. 1)WRITE(6,4)

22 IF(ME .EQ. 2)VRITE(6,5)

23 IF(ME .EQ. 3)VRITE(6,6)

24 IF(ME .EQ. 4)WRITE{6,7)

25 IFCME .EQ. 5)WRITE(6,8)

26 WRITE(6,2)JMAX,MAXEX,NMAX,TMAX,TST,GAM

27 2 FORMATC JMAX-M5,* MAXEX-M5. NMAX-*,15, 2S 1* TMAX-*,F8.2,* TST--,F5.2,* GAM-*,F5.2)

29 WRITE(6,3)S,ALPH;DELT,DELX

30 3 FORMAT(*.S ,F5.3,* ALPH - .ЕЮ.З,* DELT - ЕЮ.З,

31 1* DELX - *,E10.3, )

32 4 FORMATC DIFFUSION EQUATION: FDM-2ND*, )

33 5 FORMATC DIFFUSION EQUATION: FEM-2ND*, )

34 6 FORMAT(* DIFFUSION EQUATION: FDM-4TH*, )

35 7 FORMAT(* DIFFUSION EQUATION: FEM-4TH*, )

36 8 FORMATC DIFFUSION EQUATION: COMP*, )

37 С

38 С IMPLICIT PARAMETERS

39 С

40 BET - 0.5 + GAM

41 IFCME .EQ. 3)BET = BET - 1./12./S

42 IFCME .EQ. 4)BET - BET + 1./12./S

43 IFCME .EQ. 1 .OR. ME .EQ. 3)EMX(1) = 0.

44 IFCME .EQ. 2 .OR. ME .EQ. 4)EMXC1) - 1./6.

45 IFCME .EQ. 5)EMXC1) - 1./12.

46 БМХС2) - 1. 2.*EMXC1)

47 EMXC3) - EMXCl)

48 AD - EMXCl) * CI. +4aM) - BET*S

49 BD = EMXC2)*C1. + GAM) + 2.*BET*S

50 CD - AD

51 ELX(l) = 1.

52 ELX(2) = -2.

53 ELX(3) = 1.

54 JMAF - JMAX - 2

55 С

56 С OBTAIN INITIAL CONDITIONS FROM THE EXACT SOLUTION

57 С

5$ \Т TST - 2.*DELT

59 00 10 I 1,2

60 T = T + DELT

61 С \

62 CMiL EXTRA(JNAX,MAXEX,DELX,PI,ALPH,T,TE)

63 С Y

64 DO J 2,JMAP

65 IFd .EQ. l)TOL(J) = TE(J)/100.

66 IF(H .EQ. 2)TM(J) = TE(J)/100.

67 9 CONTINUE

68 10 CONTINUE

69 С

70 С SET BOUNDARY CONDITIONS

71 С

72 TOL(l) = 1.

73 TOL(JMAX) 1.

74 TN(1) 1.

75 TH(JMAX) 1.

76 TD(1) 100.*TN(1)

77 TD(JMAX) 100.*TN(JMAX)

78 N = 0

79 С

80 С EACH TIME STEP STARTS AT STATEMENT 11

81 С

82 11 N N 1

83 С

84 С SET UP THE TRIDIAGONAL SYSTEM OF EQUATIONS

85 С

86 DO 14 J 2,JMAP 7 JM J - 1

88 ir(N .GT. DGOTO 12

89 A(1,JM) 0.

90 A(2,JH) AD

91 A(3,JM) BD

92 A(4,JM) CD

93 A(5,JM) 0.

94 12 D(JM) 0.

95 DO 13 К 1,3

96 KJ J - 2 + К

97 D(JM) s D(JM) 4- EMX(K)*((1. * 2.*GAH)*TN(ICJ) - GAM*TOL(ICJ}>

98 D(JM) D(JM) + S*ELX(K)*(1.-BET)*TN(KJ)

99 13 CONTINUE

100 14 CONTINUE

101 D(l) D(l) - A(2,1)*TH(1)

102 D(JMAF) = D(JMAF) - A(4,JMAF)*TN(JMAX)

103 С

104 С SOLVE BANDED SYSTEM OF EQUATIONS

105 С

106 IF(N .EQ. DCALL BANFAC (A, JMAF, 1)

107 С

108 CALL BANSOL(D,DUM,A,JMAF,1)

109 С

110 DO 15 J - 2,JMAP

111 TOL(J) - TN(J)

112 15 TN(J) = DUM(J-l)

113 С

114 DO 16 J s 2,JMAP

115 16 TD(J) 100.*TN(J)

116 T T + DELT

117 IFdPR .EQ. 1)WRITE(6,17)T,(TD(J),J=1,JMAX) 11$ 17 FORMATC T \Г5.2, TD M1F6.2)

Рис. 7.7 (продолжение).

119 С

120 С XF MAXIMUM TIME OR MAXIMUM NUMBER OF TIME-STi S EXCEEDEi

121 С EXIT FROM LOOP

122 С

123 IF(N .GE. NMAX)GOTO 18

124 IF(T .LT. .TMAX)GOTO 11

125 С

126 С OBTAIN EXACT SOLUTION AND COMPARE

127 С

128 18 CALL EXTRA(JMAX,MAXEX,DELX,PI,ALPH,T,ТЕ)

129 С

130 SUM 0.

131 DO 19 J 1,JMAX

132 DMP - TE(J) - TD(J)

133 SUM = SUM + DMP*DMP

134 19 CONTINUE

135 IFdPR .NE. 1)WRITE(6,17)T,(TD(J),J 1,JMAX)

136 WRITE(6,20)T,(ТЕ(J),J=l,JMAX)

137 20 FORMATC/, T ,F5.2,* TE M1F6.2, )

138 С

139 С RMS IS THE RMS ERROR

140 С

141 AVS > SUM/(1. -I- AJM)

142 RMS - DSQRT(AVS)

143 WRITE(6,21)RMS

144 21 FORMATC RMS DIF ,D11.4, )

145 STOP 14$ END

Рис. 7.7. (окончание). Ц

начальном условии, соответствующем точному решению (3.42) при t = 4.5. Граничные условия Дирихле были наложены на Т при х = 0 и 1.0. Приближенная скорость сходимости была получена такой, как в табл. 7.3.

DIFFUSION EQUATION: ЕЕМ-4ТЯ

JMAX- 11 МАХЕХ- 10 NMAX- 500 ТМАХ 12.00 TST- 4.50 GAM- .00 S- 1.000 ALPH - .lOOE-01 DELT - .lOOE+01 DELX - .lOOE+00

T- 12.50 TD-100.00 88.55 78.22 70.02 64.76 62.94 64.76 70.02 78.22 88.55100.00 T 12.50 TE-100.00 88.54 78.21 70.00 64.74 62.92 64.74 70.00 78.21 88.54100.00

RMS DIF - .1522D-01

Рис. 7.8. Типичная выдача результатов по программе DIFIM.

Все результаты, приводимые в табл. 7.6, были получены при 5=1.0. Методы, номинально имеющие второй порядок (случаи 1 и 2), обнаруживают приблизительно второй порядок схо-