В качестве альтернативной и более эффективной стратегии по сравнению с применением формулы (6.19) можно указать на сочетание метода Ньютона с псевдонестационарной форму-

2 SUBROUTINE RESBU(U,V,R)

4 с EVALUATES RESIDUALS OF 2D STEADY BURGERS EQUATIONS

6 REAL*8 R,U,V,CX,CY,CCX,CCY,DUH

7 DIMENSION U(21,21),V(21,21),R(50)

8 COHKON DX,DY,RE,NX,NY

9 C.

10 NXP = NX - 1

11 NYP NY - 1

12 CX 0.5/DX

13 CY 0.5/DY

14 CCX l./RE/DX/DX

15 CCY l./RE/DY/DY

16 NCT 1

17 С

1$ DO 2 J * 2,NXP

19 JM J - 1

20 JP = J + 1

21 DO 1 К = 2,NYP

22 KM = К - 1

23 KP К + 1

24 DUM CX*U(K,J)*(U(K.JP)-U{K,JM)) + CY*V(K, J) * (U(KP, J)-U(KM, J)

25 R(NCT) = DUM - CCX*(U(K,JM)-2.*U(K,J)+U(K,JP)) - CCy*(U(KM,J)

26 1 - 2.*U(K,J) + U(KP,J))

27 DUM CX*U{K,J)*(V(K,JP)-V(K,JM)) + CY*V(K,J)*(V(KP,J)-V(KM,J)>

28 R(NCT+1) DUM - CCX*(V(K. JM)-2.*V(K, J)+V(K, JP)) - CCY* (V<KM, J)

29 1 - 2.*V(K,J) + V(KP,J))

30 NCT NCT + 2

31 1 CONTINUE

32 2 CONTINUE

33 С

34 RETURN 35 END

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

лировкой той же задачи (см. § 6.4). Решение данной задачи с помощью псевдонестационарного метода Ньютона излагается в п. 6.4.1.

6.L4. Квазиньютоновский метод

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

\{k+i) = v - co()H()R<), (6.20>

SUBROUTINE JACBU{U,V,AJ)

4 С EVALUATES THE JACOBIAN OF THE 2D STEADY BURGERS EQUATIONS

6 REAL*8 U,V,CX,CY,CCX,CCY

7 DIMENSION U(21,21),V(21,21),AJ(50,50)

8 COMMON DX,DY,RE,HX,NY,DTI

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

11 DO 2 J 1,N

12 DO 1 К 1,H

13 1 AJ(K,J) 0.

14 2 CONTINUE

15 С

16 NXP = NX - 1

17 NYP NY - 1

18 CX 0.5/DX

19 CY 0.5/DY

20 CXX l./RE/DX/DX

21 CYY l./RE/DY/DY

22 С

23 DO 7 J 2,NXP

24 JM = J - 1

25 JP = J + 1

26 DO 6 К = 2,NYP

27 KM К - 1

28 KP К + 1

29 MB = 2*(K-2) + 2*(NY-2)*(J-2)

30 LU = MB + 1

31 LV = MB + 2

32 ML = MB - 2*(NY - 2)

33 IF(ML .LT. 0)GOTO 3

34 AJ(LU,ML+1) = - CX*U(K,J) - CXX

35 AJ(LV,ML+2) = - CX*U(K,J) - CXX

36 3 AJ(LU,MB+1) 2.*(CXX+CYY) + CX*(U(K,JP)-U(K,JM)) + DTI

37 AJ(LU,MB+2) = CY*(U(KP,J) - U(KM,J))

38 AJ(LV,MB+1) CX*(V(K,JP) - V(K,JM))

39 AJ(LV,MB+2) 2.*(CXX+CYY) + CY*(V(KP,J) - V(KM,J)) + DTI

40 MR = MB + 2*(NY-2)

41 IF(MR .GT. N-2)G0T0 4

42 AJ(LU,MR+1) CX*U(K,J) - CXX

43 AJ(LV,MR+2) = CX*U(K,J) - CXX

44 4 MK = MB - 2

45 IF(MK .LT. 0)GOTO 5

46 AJ(LU,MK+1) = - CY*V(K,J) - CYY

47 AJ(LV,MK+2) = - CY*V(K,J) - CYY

48 5 MT = MB + 2

49 IF(MT .GT. N-2)G0T0 6

50 AJ(LU,MT+1) = Cy*V(K,J) - CYY

51 AJ(LV,MT+2) = CY*V(K,J) - CYY

52 6 CONTINUE

53 7 CONTINUE

54 С

55 . -RETURN :

56 , EIP , . , . , - .

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

где матрица Н<) представляет собой аппроксимацию величины (J*))~ систематически модифицируемую на каждом итерационном шаге, так что она приближается к (J(>)-i по мере прибли- жения к сходящемуся решению. Модификация матрицы Н(*> на каждом шаге оказывается значительно более экономичны *

NEWTONS METHOD FOR N DT .5O000E+03 EPS=

18 ITMX=400 IRD= 1

.lOOE-04 RE= .lOOE+02

IPN=20 0M= .150E+0a

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

процессом, чем факторизация матрицы J*). Уравнение (6.20> может быть переписано в форме

V(ft+i) = \(k) (oik)i{k) (6.21>

где вектор g(*> = H<*>R*) может рассматриваться как определяющий направление поиска. Скаляр (о<*) выбирают так, 4To6bL величина ijJSiV оказывалась минимальной в направлении по- иска При больших N последнее свойство обеспечивает зна-