метода Ньютона может быть символически представлен формулой

j( )Aq( +i) = --RH (6.17)

где составляющими вектора Aq являются Агг/, Au/, а вектор R имеет компоненты Ruj и и Rvj, и, оцениваемые в каждом из

2/=0

u=u(x,0), v=v[xfi)

Рис. 6.9. Вычислительная область для решения уравнений (6.12).

внутренних узлов. Верхний индекс {п) используется для обозначения номера итерации. Якобиан J содержит такие члены:

dRUf ,/duj ,, = ~ O.Suj JAx - l/(Re А:),

dRUf JdUf, = (2/Re) (1/Ajc2 + 1/Ду2) (0.5/Д;,) ( (5/? д ,/(5 y , = 0.5a;, JAx - l/(Re Ад:),

dRUfJdVfk-i = -0.5vjJAy\/{ReAy\

dRuj JdVfk, = O.SVfJAy ~ l/(Re Ay\ dRuj k/dVf k = {0.5/Ay) (wy, .i r/y, ;.,).

Аналогичные члены можно записать для Rvj\

Вышеприведенная формулировку задачи и ее решение реализуются в программе NEWTBU. Распечатки этой программы, а также подпрограмм EXBUR для построения точного решения, RESBU для расчета невязок согласно (6.13) и JACBU для вычислений по формулам (6.18) приводятся на рисунках 6.10-6.13. В табл. 6.2 расшифровывается смысл параметров, используемых в программе NEWTBU.

3 С NEWTBU APPLIES NEVTONS METHOD TO SOLVE BURGERS EQUATIONS

4 С FOR U(X,y) AND V(X,y).

5 С EXBUR EVALUATES THE EXACT SOLUTION OF BURGERS EQUATIONS

6 С RESBU EVALUATES THE RESIDUALS

7 С JACBU EVALUATES THE JACOBIAN

8 С FACT FACTORISES THE JACOBIAN INTO L.U

9 С SOLVE SOLVES THE LINEAR SYSTEM FOR DT

10 С IRD CONTROLS READ/WRITE OF STARTING SOLUTION

11 С IPN IS THE INCREMENT FOR PRINTING ITERATIVE RMS

12 С

13 REAL*8 SUM,DSQRT,R,U,V

14 DIMENSION и(21,21),V(21,21),UE(21,21),VE(21,21),A(5),R(50>

15 1,AJ(50,50),JPVT(50),RD(50)

16 COMMON DX,DY,RE,NX,NY,DTI

17 С

18 OPEN(1, FILE= NEWTBU.DAT)

19 OPEN(6, FILE* NEWTBU.OUT)

20 OPEN(7, FILE= NEWTBU.STO)

21 READ(1,1)NX,NY,ITMX,IRD,IPN,EPS,RE,0M,DT

22 READ(1,2)(A(J),J l,5),AL

23 1 FORMAT(5I5,4E10.3)

24 2 FORMAT(5F8.2,F5.2)

25 С

26 NXP = NX - 1

27 NYP = NY - 1

28 N (NX -2)*(NY - 2)*2

29 WRITE(6, 3)N,ITMX,IRD,IPN

30 3 FORMATC NEWTONS METHOD FOR N ,13,5X, ITMX M3,

31 1 IRD=,I2, IPN M2)

32 WRITE(6,4)DT,EPS,RE,OM

33 4 FORMATC DT ,E12.5/ EPS \E10.3 RE \E10.3,

34 1 OM ,E10.3)

35 WRITE{65)(A(J),J=l,5) ,AL

36 5 FORMATC A =.5F8.2, AL ,F5.2, )

37 С

38 С CALCULATE EXACT SOLUTION

39 С

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

41 С

42 С GENERATE INITIAL SOLUTION

43 С

44 IF (IRD .LT. 2) GOTO 8

45 С

46 С IRD = 2 OBTAIN STARTING SOLUTION FROM NEWTBU.STO (UNIT 7>

47 С

48 DO 7 К * 1,NY

49 READ(7,6)(U(K,J),J 1,NX)

50 READ(7,6)(V(K.J),J=1,NX)

51 6 FORMAT(10F8y5)

52 7 CONTINUE

53 GOTO 11

54 <:

Рис. 6.10. Распечатка программы NEWTBU (начало). 16 К. Флетчер, т. 1

55 8 ЛЮ 10 J = 1,NX

56 DO 9 К = 1,NY

57 U(K,J) UE(K,J)

58 V(K,J) = VE{K,J)

59 9 CONTINUE €0 10 CONTINUE €1 с

62 И DO 12 К = 1,KY

€3 12 УК1ТЕ(6ДЗ) (UE(K,J) J 1,NX)

4 13 FORHATC UE ,7ri0.4)

65 DO 14 К 1,NY

€6 14 VRITE(6,15){VE(K,J),J 1,HX)

67 15 FORHATC VE ,7F10.4)

9 AM N

70 IT = 0

71 ITP 0

72 OMH = OM

73 DTI l./DT

74 16 CONTINUE

75 с

76 с CALCULATE RESIDUALS

77 с

78 CALL RESBU(U,V,R)

79 с

0 SUM 0.

1 DO 17 I 1,N

2 RD(I) Rd)

83 17 SUM SUM + R(I)*R(I)

S4 18 RMS ж DSQRT(SUM/AN)

85 IF(IT CE. ITP)VRITE(6,19)IT,RMS,(R(J),J-1.3)

M 19 FORMATC IT M3. RMS * ,D11.4, * R-43D11.4)

87 с

$S с TEST FOR CONVERGENCE OR MAXIMUM NUMBER

:89 С OF ITERATIONS EXCEEDED Э0 С

91 IF(RMS .LT. EPS)GOTO 23

S2 IF(IT -GE. ITP)ITP ITP + IPM

93 IT IT + 1

4 IF(IT .EQ. ITMX)GOTO 23

95 IF(RMS .GT. 1.0E+03)GOTO 30

96 С

97 С CALCULATE JACOBIAN

98 С

99 CALL JACBU(U,V,AJ)

100 С

101 С FACTORISE THE JACOBIAN INTO ЬЛ

102 С

103 CALL FACT(N,AJ,JPVT):

104 С

105 IF(JPVT(N) .EQi -1)VRITE(6,20)

106 20 FORMATC ZERO PIVOT DETECTED)

107 IF(JPVT(N) .EQ. -1)GOTO 30

108 С

109 С SOLVE FOR THE CORRECTION, DT

110 С

111 CALL SOLVE(N,AJ,JPVT,RD)

112 С

113 С INCREMENT U AND Y

114 С

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