Вариант многосеточного метода с полуукрупнением сетки Текст научной статьи по специальности «Математика»

УДК 519.615.5

А. А. Буздин, С. С. Дедух

ВАРИАНТ МНОГОСЕТОЧНОГО МЕТОДА С ПОЛУУКРУПНЕНИЕМ СЕТКИ

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

A variant of multigrid method for solving large systems of linear equations with block tridiagonal matrices that have higher robustness properties is presented.

In this method the construction of coarse grid correction operators is based on approximation of the Schur complement. Numerical experiments show high efficiency of presented methods.

Ключевые слова: система линейных уравнений, многосеточный метод, крупносеточная коррекция, дополнение Шура.

Keywords: system of linear equations, multigrid method, coarse grid correction, Schur complement.

В работе рассматривается один из методов решения больших разреженных систем линейных алгебраических уравнений

где А = blocktridiag(-Li,Di,-Я), Li,Di,Я єЯпхп, і = 1,...,п — положительно определенная, блочно-трехдиагональная матрица из яп'тхп'т. Такие системы обычно получаются в результате дискретизации краевых задач методом конечных разностей или методом Галеркина.

Одним из наиболее эффективных методов решения такого рода задач является многосеточный метод. Главная его отличительная черта — независимая от числа узлов сетки скорость сходимости. Это позволяет решать задачи за время О(п) [1], т. е. затрачивая на каждый узел сетки строго определенное количество действий.

В основе алгоритма лежит снижение низкочастотной части невязки благодаря решению систем на более грубых сетках. Для этого вводится последовательность сеток Хп з Хп-1 з ... з Х0, каждой из которых соответствует своя матрица Ах, I = 0,...,п.

В классическом варианте многосеточного метода для построения матриц на грубых сетках используется метод Галеркина А;-1 = гАгр [2], где г — оператор сужения; р — продолжения. При этом в двумерном случае число

Введение

Au = f,

(1)

Вестник Российского государственного университета им. И. Канта. 2009. Вып. 10. С. 74 — 81.

узлов на каждой более крупной сетке уменьшается примерно в четыре раза. В варианте многосеточного метода с полуукрупнением сетки число узлов на более грубых сетках уменьшается в два раза, и только в одном направлении (уменьшается только количество строк или столбцов сетки). Шаблон оператора продолжения для такого варианта метода выбирается,

как и в одномерном случае, в виде р =

[2]. Оператор сужения

обычно берется сопряженным оператору продолжения, а матрицы на более грубых сетках строятся методом Галеркина. Для такого метода в работе [3] получена следующая оценка скорости сходимости:

||ММСМ|| <_^, (2)

II 11а ^ +1

где V — количество сглаживающих итераций блочного метода Гаусса — Зейделя на каждом уровне.

Однако описанный в работе [3] вариант многосеточного метода применим только для сеток размерности (2п -1). В предлагаемой работе рассмотрен вариант построения матриц на грубых сетках, основанный на исключении отдельных строк или столбцов сетки и аппроксимации полученной системы уравнений. Приводимый подход позволяет применять описанные ранее алгоритмы к сеткам произвольных размерностей, а также повысить скорость сходимости метода.

Опишем способ построения матриц на грубых сетках, являющийся основой предлагаемого метода.

Рассмотрим векторы и и / в выражении (1) как блочные вектор-столбцы из Япт : и = (и1,и2,...,ип)Т, / = (Р1,Р2,...,Рп)т . Опишем процесс построения грубой сетки, не содержащей неизвестные, образующие столбец с номером к, 1 < к < п, т. е. ик. Рассмотрим уравнения с участием ик:

-^к-2ик-2 + Бк-1ик-1 - Як-1ик = Р-1,

-^к-1ик-1 + Бкик - Якик+1 = Р,

-^кик + Бк+1ик+1 - Як+1ик+2 = Рк+1.

Выразив из второго уравнения Ык и подставив в остальные, получим:

[-^к-2Ык-2 + Фк-1 - ^к-1^кг^к-1)ик-1 - ^к-^-^кЫк+1 = ^к-1 + ^к-^к^к, 1-^кЦ/^к-1Ык-1 + (^к+1 - ^к^-г^к )Ык+1 - ^к+1Ык+2 = ^к+1 + ^к^^^к.

Мы видим, что матрица системы уравнений А;-1 получается из матрицы Ах путем удаления т строк и столбцов, соответствующих Ык, и изменению четырех блоков, соответствующих Ыкк1 и Ык+1 . Запишем эти блоки в матричном виде:

