CC BY
167
Вычислительные технологии
Том 10, № 5, 2005
МНОГОСЕТОЧНЫЙ МЕТОД РЕШЕНИЯ СИЛЬНО НЕСИММЕТРИЧНЫХ СИСТЕМ*
Е. М. Андреева, Г. В. Муратова ЮГИНФО РГУ, Ростов-на-Дону, Россия e-mail: andreeva@rsu.ru, muratova@rsu.ru
A multigrid method for solution of the systems of linear algebraic equation with strongly nonsymmetric matrix, which obtained after difference approximation of the convection-diffusion equation with dominant convection, is proposed. Modifications of the multigrid method, which use iterative methods from the class of triangular skew-symmetric methods as the smoothers are provided.
Введение
Для решения задач линейной алгебры существует множество различных численных методов, которые непрерывно усовершенствуются и модифицируются. Активно разрабатываются новые методы. В последние годы одним из эффективных и довольно универсальных итерационных методов решения систем линейных алгебраических уравнений стал многосеточный метод (MGM). Он не является строго фиксированным методом. Это скорее некий шаблон, сборная конструкция метода, эффективность которого зависит от адаптации его компонентов к решаемой задаче [1].
Одим из компонентов многосеточного метода — сглаживатель или базовый итерационный метод. В данной статье предложен многосеточный метод со сглаживателями специального вида, которые эффективно решают сильно несимметричные системы линейных алгебраических уравнений, появляющиеся после центрально-разностной аппроксимации уравнения конвекции-диффузии с преобладающей конвекцией. Данное уравнение является математической моделью многих природных процессов и явлений, в основе которых лежит конвективно-диффузионный перенос в средах с преобладающей конвекцией.
Приводится исследование предложенных модификаций многосеточного метода с различными сглаживающими процедурами [2].
1. Постановка задачи
Рассматривается система линейных алгебраических уравнений
Au = f, (1)
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 03-01-00005, № 04-01-96807) и программы "Университеты России" (грант № УР.03.01.273). ( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.
в которой матрица А является разреженной, сильно несимметричной и не имеет диагонального преобладания. Матрицу А представим в виде суммы ее симметричной и косо-симметричной части:
А = Ао + Ах,
где А0 = 1 (А + А*) — симметричная часть, а Ах = ^ (А — А*) — кососимметричная часть матрицы А. А является сильно несимметричной, что означает выполнение в некоторой матричной норме следующего неравенства:
1|Ао||* << ||Ах||*.
Кроме того, матрица А является диссипативной, т. е. ее симметричная часть положительно определена:
Ао = АО > 0, Ах = —АХ. Кососимметричную часть матрицы А — матрицу Ах — представим в виде
Ах = К1о„ + Кир и Кир = — К]*оте,
где К]о„ и Кир — соответственно нижне- и верхнетреугольные части кососимметричной матрицы Ах .
Рассматриваемая система (1) решается многосеточным методом, в котором в качестве сглаживателей используются треугольные кососимметричные методы.
2. Многосеточный метод
Алгоритмическую основу многосеточного метода составляют два приема. Первый из них заключается в последовательном решении сеточных задач, аналогичных исходной задаче, но решаемых на более крупных сетках. Процесс решения начинается с самой грубой сетки, где оно может быть осуществлено довольно экономично. Затем полученное решение интерполируется на более мелкую сетку и используется в качестве начального приближения в каком-либо итерационном процессе. Такое начальное приближение уже имеет точность, близкую к удовлетворительной. В итоге в итерационном процессе требуется существенно меньше итераций для доведения точности до требуемых пределов.
Второй прием выдвинут в работе Р. П. Федоренко [3]. Он основан на свойстве некоторых итерационных методов сходиться с высокой скоростью на первых итерациях, замедляясь в дальнейшем. Это происходит за счет быстрого подавления высокочастотных компонентов Фурье начальной ошибки в разложении по базису из собственных векторов, т.е. метод быстро гасит высокочастотные гармоники ошибки, что позволяет за несколько итераций существенно подавить их вклад. Низкочастотные гармоники сходятся гораздо медленнее и поэтому будут составлять подавляющую часть ошибки.
Одна итерация многосеточного метода включает в себя четыре наиболее важных этапа: smoothing (сглаживание), restriction (ограничение, проекция), prolongation (продолжение, интерполяция) и coarse grid correction (грубосеточная коррекция). Для погашения высокочастотных компонентов погрешности, появляющихся после интерполяции на более мелкую сетку, можно сделать несколько сглаживающих итераций на этом уровне. Это так называемое post-smoothing — последующее сглаживание. Вычисления прекращаются по достижении заданной точности.
3. Сглаживающие процедуры
Важным этапом в развитии многосеточных методов является разработка сглаживающих методов. Сглаживающий метод — это центральная компонента многосеточного алгоритма, он является наиболее зависимой частью многосеточного метода от типа решаемой задачи. Роль сглаживающих методов (методов релаксации) заключается в том, что они должны не столько уменьшать полную ошибку, сколько сглаживать ее (а именно подавить высокочастотные гармоники ошибки) так, чтобы ошибка могла быть хорошо приближена на грубой сетке. Стандартными сглаживающими методами являются линейные итерационные методы, например метод Гаусса — Зейделя. Альтернативой могут служить:
— итерационный метод Ричардсона;
— метод Гаусса — Якоби;
— симметричный метод Гаусса — Зейделя;
— метод Гаусса — Зейделя с красно-черным упорядочиванием;
— четырехцветный метод Гаусса — Зейделя;
— метод попеременных направлений Гаусса — Зейделя;
— метод неполной факторизации;
— специально адаптированный SOR.
Для решения поставленной задачи в данной работе используются три метода из класса треугольных кососимметричных методов — TKM, TKM1 и TKM2. TKM — это треугольный кососимметричный метод, TKM1 и TKM2 — две его модификации. Эти методы применимы к задачам, в которых начальная матрица A не имеет диагонального преобладания. Возможность применения методов для таких задач достигается с помощью специального выбора оператора B [2, 4]. Рассмотрим более подробно структуру метода. Любой итерационный метод можно записать в канонической форме
BVn+1 - Уп + Ayn = f, n = 0,1, 2,... (2)
т
Для TKM оператор B строится следующим образом:
B = E + 2tK¡ow или B = E + 2tKUp, (3)
для TKM1
B = aE + 2K[o w или B = aE + 2KU p, (4)
для TKM2
B = aiE + 2K[o w или B = a¿E + 2KUp, (5)
где т — скалярный итерационный параметр. Напомним, что Kiow и Kup — соответственно нижне- и верхнетреугольные части кососимметричной матрицы A1. Параметры предложенных методов a, ai > 0 выбираются по формулам
a = ||M|| ,
n
ai = ^2 \mij|, i = 0,..., n, j=0
где M = {mij }П — симметричная матрица, которая строится по формуле M = A0 + Kup — K¡ow; n — размерность матрицы A. Эти методы предложены Л.А. Крукиером. Были проведены теоретические исследования и доказана сходимость этих методов [5].
Сходимость многосеточного метода с предложенными треугольными кососимметрич-ными сглаживателями можно найти в [6, 7].
(6)
(7)
4. Функции restriction и prolongation
Рассмотрим функции, описывающие процедуру перехода с более мелкой сетки на более грубую и обратно. Для этого построим две сетки — Hi и где I — уровень сетки в
многосеточном методе. Мы считаем, что П и — квадратные сетки с шагами hi и
hi_1 = 2hi. Так же можно рассмотреть случай, когда шаги сеток в по направлениям x и y различны, но при этом должно выполняться соответствие h^/h^ = hy—i/hy = 2 [1].
Функция prolongation осуществляет переход с грубой сетки Qi_i на более мелкую сетку П1. Точки (0, 0) , (0, 2hi) , (2hi, 0) , (2hi, 2hi) G H П Hi—1. Для заданного vi—1 мы должны посчитать vi = prolongation^^ ). В точках грубой сетки значения остаются неизменными:
vi (0, 0) = vi_1 (0, 0), vi (0, 2hi) = vi_1 (0, 2hi), vi (2hi, 0) = vi_1 (2hi, 0), vi (2hi, 2hi) = vi_1 (2hi, 2hi).
Точки (0, hi) , (hi, 0) , (2hi, hi) , (hi, 2hi) считаются по формулам
vi (0, hi) = 1 (vi_1 (0, 0) + vi_1 (0, 2hi)), vi (hi, 0) = 1 (vi_1 (0, 0) + vi_1 (2hi, 0)), vi (2hi, hi) = 1 (vi_1 (2hi, 0) + vi_1 (2hi, 2hi)) , vi (hi, 2hi) = 2 (vi_1 (0, 2hi) + vi_1 (2hi, 2hi)). Для определения точки vi(hi, hi) есть две возможности:
vi (hi, hi) = 1 (vi_1 (0, 0) + vi_1 (0, 2hi) + vi_1 (2hi, 0) + vi_1 (2hi, 2hi)) ; (8)
vi (hi, hi) = 1 (vi_1 (0,0) + vi_1 (2hi, 2hi)). (9)
Аналогичным образом vi определяется в других ячейках. Интерполяция с помощью формул (6)-(8) называется девятиточечным продолжением, а с помощью формул (6), (7), (9) — семиточечным продолжением.
Функция restriction выполняет обратную процедуру — переход с мелкой сетки Hi на более грубую сетку Hi_1. Для заданного di необходимо построить di_1 = restriction(di). Общая формула имеет вид
1
restriction (di)(x,y)= ^^ di (x + ahi,y + eh).
а,в=_1
Здесь используется шаблон
0"_1,1 ^0,1 ^1,1 ^_1,0 ^0,0 ^ 1,0 ^_1,_1 ^0,_1 0"1,_1
Так же различаются девяти- и семиточечное ограничения, они соответственно задаются шаблонами
if1 2 M 1 1 1
— 242 , - 121
14l 2 J , 4 110
5. Численные результаты
В качестве модельной задачи для исследования свойств МОМ с предложенными сглажи-вателями рассматривалась стационарная задача конвекции-диффузии с преобладающей конвекцией:
1' V „ (х) ди(х) + д(^»(х)и(х)м - V —(ди(х) ^ = ^(х)
2^а(х) дх + дх дх \ дх = 1 (х)
2 \ ^д^Са ^д^Са } Р^е ^са V ^Са
ч«=1 / а=1 4
и(х)|дс = g(x). (10)
Требуется найти решение задачи в области Д = [0,1] х [0,1]. Уравнение (10) содержит малый параметр при старшей производной. Переменные коэффициенты при первых производных определяют скорость движения среды, V = — вектор скорости дви-
жения среды. Мы ограничиваемся рассмотрением движения неразрывной несжимаемой
2 дь
жидкости, для которой divV =/ 7—^ = 0. Выбирались различные коэффициенты при
' ^ ппг
дса
а= 1
первых производных и разные числа Пекле от 10 до 100 000.
Для аппроксимации данного дифференциального уравнения использовался метод конечных разностей с центрально-разностной аппроксимацией первых производных на сетке 33 х 33. Рассматривалась модельная краевая задача с известным точным решением и (х,у) = 8т(пс) вт(пу) ехр(ху). Решение находилось с заданной точностью е = 10-6.
Задача (10) решалась многосеточным методом, где в качестве сглаживателей использовались ТКМ, ТКМ1 и ТКМ2. Для метода ТКМ итерационный параметр т выбирался из условия диагонального преобладания матрицы В. Количество сглаживающих итераций в МОМ = 15.
В табл. 1 представлены четыре варианта задания коэффициентов при конвективных членах, которые были использованы при проведении численных экспериментов. В табл. 2 представлены результаты сравнения предложенных модификаций многосеточного метода с треугольными кососимметричными методами в качестве сглаживателей. Проводится также сравнение методов с многосеточным методом со сглаживателем методом Зейделя. Указано количество итераций исследуемых методов, символ А означает, что на данной задаче метод не сошелся за количество итераций больше 5000, а символ N, что метод разошелся. Приведен коэффициент кососимметрии К = Рек/2, оказывающий большое влияние на скорость сходимости многосеточного метода, здесь к =1/32 — шаг по сетке.
Более полное исследование многосеточного метода можно найти в [7, 8].
Таблица1 Коэффициенты при конвективных членах
Задача Р (х,у) Я(х,у)
1 1 -1
2 1 - 2х 2у - 1
3 х + у х - у
4 8ш2пх —2пу сов 2пх
Таблица2
Количество итераций многосеточного метода с различными сглаживателями
Pe MGM MGM MGM MGM K = Peh/2
(Seidel) (TKM) (TKM 1) (TKM 2)
Задача 1
10 13 35 30 30 0.1562
102 63 7 5 5 1.5625
103 A 13 9 9 15.625
104 A 75 58 58 156.25
105 A 535 460 430 1562.5
Задача 2
10 22 72 53 50 0.1562
102 18 24 19 14 1.5625
103 A 16 12 6 15.625
104 A 59 51 32 156.25
105 A 384 522 165 1562.5
Задача 3
10 16 43 35 35 0.1562
102 23 10 7 5 1.5625
103 N 16 12 8 15.625
104 N 74 55 36 156.25
105 N 570 441 258 1562.5
Задача 4
10 17 42 32 27 0.1562
102 N 16 12 7 1.5625
103 N 30 22 10 15.625
104 N 193 159 65 156.25
105 N A A A 1562.5
Заключение
Вычислительные эксперименты показали, что предложенные модификации многосеточного метода со сглаживателями TKM, TKM1, TKM2 эффективны для решения задач конвекции-диффузии с преобладающей конвекцией. Многосеточный метод со стандартным сглаживателем (методом Зейделя) неэффективен для решения поставленной задачи. Многосеточный метод со сглаживателями TKM1 и TKM2 более эффективен для большинства задач, чем многосеточный метод со сглаживателем TKM. Наиболее эффективным для исследуемой задачи является многосеточный метод со сглаживателем TKM2.
Список литературы
[1] Hackbusch w. Multigrid Method and Application. Berlin: Springer-Verlag, 1985. 377 p.
[2] Krukier l. Special preconditions for iterative Solution of strongly nonsymmetric linear systems // Conf. PRISM'97. Nijmegen, 1997. P. 107-119.
[3] Федоренко р.п. Релаксационный метод решения разностных эллиптических уравнений // Журн. вычисл. математики и мат. физики. 1961. Т. 1, № 5. C. 922-927.
[4] Крукиер л.а. Неявные разностные схемы и итерационный метод их решения для одного класса систем квазилинейных уравнений // Изв. вузов. Математика. 1979. № 7. C. 41-52.
[5] Krukier l.a., Chikina l.g., Belokon t.v. Triangular skew-symmetric iterativ solvers for strongly nonsymmetric positive real linear system of equations // Appl. Num. Math. 2002. Vol. 41. P. 89-105.
[6] Muratova g., Krukier l. Multigrid method for the iterative solution of strongly nonselfadjoint problems with dissipative matrix // Conf. AMLI'96. Nijmegen, 1996. Vol. 2. P. 169-178.
[7] Muratova g.v., Andreeva e.m. Solution of convection-diffusion problem by multigrid method with different smoothers // Intern. Conf. Comput. Math. Novosibirsk: IMC&MG, 2002. Vol. 2. P. 649-654.
[8] Андреева е.м., Муратова г.в. Модификация многосеточного метода для решения задач конвекции-диффузии с малым параметром // Всерос. молодежная научная школа-конф. "Численные методы решения линейных и нелинейных краевых задач". Казань: ДАС, 2001. Т. 11, № 11. С. 132-141.
Поступила в редакцию 18 мая 2005 г., в переработанном виде — 30 июня 2005 г.
