необходимость в проведении новой факторизации блочно-трех-диагональной системы на каждой итерации п. Эффективное решение блочно-трехдиагональной системы строится с помощью блочного алгоритма Томаса (п. 6.2.5).

Анализ устойчивости по Нейману применяется к схеме (10.70), (10.71) после временного замораживания значений А и т. д. Этот анализ показывает, что алгоритм схемы расщепления (10.70), (10.71) является безусловно устойчивым. На практике использование чрезмерно большого шага по времени может привести к нелинейной устойчивости.

Схема расщепления имеет ошибку аппроксимации порядка О (А/, Ал:2, Ду2) Поскольку нас интересует только решение для стационарного состояния, можно рассмотреть такие модификации вида левых частей уравнений (10.70) и (10.71), которые позволили бы достичь стационарного состояния при уменьшенном объеме вычислений. В пределе стационарного состояния RHS = О, так что никакая модификация левой части не может оказать влияния на решение для стационарного состояния, даже если она и повлияет на нестационарное решение. Предпочтительная стратегия сводится к такому выбору формы левых частей, чтобы на каждой итерации построение решения было наиболее эффективным, однако чтобы можно было избежать существенного увеличения числа итераций, необходимых для достижения стационарного состояния.

Как отмечалось в п. 6.2.5, т скалярных трехдиагональных систем могут быть решены более эффективно, чем одна блоч-но-трехдиагональная система уравнений размерностью ту, т. Уравнения (10.70) и (10.71) могут каждое в отдельности рассматриваться так, как если бы это были две последовательно решаемые трехдиагональные системы. Для этого в формулах (10.69) следует вычеркнуть члены, не попадающие на диагональ. Дополнительная экономия может быть получена за счет замены формул (10.69) на

-с:). -()(::)- о -)

Это означает, что одинаковые левые части могут быть использованы для каждой из скалярных компонент уравнений (10.70) и (10.71). Иначе говоря, левая часть должна подвергаться факторизации (подпрограмма BANFAC) лишь один раз для двух решений (подпрограмма BANSOL) с различными правыми частями ARHS? и ARHS . Формулы (10.72) вместе с уравнениями (10.70) и (10.71) являются составной частью программы TWBURG (рис. 10.13).

3 С TVBURG APPLIES APPROXIMATION FACTORISATION TO SOLVE

4 С THE STEADY 2-D BURGERS EQUATIONS FOR U(X,Y) AND V(X,Y)

6 DIMENSION RU(21,21),RV(21,21),U(21,21),V(21,21),UE(21,21),

7 1VE(21,21),Ъ{Ь,еЪ),RRU(65),DDU(65),RRV(65),DDV(65),CX(3),CY(3),

8 2CXQ(4),CYQ(4),EMX(3),EMY{3),DMX(3,3),DMY(3,3)Д(5)

9 COMMON DX,Dy,RE,NX,NY,CX,CXQ,CY,CYQ,DMX,DMY,EMX,EMY,RU,RV,U,V

10 OPEN(1.FILE= TVBURG.DAT)

11 OPEN(б, FILE= TWBURG.OUT)

12 READd ,1)NX,NY,ME, ITMAX, IPRrBET,QX,QY

13 READ(1,2)DTIM,EPS,RE,EM1,EM2

14 READ(1,3)(A(J),J=1,5),AL

15 1 FORMAT(515,3F5.2)

16 2 FORMAT(5E10.3)

17 3 FORMAT(2E10.3,4F5.2)

18 С

19 NXP = NX - 1

20 NXPP = NXP - 1

21 NYP = NY - 1

22 NYPP = NYP - 1

23 AN = NXPP*NYPP

24 С

25 CALL EXBUR(UE,VE,A,AL)

26 С

27 CX(l) = -0.5/DX

28 CX(2) = 0.

29 CX(3) = 0.5/DX

30 CY(1) = -0.5/DY

31 CY(2) = 0.

32 CY(3) = 0.5/DY

33 CXQ(l) = QX/DX/3.

34 CXQ(2) = -3.*CXQ(1)

35 CXQ(3) = -CXQ(2)

36 CXQ(4) = -CXQ(l)

37 CYQ(l) = QY/DY/3.

38 CYQ{2) = -3.*CYQ(1)

39 CYQ(3) = -CYQ(2)