Ак =

Л1-1 =

^ Бк-1 Як-1 Бк 1 ^к-1 Як-1Бк 1 Як

V-LkDk "^к-1 Бк+1 - ^к^к 1 Як )

^ Бк-1 (-1; Як)

В общем случае полученные блоки матрицы будут полностью за-

Если удастся за счет выбора параметров а1, а2, Р1 и Р2 приблизить второе слагаемое к нулю, что позволило бы пренебречь им, то после исключения аппроксимированные блоки будут трехдиагональными, если таковыми были соответствующие исходные матрицы, и будут рассчитываться по следующим формулам:

При исключении приграничных столбцов используются те же формулы, что и для внутренних, но при этом блоки, относящиеся к нулевому и (Ы+1)-му столбцу, полагаются равными нулю.

В идеале получаем четыре задачи на нахождение параметров:

Эти задачи в общем случае для Уф е Ят решения не имеют. Однако для многосеточного метода главным является уменьшение низкочастотной части невязки, следовательно, нам достаточно решить задачи (5) только для низкочастотных ф , что, в свою очередь, соответствует решению обобщенных спектральных задач о поиске наименьшего действительного собственного значения. Эти задачи можно решать одним из численных методов, например методом обратных итераций (см., например, работу [4]), однако, как оказывается, вполне достаточно использовать некоторое приближение к реальным собственным значениям.

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

Вышеописанный метод построения матриц на более грубых сетках для задач с симметричной матрицей будет корректным, поскольку матрица остатка в выражении (3) положительно полуопределена. Для задач с несимметричными матрицами вопрос о корректности описанного метода для модельных задач рассмотрен в работе [3], а в общем случае требует дополнительного исследования и в данной статье не рассматривается.

полненными. Перепишем АІ-1 следующим образом:

(3)

+

а1Цк-1 - в1 Як-1 + а1$1Бк -а1 Як - в2 Як-1 + а1в2 Бк -а2Цк-1 - РА + а2 РА -а2 Як - в2Цк + а2 в2 Бк

Бк-1 = Бк-1 а1Цк-1 в1 Як-1 + а1в1 Бк, К™ = а 1 Як +Р2 Як-1 -а 1Р2 Бк,

ЦТ = а2Ік-1 + РАк -а2РА,

Бп+1 = Бк+1 -а2Як - в2Цк +а2Р2Бк.

(4)

Як-1Ф = а1Бк Ф, Цк Ф = а2 Бк Ф, Цк-1Ф = РА Ф, Як Ф=Р2 Бк ф. (5)

Применение метода к задачам с симметричной матрицей

Рассмотрим задачу (1) с симметричной, положительно определенной блочно-трехдиагональной матрицей А из яп'тхп'т :

