чительно больший радиус сходимости, чем при применении обычного метода Ньютона. Характерный алгоритм, известный как формула BFGS [Fletcher, 1980], позволяет строить следующее приближение у(*+> с помощью формулы (6.21). В результате имеем

y()=R(+0 -

+ р() (1 + р(/)у(л)гн()у()) р(Л)р(Л) г. (6.22)

Несмотря на то что матрица J является, как правило, разреженной, структура выражений (6.22) свидетельствует о том, ято матрица Н оказывается плотной. При k = 0 матрица полагается равной единичной матрице I. Эффективность квазиньютоновских методов зачастую зависит от того, обладает ли матрица J некоторыми специальными свойствами, такими, как положительная определенность [Jen-tiings, 1977а]. Следовательно, применимость этих методов к задачам вычислительной гидроаэродинамики должна исследоваться для каждого случая отдельно. Значительная часть литературы по квазиньютоновским методам связывает их приме- ение с безусловной минимизацией (см., например, [Powell, 1976; Fletcher, 1980]). Обсуждение квазиньютоновских методов в применении к специальным задачам можно найти в работах [Broyden, 1965; Ortega, Rheinboldt, 1970; Shanno, 1983]. Характерный вариант применения квазиньютоновского метода к задаче гидроаэродинамики дается в работе [Engelman et al., 3981].

§ 6.2. Прямые методы для линейных систем

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

AV = В (6.23)

-относительно компонентов вектора V, где все элементы матрицы А и вектора В известны. Метод исключения по Гауссу [Dahlquist, Bjorck, 1974] является предпочтительным средством решения уравнений (6.23), однако эффективность реализации

ЭТОГО метода зависит от структуры матрицы А. Заслуживают выделения три категории подобного рода матриц:

(a) Матрица А содержит лишь несколько или не содержит совсем нулевых элементов. Тогда говорят, что матрица А полная или плотная.

(b) Матрица А содержит много нулевых элементов. Такую матрицу называют разреженной.

(c) Матрица А содержит много нулевых элементов, причем ее ненулевые элементы группируются вблизи главной диагонали. Тогда матрицу А называют разреженной и ленточной.

Компьютерные программы, основанные на предположении о том, что матрица А является полной, оказываются более универсальными, но и более дорогостоящими в вычислительном смысле. Такие программы могут использоваться в применении к разреженным и/или ленточным матрицам, однако для последних категорий существуют специальные и притом более экономичные процедуры. Какой бы ни была структура матрицы А, метод исключения по Гауссу должен применяться в два этапа. Сначала матрица А факторизуется, превращаясь в нижнетреугольную матрицу L и верхнетреугольную матрицу U в качестве двух сомножителей. После этого факторизованная форма уравнения (6.23) решается как

UV = L B. (6.24>

Вследствие структуры матрицы U процесс решения (6.24) сводится к повторным подстановкам, осуществляемым после формирования комплекса L-B.

6.2.1. FACT/SOLVE: решение систем с плотными матрицами

Если матрица А плотная, то все ее элементы участвуют в образовании матриц L и U. Подпрограмма FACT (рис. 6.15) реализует требуемую факторизацию по формуле LU = А. Последующая серия подстановок согласно (6.24) осуществляется с помощью подпрограммы SOLVE (рис. 6.16).

На каждом шаге факторизации множество элементов иа текущего ряда k вычитается из элементов последующего ряда / с целью исключения и. При этом вводится множитель

л/лл, ky а величина а, k называется eedyuiUM элементом. Если ведущий элемент равен нулю или очень мал, то исключение по Гауссу либо не проходит вообще, либо дает очень неточные результаты. Во избежание этого вводится частичный отбор ведущего элемента. Для исключения Vh проводится поиск в k-u столбце в интервале i = *+l, Л/, ставящий целью найти

1 SUBROUTINE FACT{N,A,JFVT)

3 С FACTORISES А INTO PERMUTED L.U SO THAT PERM*A L*U

4 С JPVT(К) GIVES INDEX OF KTH PIVOT ROW

5 С SETS J?VT(N) = -1 IF ZERO PIVOT OCCURS

7 DIMENSION A(50,50),JPVT(50)

8 NM = N - 1

10 С GAUSS ELIMINATION WITH PARTIAL PIVOTING

11 С

12 DO 5 К = 1,NM

13 KP = К + 1

14 С

15 С SELECT PIVOT

16 С

17 L = К

18 DO 1 I = KP,N

19 IF(ABS(A(I,K)) .GT. ABS(A(L,K)))L = I

20 1 CONTINUE

21 JPVT(K) = L

22 S = A(L,K)

23 A{L,K) = A(K,K)

24 А(КД) = S

25 С

26 С CHECK FOR ZERO PIVOT

27 С

28 IF(ABS(S) .LT. 1.0E-15)GOTO 6

29 С

30 С CALCULATE MULTIPLIERS

31 С

32 DO 2 I = KP,N

33 A(I,K) = -A(I,K)/S

34 2 CONTINUE

35 С

J6 С INTERCHANGE AND ELIMINATE BY COLUMNS

37 С

38 DO 4 J = KP,N

39 S = A(L,J)

40 A(L,J) = A(K,J)

41 A(K,J) = S

42 IF(ABS(S) .LT. 1.0E-15)GOTO 4

43 DO 3 I = KP,ir

44 A(I,J) = A{I,J) + A(I,K)*S

45 3 CONTINUE

46 4 CONTINUE

47 5 CONTINUE

48 RETURN

ЛЭ 6 JPVT(N) = -1

50 RETURN

51 END

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