после чего с помощью обратных подстановок можно рассчитать V. Ясно, что чем более широкую ленточную структуру будет иметь матрица А, тем менее экономичным будет вышеописанный алгоритм.

Если внутри самой ленты вкраплено достаточно большое число нулей, то, как правило, более экономичными оказываются некие альтернативные приемы обработки разреженных матриц [Jennings, 1977а], хотя они вызывают больше трудностей при программировании.

6.2.5. Блочные трехдиагональные системы

Описание алгоритма Томаса для решения трехдиагональных систем дается в п. 6.2.2. Как правило, этот алгоритм следует применять в тех случаях, когда единственное исходное уравнение подвергается дискретизации с помощью неявного алгоритма (§ 7.2). Однако многие задачи гидроаэродинамики описываются системами уравнений (§ 10.2 и гл. И). Попытки реализовать неявные алгоритмы приводят тогда, как правило, не к скалярным трехдиагональным системам уравнений, а к блочным трехдиагональным структурам. Однако алгоритм Томаса без: труда распространяется и на действия в таких ситуациях.

Блочно-трехдиагональная система уравнений, эквивалентная системе (6.28), может быть представлена в виде

bi ci 2 Ьг 02

а/ bi Ci

ал 1 Ьл1 C/v i

.d;v

(6.36)

где a,-, b;, Ci - субматрицы размера М X Л1, a V/ и d, суть М-компонентные субвекторы. Число М соответсвует числу уравнений, записываемых в каждой точке сетки; например, для трехмерного течения вязкой сжимаемой жидкости M = 5 (п. 11.6.3). Отсюда следует, что V/ -субвектор решения, связанный с конкретной точкой сетки. Уравнение (6.36)-система состоящая из N блоков уравнений, причем каждый блок связанный с конкретной точкой сетки, включает в себя М уравнений.

Решение системы (6.36) идет по методике, очень близкой методике решения (6.28). Сначала трехдиагональная матрица

блоков из (6.36) преобразуется к верхнетреугольной форме за счет исключения субматриц а/. По аналогии с (6.29) первый блок уравнений дает

c; = (bi)-c d[ = (bi) (6.37)

тогда как для блока общего вида имеем

ь; = ь,-а,с; с; = (ь;)-с d; = (b;)-{d,-a,d; ,}. (6.38)

В выражениях (6.37) и (6.38) фигурируют явные выражения обратных матриц. На практике более экономичный путь связан с решением М-компонентной субсистемы; например, уравнение Ьс. = с. решается относительно с. После выполнения операций (6.37) и (6.38) уравнение (6.36 ) приобретает верхнетреугольную форму при замене с, и d/ нас и и при замене Ь-на единичную матрицу I.

Второй этап, эквивалентный применению формул (6.31), требует обратных подстановок по формулам

V;. = d;, v, = d;-c;v,,. (б.зэ)

Как показывает подсчет, для реализации блочного алгоритма Томаса требуется 0{5NM/3) операций, что явно предпочтительно числу 0((ЛЛ1)/3), требуемому для реализации полного исключения по Гауссу. Однако если бы оказалось возможным расщепить блочно-трехдиагональную систему (6.36) на М скалярных трехдиагональных систем, то полный подсчет показал бы 0(5ЛЛ1) операций, т. е. грубо говоря, уменьшение времени счета в раз при М. В результате возникает

стимул к тому, чтобы построить неявные алгоритмы, которые позволили бы части уравнений отделиться от системы в целом. Вариант блочного алгоритма Томаса, описанный выше, излагается по книге [Isaacson, Keller, 1966].

6.2.6. Прямые алгоритмы решения уравнения Пуассона

Уравнение Пуассона, так же как и связанное с ним уравнение Лапласа, достаточно часто встречается в гидроаэродинамике, чтобы имело смысл исследовать возможности разработки специальных процедур для решения дискретизированного уравнения Пуассона. Трехточечная конечно-разностная дискретизация двумерного уравнения Пуассона приводит к системе уравнений типа (6.23), для которой матрица А имеет форму, показанную на рис. 6.1.

Для случая однородной сетки (Ajc = A(/), покрывающей прямоугольную область Ojcl, Oyl, один из блоков

системы уравнений может быть записан в виде

IV, , + gV, + IV , = h /=1, (6.40>

где М - число активных точек в направлении оси у. Каждый вектор \k содержит N неизвестных значений У/, и, связанных с одной из линий сетки {у = Ук)у направленной вдоль оси х. Матрица g является трехдиагональной, имеет порядок и отличные от нуля элементы 1, -4, 1, относящиеся соответственно к точкам (/- 1, fe), (/, k) и (/-+- I, k). Вектор hk содержит все те вклады в В из уравнения (6.23), которые вносит в правук> часть единственная линия сетки, идущая параллельно оси х, т. е. у = ук.

При решении уравнения (6.40) особенно эффективными оказываются два метода - метод циклической редукции и метод разложения в ряд Фурье. Предпочтительная стратегия [Swartztrauber, 1977] сводится к тому, чтобы вначале с помощью циклической редукции уменьшить размер М в уравнениях (6.40), после чего строить решение редуцированной системы с помощью разложения в ряд Фурье.

В процессе циклической редукции блок, расположенный в k-Pi (четной) строке уравнения (6.40), умножается на -g и складывается с блоками в {k-1)-й и {к+1)-й строках с целью исключения V-i и V+i. В результате получим

IV, 2 + gv, -f IV,2 = (6.41>

g() = 21 - gg = (V2 I -f g) (V21 ~ g), 4> = h, ,-gh,-fh,,.

Этот процесс можно повторить требуемое число раз. После осуществления / редукций уравнение (6.40) приобретает форму

IV, 2/ + gV + W,,i = h (6.42>

= - Е is - т = 2 cos l{2k - 1) я/2+1].

В принципе редукция может быть определена до тех пор, пока граничные условия не позволят определить Wk-2i и [Dorr, 1970]. Затем с помощью обратных подстановок будут определены все промежуточные векторы V. На практике, однако, оказывается более эффективным завершить данный процесс после / редукций, переключаясь на применение рядов Фурье