А = blocktridiag{-{ Ц, -Ц}, В{,Ц е Ятхт, I = 1,...,п, (6)

причем для матриц В{, Ц выполняются следующие условия:

В1 > 0; В1 = В?; Ц = I?. (7)

Эта задача является частным случаем описанной выше задачи. При этом заметим, что так как Як = Цк для любого к, то из выражений (5) получаем, что а1 = Р1, а2 = Р2 . Таким образом, для симметричной задачи вместо нахождения четырех параметров вычисляются всего два как приближения к наименьшим собственным значениям следующих задач:

1к-1Ф = а1Вк Ф, 1к Ф = а2 Вк Ф. (8)

При этом в соответствии с формулами (4) блоки матрицы будут пересчитываться по формулам

В^ = Вк-1 - 2а11к-1 +а 2 Вк,

■ЬП = 1^1 = а^ +а21к-1 -а1а2 Вк, (9)

ВП+г = Вк+1 - 2а21к + а2Вк.

В том случае, когда п = 21 -1, I е N, а а; = 0,5, формулы (9) будут

соответствовать варианту многосеточного метода с полуукрупнением описанному в работе [3], для скорости сходимости которого справедлива оценка (2).

Опишем другой вариант построения оператора крупносеточной коррекции, аналогичный негалеркинскому случаю из работы [3], который позволяет в случае, когда матрицы Ц диагональные, получать диагональные аппроксимации матриц на более грубых сетках Ц™, как мы увидим в дальнейшем, практически без изменения скорости сходимости.

Пусть для задач (6) —(7) матрицы на более грубых сетках строятся согласно описанному методу. Тогда в соответствии с формулами (9) изменяемые блоки матрицы будут выглядеть следующим образом:

' Вк-1 - 2а1Цк-1 + а1 Вк -(а1Цк + а2 Цк-1 -а1а2 Вк ) ^ ч-(а1Цк + а21к-1 - а1а2Вк ) Вк+1 - 2а21к + а2Вк у

Добавим к этому выражению положительно определенную симметричную матрицу Q следующего вида:

Г '

Q =

V 2 2 2 2 У

Заметим, что в соответствии с формулами (8) блоки этой матрицы будут близки к нулю на низкочастотных гармониках. При этом преобразовании

а

-2-(В - 1к-1 ) (Вк - 1к) -^(а1Вк - 1к-1 )(а2Вк - 1к)

-^(В - 1к-1 ) (а2Вк - 1к) 01 (а1Вк - 1к-1 ) (а2Вк - 1к)

а

а

недиагональные матрицы Бк исключаются с побочных диагоналей и блоки результирующей матрицы будут выглядеть следующим образом:

т^ПвШ ик-1

= Цк-1 2а 1 +0г) 1к-1 -0г 1к + а1(а 1 + а2)ик /

т пеш т пеш а2 т а1 т /1А\

Ьк-1 = Ьк+1 = Ьк-1 "у Ьк, (10)

т-лпеш ик+1

= ик+1 ^ Ьк-1 -(^2а 2 +0;1| 1к +а 2(а1 +а 2)ик •

Отметим, что при п = 21 -1, І є N и а, = 0,5 формулы (10) приводят к формулам второго варианта многосеточного метода с полуукрупнением, описанного в работе [3], для скорости сходимости которого вновь справедлива оценка (2).

Численные эксперименты

Рассмотрение численных примеров применения метода начнем с задачи Дирихле для уравнения Пуассона в единичном квадрате с однородными граничными условиями

-Дн = /, м|эо= 0, 0 = [0,1] х [0,1]. (11)

Разностная аппроксимация краевой задачи (11) приводит к системе линейных алгебраических уравнений вида (1) с матрицей вида

А = ЫоскМй1а2[-Ь, и, -Ь|, и, Ь є .

где и = tridiag{ -1,4, -1|, Ь = diag{1|•

Матрица данной задачи является симметричной и может также быть записана в виде звезды-шаблона [2]:

" 0 -1 0"

А = -1 4 -1

0 -10

Определим половину итерации блочного метода Гаусса — Зейделя как процесс, при котором значения неизвестных на нечетных (четных) столбцах находятся по значениям неизвестных на четных (нечетных) столбцах сетки. Таким образом, например, термин «полторы итерации» означает следующую последовательность половинных итераций: по исключаемым, по неисключаемым, затем опять по исключаемым столбцам. Преимущество такого подхода состоит в том, что, применяя завершающую итерацию по исключаемым столбцам, мы также обнуляем значение невязки на этих столбцах, что позволяет упростить пересчет правой части.

В таблице 1 приведены значения спектрального радиуса многосеточной итерации при различных значениях числа узлов N и различных значениях параметров аі (значение аі = 0,5 соответствует классическому варианту метода Галеркина для сужения матриц; аопт — параметры, находящиеся по формулам (12)):

а = (ь.-1ф- ф), а, = Ьф^. (12)

(икф, ф) ф)

Таблица 1

Зависимость спектрального радиуса многосеточной итерации от мелкости разбиения и выбора параметров аі

N 99 401 402 lll

ai 0,5 а ^опт. 0,5 а ^опт. 0,5 а ^опт. 0,5 а ^опт.

P 0,049 0,046 0,053 0,052 0,191 0,052 0,053 0,052

В формулах (12) тестовый вектор ф задавался как вектор с компонентами ф; = sin I — I.

IN J

В качестве сглаживающих итераций здесь и в дальнейшем использовались полторы итерации блочного метода Гаусса — Зейделя при пред- и послесглаживании. Норма, скорость уменьшения невязки, а также спектральный радиус оператора перехода многосеточной итерации определялись как

»N-1N-1 ...........

= I I r.I u. = \\rk\\/ rk-1 p = limo.,

;=i j=i1 1 ........... k^“

где rk — невязка на k-й итерации. Значение скорости сходимости p вычислялось в момент, когда ошибка уменьшилась не менее чем в

1010 раз.

Полученные значения подтверждают основное свойство многосеточного метода — сходимость практически не зависит от мелкости разбиения (числа уравнений). Кроме того, из приведенных результатов видно, что выбор параметров по формулам (12) позволяет избежать снижения скорости сходимости на некоторых сетках (в частности, на сетке с N = 402 и др.). Во всех дальнейших примерах параметров результаты приведены при значениях параметров, вычисленных по формулам (12). В таблице 2 даны результаты численного эксперимента для сравнения скоростей сходимости многосеточного метода в случае применения галеркинского и негалеркинского способа построения матриц на более грубых сетках.

Таблица 2

Спектральный радиус многосеточной итерации при использовании аналога метода Галеркина рх и негалеркинского матричного сужения Р2

N 99 25l 402 lll

Pi 0,046 0,051 0,052 0,052

P2 0,050 0,053 0,054 0,055

Как видно из таблицы 2, скорость сходимости многосеточного метода при применении негалеркинского способа построения матриц незначительно ниже классического метода Галеркина. Однако при использовании этого метода матрицы Ц на грубых сетках остаются диагональными, что позволяет снизить примерно в полтора раза количество операций для выполнения блочного метода Гаусса — Зейделя и упростить задачи (8) по нахождению параметров метода.

В таблице 3 показана зависимость спектрального радиуса многосеточной итерации от значения коэффициента е, мелкости разбиения и использованного метода матричного сужения для стационарного анизотропного уравнения диффузии:

еихх + иуу = /, и1за= 0, п = [ОД] х[0,1].

Таблица 3

Зависимость спектрального радиуса многосеточной итерации от значения коэффициента диффузии Б, мелкости разбиения и использования аналога метода Галеркина рг или негалеркинского матричного сужения р2

N 99 777

є 0,1 10 100 1000 0,1 10 100 1000

Р1 0,037 0,049 0,048 0,033 0,052 0,053 0,053 0,052

Р2 0,038 0,053 0,053 0,046 0,053 0,055 0,055 0,055

Матрица системы записывается в виде звезды-шаблона:

0 -1 0"

А = -є 2(1 + є) -є

_ о -1 0 _

Таким образом, мы видим, что скорость сходимости для рассмотренной задачи допускает оценку, не зависящую от выбора коэффициента є как при использовании аналога метода Галеркина, так и с использованием негалеркинского матричного сужения.

В завершение рассмотрим задачи для уравнений с переменными коэффициентами:

д_

ду

р(х, у) = /,

ду

Ш

= 0, 0 = [0,1] X [0,1].

Приведем результатні численных экспериментов по определению скорости сходимости предложенного варианта многосеточного метода для различных видов функции р(х, у).

В приведенных результатах экспериментов коэффициент є полагался равным 0,5, X — равным 10.

Таблица 4

Спектральный радиус многосеточной итерации для различных задач при использовании аналога метода Галеркина рг и негалеркинского матричного сужения р2

р( у) = ■

Функция p(x,y) N 99 257 402 777

p( x, y) 1 - xy P1 0,048 0,052 0,053 0,053

P2 0,051 0,054 0,054 0,055

1+ X sin [14 nx )sin (14 ny) P1 0,050 0,05 0,052 0,054

P2 0,061 0,061 0,062 0,063

= 1 + e(x(1 - x) + y(1 - у)) P1 0,048 0,052 0,052 0,053

P2 0,051 0,053 0,054 0,055

X, x e " 1.3" 4'4 , у e " 1.3" 4'4 P1 0,066 0,083 0,164 0,254

1, x € " 1.3" _ 4. 4 _ , у € " 1.3" _ 4'4 _ P2 0,058 0,067 0,069 0,069

Приведенные результатах показывают, что скорость сходимости описанного варианта многосеточного метода практически не зависит от размерности сетки и от коэффициентов уравнения (для широкого класса задач).

Работа выполнена при поддержке гранта РФФИ № 08-01-00431.

Список литературы

1. Ольшанский М. А. Лекции и упражнения по многосеточным методам. М., 2005.

2. Hackbusch W. Iterative solution of large sparse systems of equations. New York, Springer-Verlag, 1993.

3. Белякова О. В., Буздин А. А. Многосеточный метод с полуукрупнением для решения систем с блочной трехдиагональной матрицей // Методы вычислений. 2005. № 21. С. 5-19.

4. Калиткин Н. Н. Численные методы. М., 1978.

5. Hackbusch W. Multi-grid method and applications. Springer, 1985.

6. Шайдуров В. В. Многосеточные методы конечных элементов М., 1989.

Об авторах

А. А. Буздин — канд. физ.-мат. наук, доц., РГУ им. И. Канта.

С. С. Дедух — асп., РГУ им. И. Канта, e-mail: [email protected]

Authors

Dr A. A. Buzdin — assistant professor, IKSUR.

S. S. Dedukh — PhD student, IKSUR, e-mail: [email protected].