CC BY
55
Вычислительные технологии
Том 4, № 1, 1999
ОБ ОДНОМ ИТЕРАЦИОННОМ АЛГОРИТМЕ РЕШЕНИЯ РАЗНОСТНЫХ ЭЛЛИПТИЧЕСКИХ
УРАВНЕНИЙ*
В. Г. Зверев Томский государственный университет, Россия e-mail: vgz@niipmm.tsu.tomsk.su
Within the framework of a scalar technology of the alternating direction method an iterative algorithm for the solution of difference elliptical equations is offered. It is based on the using of two-point recurrence relation coefficients to calculating the influence of the adjacent nodes and therefore the boundary condition along the direction transverse to the grid line. The test elliptic problems show the high effectiveness of the algorithm, comparison with other methods is carried out.
Введение
Описание большинства интересных с практической точки зрения задач в области механики вязкой жидкости и тепломассообмена основывается на многомерных дифференциальных эллиптических уравнениях. Их разностная аппроксимация приводит к системам линейных алгебраических уравнений (СЛАУ), матрицы которых имеют высокую размерность 104 в двумерном случае), упорядоченную и сильно разреженную структуру. Решение СЛАУ представляет собой одну из основных трудностей, связанную с большими затратами времени ЭВМ [1]. Поэтому актуальным является построение экономичных и эффективных алгоритмов, предназначенных для численного анализа таких систем.
Итерационные методы [1-4] имеют значительные преимущества перед прямыми подходами и являются основным инструментом при численном решении СЛАУ. В практических расчетах задач гидродинамики и теплообмена широкое распространение получил метод переменных направлений и его многочисленные модификации [5-7]. Его основу составляет расщепление сложного разностного оператора на простые одномерные [8, 9] и использование скалярных прогонок для решения трехточечных уравнений вдоль координатных направлений.
В данной работе рассматривается один из вариантов такой технологии — полиней-ный метод [5], в котором используется подход Зейделя в поперечном к сеточной линии направлении. Усложнение задач математического моделирования зачастую приводит к тому, что скорость сходимости данного метода оказывается недостаточной. В наибольшей степени это проявляется при увеличении числа узлов расчетной сетки, что влечет за собой неоправданные затраты расчетного времени.
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-00971.
© В. Г. Зверев, 1999.
Известно, что улучшение сходимости итерационных методов может быть достигнуто при помощи разных приемов [1, 10]. Сильное отличие исходного и ускоряющего алгоритмов технологически является не всегда удобным. Желательно для этой цели использовать информацию, получаемую в ходе итерационного процесса. С этих позиций, опираясь на результаты исследований [11, 12], предлагается достаточно простой способ повышения эффективности полинейного метода [5], основанный на усилении эллиптичности алгоритма при реализации всех его этапов. При этом не меняется технология исходного метода и не требуется информация о свойствах и границах спектра разностного оператора.
1. Разностный аналог и процедура расчета полинейного метода
Обобщенное дифференциальное эллиптическое уравнение второго порядка, выражающее закон сохранения переменной Ф в стационарном случае в декартовой системе координат (х, у) может быть записано в виде (обозначения [5])
дриФ друФ д / дФ \ д / дФ \ п ^ п . . ^
+ -4— = тт- Г— + — Г— + БрФ + Ба, (х, у) е С (1)
дх ду дх \ дх) ду \ ду) с краевыми условиями
дФ
01 дП = У2Ф + 0-3, (х,у) е дС.
Здесь (и, у) — компоненты поля скорости, Г — коэффициент диффузии, р — плотность среды, Бр(< 0), Ба — источники порождения Ф, п — нормаль к поверхности дС, 91, 02, 93 — кусочно-непрерывные функции х, у.
Разностная форма эллиптического уравнения (1) может быть получена различными способами [1-4]. Проблемы построения разностных схем для задач конвективного тепло-массопереноса рассматривались, например, в [13]. Здесь, с целью анализа итерационного алгоритма, кратко укажем схему метода контрольного объема [5], применение которого обеспечивает выполнение интегрального закона сохранения. Интегрируя (1) по окружающему узел Р (рис. 1) контрольному объему А У, получим балансовое соотношение, которое так или иначе лежит в основе многочисленных аппроксимаций уравнения (1):
Л - Л + Ъ - Ъ = (БрФ + Ба)АУ, Зе = (риФ - ГдФ/дх)еЛе, (2)
где Л — полный поток (конвективный и диффузионный) переменной Ф, Л — площадь грани, индексы п, в, е, указывают на ее положение. Определим Л в г- и ]-направлениях как разность между центральной (г,]) и соседними по шаблону точками [5]:
Ле = ав(Фу - Фг+1,у), Лт = аш(Фг-1,у - Фу),
Зп = аИ(Фу - Фл), Л = ав(Ф^-1 - Фу).
Подставляя данные выражения в (2), получим пятиточечную аппроксимацию интегрального баланса
арФу = авФг+1,у + ашФг-1,у + аИФг,у+1 + авФг,— + Ь, 1 < г < т, 1 < ] < п, (3)
Рис. 1. Схема разностного шаблона [5]. а У — контрольный объем, J — полный поток переменной Ф, А — площадь грани.
ap — ад + aw + + as — Sp AV, b — Sc AV, aEmj — aw 1 j — 0, aNin — asu — 0, ар > ад + aw + aN + as. (4)
где индексы E, W, N, S относятся к узлам шаблона (см. рис. 1), aE, aw и aN, as — неотрицательные коэффициенты, которые определяют конкретный вид разностной схемы [5, с. 83], b — источниковый член, нижние индексы i, j у этих величин для краткости опущены. Будем считать, что область G является прямоугольной и система уравнений (3) записана с учетом (4) разностного аналога граничных условий (1). Таким образом, матрица СЛАУ (3) имеет диагональное преобладание, является положительно определенной и в общем случае несимметричной, с максимальным числом ненулевых элементов в строке равным пяти.
Следуя [5-7], рассмотрим процедуру расчета полинейного метода [5], которая состоит в последовательном решении одномерных задач. На первом этапе неявно рассматривается трехдиагональная система уравнений, например для горизонтальных линий (j — const, 1 < i < m), перебором которых от j — 1 до j — n покрывают всю расчетную область
—aw Фк+ij2 + ap фк+1/2 — aE Ф^2 — as Ф jf2 + aN Ф^.+ + b, i — 1~m, j — (5)
Затем аналогичным образом решается другая система для вертикальных линий (i — const, 1 < j < n)
—as Ф^+Л + ap Ф^+1 — aN — aw ФJ+jj + aE ФЦ1,/2 + b, j — 17n, i — 17m. (6)
Верхние индексы k, k +1/2, k +1 относятся соответственно к предыдущей итерации, первому и второму этапам текущей итерации. В итоге расчетная область "линия за линией" просчитывается вдоль каждого координатного направления. Для функции Ф с поперечного направления используются последние в итерационном смысле значения (подход Зейде-ля), так как Ф^^ в (5) и Ф^^ в (6) уже известны из решения для ниже расположенной j — 1-й строки и i — 1-го столбца соответственно.
Преимуществом данного алгоритма является применение прямого метода вдоль сеточной линии, который реализуется экономичной скалярной прогонкой [1]. В этом случае
влияние граничных условий распространяется на все внутренние узлы линии независимо от их количества [5].
В поперечном к сеточной линии направлении имеет место иная картина. Явная форма учета Фг,у+1 в (5) и Фг+1у в (6) с предыдущего этапа итерации приводит к тому, что здесь скорость передачи информации остается такой же, как в явном поточечном методе Зей-деля. Принципиальным является тот факт, что при решении уравнения (5) отсутствует влияние краевого условия при ] = п, а для уравнения (6) — с правой стороны при г = т, которое детерминирует итерационное поведение сеточной функции.
Известно [12], что при вычислении решения эллиптического уравнения (1) используется информация о всех внутренних источниках в области и граничных условиях. Поэтому способ реализации алгоритма (5), (6) [5] в недостаточной степени отражает фундаментальную природу данного уравнения, что приводит к снижению скорости сходимости. Одним из путей повышения эффективности алгоритма является усиление неявной формы учета влияния соседних узлов и граничного условия с поперечного направления.
2. Описание итерационного алгоритма
Центральным моментом, лимитирующим скорость сходимости, как следует из изложенного выше, является выбор значений Фг,у+1 и Фг+1у с той части расчетной области, в которой при выполнении этапа итерации еще не получено решение. Применение формул экстраполяции для улучшения, например, Ф гу+1 по узлам Р (г,]), Б (г,] - 1) (см. рис. 1), может усилить сходимость при монотонном изменении функции, однако не снимает вопроса о влиянии граничного условия при ] = п на решение уравнения (5) для текущей ]-линии.
Рассмотрим интегральный баланс (2) для узла Р(г,]). Из него следует, что искомое влияние проявляется через разностный поток Лп на верхней грани контрольного объема из области с линиями (] +1, ... , п). Таким образом, уравнение (5) может рассматриваться независимо от указанной области при условии определения потока. Так как Лп представляет собой двухточечную связь, то она от границы ] = п может быть передана на ]-линию. Следуя [11], определим ее в ]- и г-направлениях в виде
Фг,у+1 = & Фг3 + Пг3 , (7)
Фг+1у = 1г3 Фгу + ^ , (8)
где 7, d и п — искомые коэффициенты. Подставим (7), (8) в (5), (6) и учтем неявно Фу. В результате получим следующий вид разностных уравнений двухэтапного алгоритма:
-ашФ^2 + (ар - ам%) Ф^2 - авФ^Ц? = а*Ф^2 + (ампк3 + Ь) , (9)
-авФк+-1 + (ар - ав%+1/2) Ф^1 - аиФ^ = ашФ^- + (ав1/2 + ь) . (10)
Реккурентные коэффициенты могут быть найдены различными способами: например, на основе соотношений, в которых выделены потоки в ]- и г-направлениях соответственно
Лп - Л - БрФАУ = -(Ле - Зт) + БсАУ,
Л - Лт - Бр ФАУ = -(.1п - Л) + Бс ФАУ.
(11)
Разностный вид (11) формально следует из уравнения (3), если к его левой и правой частям, аналогично методу неполной факторизации [12], добавить соответствующие выражения:
Уц (Ф) = -0(ая + ащ )Фу, ¥%3 (Ф) = -в(аж + а5 )Фу, где в — коэффициент нижней релаксации, 0 < в < 1. В результате получим:
-аяФ^-1 + (ар - в(аЕ + ащ))Ф^ - а^Фг,7+1 = Ь + а^Фг+1,.7 + ащФг-1,7 - в(аЕ + ащ)Ф^, (12)
-ащФг—1,7 + (ар - в(аж + ая))Фу - а^Фг+1,7 = Ь + а^Фм+1 + аяФм—1 - в(а^ + ая)Фу. (13)
Таким образом, матрица системы (3) заменяется на приближенную, трехдиагональную по каждому направлению. Это позволяет по формулам прогонки (левой) найти коэффициенты искомой связи (7), (8):
^¿,7—1 = ая/ (ар - аж^), Пм—1 = (Ь' + а^)/ (аР - а^^), (14)
Ь' = Ь + аЕФг+1,7 + ащФг—- в(аЕ + ащ)Фу, аР = аР - в(аЕ + ащ), ^ = п, 2, г = 1, т.
7г—1,7 = ащ/(аР - аЕТу), ^г—1,7 = (Ь'' + а^^)/(аР - а^Ту), (15)
Ь'' = Ь + ажФг,7+1 + аяФг,7—1 - в(аж + ая)Фг7-, аР = аР - в(аж + а5), г = т, 2, ^ = 1, п.
Отметим, что Ь' и Ь'' в (14), (15) включают в себя диагональную компенсацию значений функции с соседних сеточных линий. Согласно [12], это уменьшает норму итерируемого выражения и при линейном изменении функции в расчетной области позволяет получить точные значения коэффициентов связи (£,п), (т,^).
Рис. 2. Схема итерации модифицированного метода: а — 1 этап, б — 2 этап.
Схема итерационного метода представлена на рис. 2. Приведем алгоритм его реализации при известном приближении Фк :
а) левой прогонкой по формулам (14) находятся , с использованием которых
аналогичным образом решаются уравнения (9) для всех линий ] = 1,п, определяется поле
Ф
к+1/2.
ния
б) по Фк+1/2 с помощью формул (15) вычисляются (7,^)к+1/2, затем решаются уравне-(10) для линий г = 1,т, находятся Фк7+1.
Эти два этапа составляют одну итерацию и обеспечивают чередование направлений. В случае начала прохода области с верхней горизонтальной (j = n) и с правой вертикальной линий (i = m) формулы (7), (8) и (14), (15) следует заменить на формулы правой прогонки [1]. Нетрудно видеть, что алгоритм полностью сохраняет экономичную технологию метода переменных направлений [5]. Дополнительным является определение реккурентных связей с поперечного направления, которые возмущают коэффициент ар при центральном узле P(i,j) и источниковый член Ь. Реализация такого подхода требует хранения двух массивов.
Параметр нижней релаксации 0 < 9 < 1 в уравнениях (12), (13) обеспечивает сохранение запаса устойчивости, так как уменьшение нормы выражения b' и Ь" в то же время приводит к ослаблению диагональных коэффициентов а'р и ар соответственно [12]. Предельный случай 9 = 0, имеющий абсолютную устойчивость, рассматривался в [14].
3. Тестовые задачи
Оценка практической скорости сходимости и эффективности предложенного алгоритма для решения СЛАУ (3) проводилась на основе решения тестовых эллиптических задач.
Пример 1. Рассмотрим разностную систему (3) при ар = 4, аn = as = аЕ = aw = 1, которая является аппроксимацией на равномерной сетке в единичном квадрате задачи Дирихле для уравнения Лапласа
ДФ = 0, 0 <x,y< 1, Ф|г = 1. (16)
Точное решение (16) есть Ф*, = 1, начальное приближение полагалось равным Ф0, = 0, а итерации прекращались при выполнении условия [4]
Ф - ФкГЧ/ф-1< е = 2-16, yi,j.
Для данной задачи известны теоретические оценки скорости сходимости многих итерационных методов [1, 4].
Пример 2. Рассмотрим задачу Дирихле для уравнения Пуассона с переменными коэффициентами [1]
д ( дФ\ д ( дФ\ .чл,,
аХ ("1дх) + ду(а2ду) = ^ х-у е (0Л)- Ф|г = »• (17)
а1(х, у) = 1 + C[(x - 0.5)2 + (y - 0.5)2], а2(х,у) = 1 + C[0.5 - (х - 0.5)2 - (у - 0.5)2], 1 = 1 < аа(х,у) < 2 = 1 + 0.5C, C = const, а = 1, 2.
Точным решением (17) является пробная функция Ф(х,у) = х(1 - х)у(1 - у), по ней подбирается правая часть (р(х,у). Для аппроксимации (17) использовались равномерная сетка и = {х, = (ih,jh), 0 < i,j < N, h = 1/N} и разностная схема (3), которая приведена в [1]. Константой C варьировался диапазон изменения коэффициентов а1, а2. Начальное приближение полагалось равным Ф0, = 1.
Сходимость итерационного процесса отслеживалась по изменению нормы векторов невязки \ \ Як||/1|Я01| и погрешности \ \ гк\ \ / \ \ £011, вычислялась величина qk средней скорости сходимости [4, 15]:
\\Rk
1/2
к\2
Y^ (aEФк+и + aw$k-i;i + амФ^+1 + asФ^- + b - apфк.)
г,]
1/2
- ф*.
г] г]
i,j = 1,N - 1, qk = - ln(\\zk\\/\\z0\\)/k.
г,3
(18)
где верхний индекс k относится к номеру итерации, а * — к точному численному решению, которое было заранее найдено с высокой точностью (при \\Rk\\/\\R0\\ = 10-12).
Реализация алгоритма состояла в решении уравнений (9) вдоль координаты x, а затем (10) — вдоль у. Для сравнения использовался метод SIP [16], который совпадает с явной схемой неполной факторизации Н. И. Булеева [12, с. 163] с параметром 9 = 0. В этом случае итерация состояла в проходе области сначала по диагонали "SW - NE", затем "SE -NW" (см. рис. 1). В рассматриваемых алгоритмах не ставилась цель применять оптимальные итерационные параметры, поэтому коэффициент 9, если не оговорено специальным образом, в (14), (15) полагался равным предельному значению 1.
2
4. Результаты численных расчетов и их анализ
В таблице 1 приведены данные по числу итераций, необходимых для решения первого примера, в скобках указана величина скорости сходимости qk (18). Видно, что предложенный алгоритм позволяет на любой сетке получить решение Ф = const за одну итерацию. Значительное сокращение числа итераций (см. также [4] и тестированные там методы) обусловлено тем, что трехточечные уравнения (12), (13) удовлетворяют точному решению в отличие от исходного полинейного метода [5], где имеет место прямая зависимость правой части от значений функции с соседних сеточных линий. Поэтому точными получаются рекуррентные коэффициенты и, следовательно, разностный поток с поперечного сеточного направления.
Таблица 1
Расчетная Метод SIP [16] Исходный Данная
сетка N х N Зейделя [4] метод [5] работа
20 х 20 319 (0.025) 60 (0.17) 96 (0.10) 1
50 х 50 1531 (0.004) 297 (0.028) 473 (0.016) 1
Изменение векторов невязки \\Rk \\ и погрешности \\zk\\, а также скорости сходимости qk в зависимости от числа k итераций при решении уравнения Пуассона для более сложной зависимости Ф = (,) представлено на рис. 3. Параметры задачи полагались равными c2/c1 = 32, N = 32. Принята следующая нумерация кривых: 1 — исходный полинейный метод (ПМ) [5], 2, 3 — алгоритм данной работы при 9 = 0 и 9 = 1 (пунктирная кривая — для случая 9 = 9opt), 4 — метод SIP [16], 5 — поточечный метод Зейделя [15]. Для кривых 3 в отличие от других траекторий имеет место резкое падение невязки (рис. 3, а) и погрешности (рис. 3, б) после одной итерации более чем на два порядка. Сравнение с
б
lg|kfc 11/11 К
О 40 80 120 160 т 200
к
Рис. 3. Динамика сходимости итерационного процесса в зависимости от числа к итераций. \\Rk|| — невязка (а), \\zk\\ — погрешность, qk — скорость сходимости (б). 1 — ПМ [5], 2, 3 — данная работа при в = 0 и в = 1, (пунктир — в = #0pt), 4 — SIP [16], 5 — метод Зейделя [15];
\\R0\\ = 197, \\z0\\ = 30.1, c2/ci = 32 (сетка 32 х 32).
кривой 2 (в = 0) показывает, что полученный результат целиком обусловлен уменьшением нормы итерируемого выражения при вычислении источникового члена Ь' в (14) и Ь'' в (15). Скорость сходимости (кривая 3' рис. 3, б) данного алгоритма выше 0.7, что на порядок больше значения = 0.073 исходного метода [5] (кривая 1'). Приведем асимптотическое значение для других кривых 2', 4', 5': 0.13, 0.11 и 0.0093 соответственно. Последнее значение (метод Зейделя) совпадает с теоретической оценкой ~ п2Л,2 [1, 4]. Видно, что применение этого поточечного метода крайне неэффективно для решения эллиптических уравнений.
На последовательности сеток было исследовано влияние количества узлов расчетной сетки и анизотропии коэффициентов ах, а2 уравнения (17) на сходимость итерационного процесса. В таблице 2 приведены сравнительные данные по числу итераций, необходимых для уменьшения невязки в 104 раз, в скобках указано значение . Из таблицы следует, что предложенный алгоритм, в отличие от полинейного метода [6], практически не зависит от отношения с2/сх. По числу итераций его эффективность выше в 10-30 раз в зависимости от шага расчетной сетки. В наибольшей степени это проявляется при равных коэффициентах ах и а2, что соответствует сильной эллиптичности задачи. Важным является то обстоятельство, что скорость его сходимости зависит от числа узлов сетки как 0(1/Ж).
Введенный в уравнения (12), (13) параметр нижней релаксации 0 < в < 1 позволяет получить весь спектр сходимости в области между кривыми 2 (в = 0) и 3 (в = 1) (см. рис. 3). Однако его влияние является немонотонным, особенно при 0.9 < в < 1. Это показывают кривые 1, 2 на рис. 4, построенные по данным табл. 2. Оптимальное значение в для данной задачи на сетке 32 х 32 составляет = 0.992. Применение позволяет еще более усилить сходимость, это хорошо показывают пунктирные кривые 3 рис. 3, а также данные последней колонки табл. 2 для других расчетных сеток. Определение в^ является самостоятельным вопросом, в практических расчетах может использоваться широкий диапазон 0.5 < в < 1. Следует заметить, что экспериментальный подбор в не представляет особых трудностей.
В таблице 3 представлено сравнение эффективности рассмотренных алгоритмов. Число итераций (по данным табл.2), процессорное время выполнения одной итерации и задачи в целом отнесены к соответствующим параметрам исходного полинейного метода [5]
Таблица 2
Количество итераций, необходимых для выполнения условия \\Кк||/||Я°|| = 10 4,
при решении уравнения Пуассона [1]
Расчетная С2/С1 Я1Р [16] Исходный Данная работа
сетка N х N метод [5] в = 1 в - ворЬ
2 88 149 (0.042) 12 (0.66) 4 (2.44) в = 0.992
32 х 32 32 55 (0.114) 87 (0.073) 9 (0.91) 5 (1.96)
512 51 80 (0.081) 9 (0.95) 5 (1.96)
2 290 490 (0.011) 20 (0.35) 6 (1.60) в = 0.998
64 х 64 32 181 (0.029) 285 (0.019) 16 (0.47) 6 (1.55)
512 166 261 (0.021) 16 (0.49) 6 (1.53)
2 916 1544 (0.003) 32 (0.19) 8 (1.15) в = 0.9994
128 х 128 32 576 (0.007) 907 (0.005) 27 (0.25) 8 (1.08)
512 529 833 (0.005) 27 (0.26) 8 (1.06)
2
О 0.2 0.4 0.6 0.8 0.9 1.0
Рис. 4. Влияние параметра в на сходимость итерационного процесса при решении уравнения
Пуассона [1].
к - число итераций (1), - скорость сходимости (2). \\Ек ||/||Я°|| = 104, ||Я°|| = 197, с2/с1 = 32, (32 х 32).
Таблица 3
Коэффициент Исходный SIP [16] Данная
метод [5] работа
Kt 1 1.10 1.15
Kiter 1 0.63 0.10/0.025
Kef = Kt ■ Kiter 1 0.70 0.12/0.03
(коэффициенты Kiter, Kt, Kef). Нетрудно видеть, что трудоемкость одной итерации предложенного алгоритма увеличивается примерно на 15% (Kt). Высокие значения коэффициента Kiter (до 0.025 на сетках 128 х 128) обеспечивают более чем 10-кратное (1/Kef) уменьшение времени расчета. Причем наибольший эффект достигается на подробных сетках, что очень важно при практическом моделировании. Как видно из табл. 2 и 3, значительное преимущество по числу итераций имеет место и по отношению к методу SIP [16].
5. Заключение
Таким образом, проведенные расчеты показывают высокую эффективность предложенного алгоритма. Он полностью сохраняет скалярную технологию метода переменных направлений и область его применимости, не требует дополнительной информации о спектре разностного оператора. Значительное ускорение сходимости достигнуто за счет усиления "эллиптичности" алгоритма посредством выделения разностного потока и учета влияния краевого условия с поперечного к сеточной линии направления, а также уменьшения нормы итерируемого выражения при расчете рекуррентных коэффициентов. Его применение позволяет на порядок сократить число необходимых итераций по сравнению с полинейным методом [6], при этом на выполнение одной итерации требуется незначительное увеличение затрат времени и памяти ЭВМ. Полученные результаты расширяют возможности метода переменных направлений при решении разностных эллиптических уравнений.
Автор признателен проф. А. М. Гришину за внимание к работе, а также благодарит
рецензента за сделанные замечания и рекомендации, которые привели к улучшению содержания публикуемой статьи.
Список литературы
[1] Самарский А. А., Николаев е. С. Методы решения сеточных уравнений. Наука, М., 1978.
[2] Марчук Г. И. Методы вычислительной математики. Наука, М., 1989.
[3] Хейгеман Л., Янг Д. Прикладные итерационные методы. Мир, М., 1986.
[4] Ильин В. П. Разностные методы решения эллиптических уравнений. Изд-во НГУ, Новосибирск, 1970.
[5] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. Энергоатомиздат, М., 1984.
[6] Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Т. 1, Мир, М., 1990.
[7] Флетчер К. Вычислительные методы в динамике жидкостей. Т. 1, Мир, М., 1991.
[8] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Наука, Новосибирск, 1967.
[9] Марчук Г. И. Методы расщепления. Наука, М., 1988.
[10] СлЕпцов А. Г. Об ускорении сходимости линейных итераций. Моделирование в механике, 3(20), №5, 1989, 118-125.
[11] Четверушкин Б. Н. Об одном итерационном алгоритме решения разностных уравнений. ЖВММФ, 16, №2, 1976, 519-524.
[12] Булеев Н. И. Пространственная модель турбулентного обмена. Наука, М., 1989.
[13] Вабищевич П. Н., Самарский А. А. Разностные схемы для нестационарных задач конвекции-диффузии. ЖВММФ, 38, №2, 1998, 207-219.
[14] Зверев В. Г., Зверева Л. В. Модифицированный полинейный метод решения двумерных сеточных уравнений. Деп. ВИНИТИ, №76-В92 от 08.01.92.
[15] Самарский А. А., Гулин А. В. Численные методы. Наука, М., 1989.
[16] Stone H. L. Iterative solution of implicit approximation of multidimensional partial differential equations. SIAM J. Numer. Anal., 5, 1968, 530-558.
Поступила в редакцию 23 сентября 1997 г., в переработанном виде 15 октября 1998 г.
