St с

Sf с SOL?K ГОЯ TIE COIKICTIOir, ОТ

и CALL SOLVC(N,AJ,JPVT) 2 С

(3 С ХЯСКЕМСМТ Т *4 с

S DO 11 I 1,11

6 DT - -R(i) 7 11 T(I) - Td) ♦ DT

€й GOTO $

TO с GSMCRATE OUTPUT

71 с

72 12 WRITE($,13)IT. HST

73 13 rORMAK , AFTER M3,* ITBRATIOMS TK RXS RESIDVAb X$*

74 WRITE(S,14){T(I),I 1,II)

75 14 rORMATC T 7Г10.7) 7 15 CONTIIIUE

77 STOP

ti WO

Рис. 6.5. Распечатка программы NEWTON.

квадратичная невязка оказывается меньше выбранного критерия или если превышено максимально допустимое число итераций. Подпрограмма JACOB (рис. 6.7) предназначена для расчета элементов якобиана матрицы (6.11). Линейная система (6.10) решается в два этапа. Сначала якобиан факторизуется,

1 SUBROUTINE RESID(N,T,R)

3 С EVALUATES RESIDUALS REQUIRED BY NEVTONS METHOD

5 DIMENSION T(50),R(50),T4(50)

7 DO 1 J = 1,Н

8 DUM = T{J)

9 1 T4(J) DU}1*DUM*DUM*DUM

10 С 0.05348

11 DUM = T4(2) + C*T(2)

12 DAM = T4(3) + C*T(3)

13 R(l) T4{1) + 0.06823*T(1) - DUM - 0.01509

14 R(2) T4(l) + C*T{1) - 2.*DUM + DAM

15 R(3) = DUM - 2.05*T4{3) - 0.2534*T(3) + 0.06698

16 RETURN

17 END

Рис. 6.6. Распечатка подпрограммы RESID.

преобразуясь к форме LU с помощью подпрограммы FACT, а затем преобразованная система решается согласно подпрограмме SOLVE. Подпрограммы FACT и SOLVE описываются в п. 6.2.1.

Форма представления решений уравнений (6.9) показана на рис. 6.8. Как свидетельствует решение, температуры поверх-

ности поглотителя и двух стеклянных крышек равны соответственно 415 К (142°С), 379 К (106°С) и 334 К (6ГС). Для среднеквадратичных значений ошибки и невязки характерна

2 3 4 5 6 7 8 Э 10 11 12 13 14 15 16 17 18 19 20 21 22

SUBROUTINE JACOB(N,T,AJ)

EVALUATES JACOBIAN REQUIRED BY NEWTONS METHOD

DIMENSION T(50),AJ(50,50)Д3(50)

DO 2 I 1,K DO 1 J 1,M

1 AJ(I,J) 0. DUM T<I)

T3(I) 4.*DUK*DUK*DUM

2 COMTIHUE

С 0,05648

AJ{1,1) T3{1) + 0.06823 AJ(1,2) - T3(2) - С T3(l) + С

- 2.*T3(2) - 0Л16Э6 T3(3) + С T3(2) + С

- 2.05*T3(3) - 0.2534

AJ(2,1) AJ(2,2> AJ(2,3) AJ(3,2) AJ(3,3) RETURN END

Рис. 6.7. Распечатка подпрограммы JACOB.

примерно квадратичная сходимость по мере приближения к сходящемуся решению, соответствующая формуле (6.7). Этот

NEWTONS METHOD FOR N = INIT SOLN = .30 .30

3 ITMX= .30

10 EPS= .lOOE-04

RMS RH-RNS RH RMS RH RMS RH-RM8 RH

.7023D-02 RMS T

.9400D-02 RMS T

.1317D-02 RMS T

.2372D-04 RMS T

.5353D~08 RMS T

.8307D-01 .2885D-01 .3148D-02 .4647D-04 .4867D-07

-.1216D-01 .1344D-01 .1738D-02 .2928D-04 .6588D~08

.6971D-09 -.1007D-05 .7722D-02 .4993D-02 .1449D-02 .2885D-03 .2881D-04 .4479D-0& .5570D-08 .3399D-0S

AFTER 4 ITERATIONS THE RMS RESIDUAL IS T - .4151283 .3794904 .3335792

.48667E-07

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

пример показывает, насколько эффективным является метод Ньютона при решении небольших систем сильно нелинейных алгебраических уравнений.

/+l.fe- /-l.fe

(6.14)

LxxUj, к =--.

\Ау AyJ

I -(-L- LV

уу - \\у hy yЧ

Решение ищется в области, изображенной на рис. 6.9.

Граничные значения взяты из точного решения [Fletcher, 1983]

= ai -f Ogjc -f -f аху -f a ((-o) e -(-o) cos) (Xy), (6.16)

тогда как постоянные ai, аг, аз, 4, as, Х и Хо подбираются так, чтобы разнообразить поведение точного решения.

Если реализовать представление уравнений (6.13) во всех внутренних узлах, то получится 2{NX - 2) {NY -2) уравнений, которые должны быть разрешены относительно Uj, k и у/, k для всех внутренних узлов. Процесс решения с использованием

6.1,3. NEWTBU: двумерные нестационарные уравнения Бюргерса

Нестационарная форма вышеупомянутых уравнений более детально рассмотрена в § 10.4. Здесь же нам необходимо решить стационарные двумерные уравнения Бюргерса

. дй , - да \ (дЧ , да \

с граничными условиями Дирихле для и и v, применяя с этой целью метод Ньютона после дискретизации с помощью трехточечных центрированных разностных формул. Уравнение (6.12) можно записать в следующей дискретной форме:

Щ, к = /, кхЩ, к + i, куЩ, к - -{LxxU,; к + L Uj к) = О,

J (6.13)

/. * == klxV/, к + V,; kLyVj, к - r7 {LxxVj, к + LyyV, ft) = О,