где а/ = -(1-f0.5/?ceii), = 2, с/= - (1 - 0.5/?ceii). Ненулевые значения dt связаны с членами типа источника или, как с di и dN, с граничными условиями. Все элементы матрицы А, кроме показанных выше, равны нулю. Можно отметить изменение некоторых обозначений по сравнению с п. 6.2.1, в частности относящихся к bi.

Алгоритм Томаса для решения уравнений (6.28) состоит из двух частей (см. схему на рис. 6.17). Сначала уравнения (6.28)

X X XXX

XXX X X

Прогонка назад

Рис. 6.17. Алгоритм Томаса для решения трехдиагональной системы уравнений.

преобразуются к виду 1 с[

1 <

1 с;

1 с:

т. е. из них исключаются коэффициенты а/, тогда как коэффициенты bi нормализуются к единичным значениям. Для первого уравнения этой цели служат формулы

с[ = с,1Ь d\ = d,lb (6.29)

а для произвольно расположенного уравнения -

(6.30)

Модификация уравнений по формулам типа (6.30) проводится в процессе прямой прогонки (рис. 6.17). Второй этап алгоритма сводится к обратной подстановке (обратная прогонка на рис. 6.17), т. е. к использованию формул

N=4. = (6.31)

Алгоритм Томаса чрезвычайно экономичен; он требует для своего исполнения всего 5Л - 4 операций (умножений и делений). Однако во избежание плохой обусловленности (и, следовательно, накопления ошибок округления) необходимо, чтобы &/ >+ Как правило, использование расщепления при решении многомерных задач (§ 8.2) приводит к трехдиа-гональным системам уравнений, которые могут эффективно решаться при помощи алгоритма Томаса.

6.2,3, BANFAC/BANSOL: исключение по Гауссу для узколенточных матриц

В тех случаях, когда матрица А имеет узколенточную структуру, для реализации исключения по Гауссу подходят подпрограммы BANFAC и BANSOL. Подпрограмма BANFAC (рис. 6.18) выполняет прямую прогонку или несколько повторяющихся прямых прогонок, чтобы придать матрице А верхнетреугольную форму. Подпрограмма BANSOL (рис. 6.19) преобразует правую часть (6.23), а именно вектор В, и полученная система уравнения решается путем обратной подстановки.

Для случая INT = 2 подпрограммы BANFAC и BANSOL составлены так, чтобы воспользоваться особым видом структуры матрицы А, полученным при использовании одномерных квадратичных элементов. А именно матрица А оказывается попеременно то трехдиагональной (три соседних ненулевых элемента в одной строке), то пятидиагональной (пять соседних ненулевых элементов в одной строке), что соответствует уравнениям, связанным с узлами в серединах или в углах элементов соответственно (см. рис. 5.10).

В этом случае процесс факторизации (BANFACJ состоит из двух прогонок. Первая прогонка воздействует только на пяти-диагональные строки с целью исключения лишнего коэффициента ak, Вторая прогонка приводит матрицу к верхнетреугольной форме. Конкретный вид реализации зависит от того, связаны ли первое и последнее уравнения из системы (6.23) с узлами в серединах или в углах элементов. Граничные условия Дирихле соответствуют срединным узлам; граничные же условия Неймана соответствуют угловым узлам. Читателю может оказаться полезным изучить часть программы, относящуюся

2 SUBROUTINE BANFAC(В,N,INT)

4 С FACTORISES BAND MATRIX ARISING FROM LINEAR OR QUADRATIC ELEMENTS

5 С INTO L.U

7 DIMENSION В(5,65)

8 IFdNT .EQ. 2)G0T0 2

10 С INT - 1, LINEAR ELEMENTS TRIDIAGONAL SYSTEM

11 С

12 NP N - 1

13 DO 1 J 1,NP

14 JP J + 1

15 B(2,JP) = B(2,JP)/B(3,J)

16 B(3,JP) B(3,JP) - B(2,JP)*B{4,J)

17 1 CONTINUE

18 RETURN

19 С

20 С INT > 2, QUADRATIC ELEMENTS PENTADIAGONAL SYSTEM

21 С ASSUMES FIRST EQUATION FORMED AT MIDSIDE NODE

22 С

23 2 NH N/2

24 DO 5 I 1,2

25 JS 3 - I

26 DO 4 J JS,NH

27 JA 2*(J-l)

28 IFd .EQ. 2)G0T0 3

29 С

30 С I 1, FIRST PASS, REDUCE TO TRIDIAGONAL

31 С

32 JB = JA + 2

33 B(1,JB) Bd,JB)/B(2,JB-l)

34 B(2,JB) B(2,JB) - B(1,JB)*B(3,JB-1)

35 B(3,JB) B(3,JB) - B(1,JB)*B(4,JB-1)

36 GOTO 4

37 С

38 С 1=2, SECOND PASS, REDUCE TO UPPER TRIANGULAR

39 С

40 3 JB = JA + 3

41 B(2,JB-1) B(2,JB-l)/B(3,JB-2)

42 B(3,JB-1) B(3,JB-1) - B(2,JB-l)*B(4,JB-2)

43 IF(JB .GT. N)GOTO 4

44 B(2,JB) B(2,JB)/B(3,JB-1)

45 B(3,JB) = B(3,JB) - B(2,JB)*B(4,JB-1)

46 B(4,JB) B(4,JB) - B(2,JB)*B(5,JB-1)

47 4 CONTINUE

48 5 CONTINUE

49 RETURN

50 END

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

К высказанным положениям и начинающуюся с INT = 2, после прочтения п. 6.2.4.

Примером использования подпрограмм BANFAC/BANSOL

как в случае трехдиагональных, так и попеременно трехдиаго-

16 к. Флетчер, т. 1