40 CYQ(4) = -CYQ(l)

41 CCX = l./DX/DX/RE

42 CCY = l./DY/DY/RE

43 IF(ME .NE. 2)QX = 0.

44 IF(ME .NE. 2)QY = 0.

45 EMXd) = 0.

46 EMY(l) = 0.

47 IF(ME .EQ. 3)EMX(1) = EMI

48 IFIME .EQ. 3)EMY(1) = EM2

49 EMX(2) = 1. - 2.*EMX(1)

50 EMX(3) = EMXd)

51 EMY(2) = 1. - 2.*EMY(1)

52 EMY(3) = EMYd)

53 DO 4 I = 1,3

54 DMX(I,1) = CCX*EMY(I)

55 DMX(I,2) = -2.*DMX(I,1)

56 DMX(I,3) = DMX(I,1)

57 DMYd,I) = CCY*EMX(I)

58 DMY(2,I) = -2.*DMY(1,I)

59 4 DMY(3,I) = DMY(1,I)

60 1F(ME .EQ. l)VRITE(6,5)NX,My,HE,ITHAX

61 IF(ME .EQ. 2)VRITE(6,6)NX,NY,HE,IT1!AX

62 IF(HE .EQ. 3)VRITE(6,7)NX,Hy,HE,ITKAX

63 5 FORMATC 2-D BURGERS EQUATION: AF-FDM, NX,KY ,213, HK

64 112/ ITMAX M3)

65 6 FORMATC 2-D BURGERS EQUATION: AF-4PU, HX,NY ,2I3/ ME. 4

66 112/ ITMAX M3)

67 7 FORMATC 2-D BURGERS EQUATION: AF-MO, NX,HY ,2I3/ KB

68 112/ ITMAX M3)

69 VRITE(6,8)BET,DTIM,RE,QX,Qy,EMX(l),EHy(l)

70 8 FORMATC BETA \F5.2/ DTIM \F5.3/ КЕ \Г5Л/ QX M4.a

71 1/ 0У \Г4.2/ ЕМХ \Г5.3/ ЕМТ *,Г5.3)

72 WRITE(6,9)(A(J),J l,5),AL,DX,Dy

73 9 FORMATC A \2E10.3,3F5,2/ AL ,F5.2/ DX,Dy * ,2Г8.§,/)

74 С

75 С GENERATE INITIAL SOLUTION

76 С

77 DO 10 J 1,NX

78 V(1,J) VE(1,J)

79 V(NY,J) VE(Ny.J)

80 Uil,J) UE(1,J)

81 10 U(Ny,J) UE(Ny,J)

82 DO 12 К Ш 2,NyP

83 U(K,1) UE(K,1)

84 U(K,NX) UE(K,NX)

85 V(K,1) VE(K,1)

86 y(K,NX) VE(K,NX)

87 DO 11 J 2,NXP

88 AJ J - 1

89 U(K,J) U(K,1) + 0.5*AJ*DX*(U(K,NX)-U(K,1)>

90 11 V(K,J) Y(K,1) + 0.5*AJ*DX*(y(K,NX)-V(K,U)

91 12 CONTINUE

92 IFdPR .LE. DGOTO 17

93 DO 13 К l,Ny

94 13 VRITE(6,14)(UE(K,J),J 1,NX)

95 14 FORMAT(4H UE .7F10.4)

96 DO 15 К l,Ny

97 15 WRITE(6 Д6) (УЕ(К, J), J>1,NX)

98 16 F0RNAT(4H VE ,7F10.4)

99 17 PO 20 J 1,NX

100 DO 18 К l,Ny

101 RU(K,J) 0.

102 18 RV(K,J) 0.

103 DO 19 К 1.5

104 19 B(K,J) 0.

105 20 CONTINUE

106 ITER 0

107 с

108 21 DBT = BET*DTIM

109 CXA 0.5*DBT/DX

110 СУА 0.5*DBT/Dy

111 CCXA = DBT*CCX

112 ССУА DBT*CCy

113 с

114 CALL RHSBU(RMSU,RMSy,ME,DTIM)

115 с

116 IF(RMSU .LT. EPS)GOTO 34

117 IF(RMSU .GT. i:0E+04)GOTO 34

118 с

119 с TRIDIAGONAL SYSTEMS IN THE X-DIRECTIOK

120 с

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