Многосеточный метод с шахматным исключением неизвестных без сглаживающих итераций Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 9, № 2, 2004

МНОГОСЕТОЧНЫЙ МЕТОД С ШАХМАТНЫМ ИСКЛЮЧЕНИЕМ НЕИЗВЕСТНЫХ БЕЗ СГЛАЖИВАЮЩИХ ИТЕРАЦИЙ*

В. В. Ефремов, В. В. ШАйдуров Институт вычислительного моделирования СО РАН,

Красноярск, Россия e-mail: shidurov@icm.krasn.ru

In paper the variant of V-cycle is considered with red-black elimination of unknowns to solve model grid analog of Poisson equation. The more effective operator is presented to project right-hand side on rough grid, which accelerates convergence of multigrid method. As a result, it is possible to dispense with smoothing iterations and nevertheless to attain constriction of error 5.7 times less by one V-cycle which has approximately the same number of algebraic operations as 5.5 simple iterations.

Введение

Начиная с основополагающих работ [1-3], стандартный многосеточный метод решения больших систем сеточных уравнений содержит сглаживающие итерации (предварительные и заключительные), проектирование системы на более крупную сетку и коррекцию решения на мелкой сетке с использованием поправки, интерполируемой с грубой сетки [4-6]. В данной работе предложен вариант многосеточного метода, в котором сглаживающие итерации отсутствуют, но сходимость по-прежнему имеется. Для решения сеточной аппроксимации уравнения Пуассона используется последовательность сеток, получаемая исключением каждого второго узла сетки с пятиточечными шаблонами-крестами, повернутыми на 45° на соседних уровнях. Известно, что отсутствие сглаживающих итераций в стандартном методе действительно не дает сходимости. В предлагаемой модификации сходимость обеспечивается более подходящей процедурой проектирования правой части на грубую сетку, не совпадающей с обычно рекомендуемыми процедурами [4-6]. Она выбрана на основании фурье-анализа, применяемого в ряде работ [4, 5] для количественной оценки сходимости многосеточных методов. Отметим, что обычно фурье-анализ не дает строгого обоснования сходимости многосеточного метода и является эвристическим для рекомендаций по выбору различных параметров, влияющих на сходимость. Для рассмотренного алгоритма вычислительный эксперимент дал хорошее совпадение скорости сходимости двухсеточного метода с ее предсказанием на основании фурье-анализа, в том числе наихудшего множителя подавления ошибки величиной 0.15 для рассмотренного алгоритма.

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант № 02-01-00523.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

Шахматное исключение приводит к уменьшению числа неизвестных в два раза и увеличению шага грубой сетки (равномерной квадратной, но повернутой на 45°) в \[2 раз. Более распространенный прием увеличения шага вдвое приводит к уменьшению числа неизвестных в 4 раза, но различие между операторами на грубой и мелкой сетке становится более существенным, что обычно приводит к увеличению числа итераций сглаживания между сетками. Здесь удалось совсем отказаться от сглаживающих итераций.

1. Постановка задачи

Пусть П = {г = (ж, у) : 0 < ж < 1, 0 < у < 1} - единичный квадрат с границей Г. Рассмотрим задачу Дирихле для уравнения Пуассона

^ ^ / в П, (1)

дж2 ду2

и = 0 на Г (2)

с непрерывной правой частью /(ж, у) £ С(П). Построим равномерную разностную сетку

П = = (ж^у,-) : Ж = ¿Л, г = 0,..., п; у,- = ^'Л, ] = 0, ...,п}

с шагом Л = 1 /п. Введем множество внутренних и граничных узлов П = П П П и Г^ = П П Г. Для дискретизации задачи (1), (2) используем стандартную разностную схему

= , ^ £ П; (3)

< = 0, ^ £ Г^ (4)

с пятиточечным шаблоном

(X и = —— + + +

Здесь и далее для произвольной функции ^(ж,у) принято обозначение = ^(ж^у^-). Известно, что задача (3), (4) невырождена и имеет единственное решение Ц", заданное на

П.

В дальнейшем систематически будет использоваться следующий прием исключения половины неизвестных. К уравнению с четной суммой индексов г + ] прибавим четыре соседних уравнения с нечетной суммой г + ] с целью исключения неизвестных с этими индексами и разделим левую и правую части уравнения на 2. В итоге получим сеточное равенство

(Ь^у = (М/ )у (5)

с операторами

/г. 3 к 1 к к к к

и ^ = (^2л)2 ^ - 2(^2Л)2 (и^-1'^-1 + иг—1,^+1 + иг+1,^-1 + иг+1,^-+1)-

(иг—2,^' + иг+2,^' + иг,^—2 + иг,^+2)

и

/. 1

(МЛ/ )у = ^ + + /¿-1,3 + /¿+ц +

Замечание 1. Для алгоритмической определенности в граничных узлах £ Г^ положим сеточное значение функции /¿3- равным нулю. С теоретической точки зрения это согласуется с нечетным продолжением сеточной функции пн через прямолинейные участки границы. Такой способ продолжения приводит к периодичности с периодом 2 по х и по у продолженной сеточной функции на

К1 = : г)3 = 0, ±1, ±2,...} .

В дальнейшем это делает оправданным применение фурье-анализа.

Для упрощения описания сеточных операторов мы будем далее использовать схематическую запись их шаблонов. Например, для V1 шаблон имеет вид

х

<

-1

-2 0 -2 -1 0 12 0 -1 -2 0 -2 -1

< >

г - 2 г - 1 г г + 1 г + 2

3 + 2 3 + 1 } 3

3-1

3 - 2

Для исключения новых появившихся неизвестных в углах шаблона V1 используем сеточный оператор

(БН<р)ч = ^ (<#-1>3 + &+1,3 - &,з-1 - ^+1).

На гладких функциях (р он аппроксимирует выражение

02<р 02<р

дх2 ду2

со вторым порядком аппроксимации. Используя приближенное равенство Бнпн ~ 0 четыре раза в узлах и приходим к приближенному равенству

(Ь^п% « (Мн/)з (6)

с оператором шаблон которого имеет вид

1 ( -1 0 -1 ^— х < 0 4 0

(л/2^)2 \ -10 -1

Итак, вновь получено разностное уравнение, на этот раз с пятиточечным шаблоном "наклонный крест". Это сеточное уравнение аппроксимирует исходное уравнение (1) также со вторым порядком аппроксимации, что легко проверить путем разложения функции в ряд Тейлора относительно узла . Применение оператора Бн не нарушило порядка аппроксимации ввиду симметрии коэффициентов, давших добавку (Бн(Бнпн))з, которая на гладких функциях имеет вид

3 \дх2 ду2) \дх2 ду2

1

Обозначим

Q^2h = {zij е Qh : i + j — четное},

Q^2h = Q^2h П Q, r^2h = Q^2h П Г-

Взяв равенство (6) в узлах Q^h как уравнение для определения новой сеточной функции uV^h, добавим краевое условие в узлах Г^^. В итоге получается невырожденная система

(L^V52^ = (Mhf на Q^, (7)

= 0 на Г^-

Эта система имеет единственное решение U^22, заданное на Q^h и, вообще говоря, не совпадающие с uh.

Это исключение неизвестных в шахматном порядке в зарубежной литературе носит название красно-черного исключения [4, 6] и является основой ряда многосеточных методов. Сформулируем один из них в двухсеточном варианте. Многосеточный метод получается рекуррентным применением двухсеточного.

Итак, пусть задано приближение vh к решению uh на сетке Qh, совпадающее с ним на границе Г2. Тогда процедура Increase (h, vh, fh,v^) уточнения решения vh ^ v^ системы (3), (4) на сетке h с правой частью fh реализуется в четыре шага.

1. Вычисляем невязку

rh = fh - Lhv h на Qh- (8)

2. Решаем систему

= M hrh на Q^h, (9)

= 0 на Г^-

3. Полученную поправку добавим к приближенному решению

$ = vh + на (10)

4. По приближенному решению в узлах Q^h построим новое приближенное решение Vi в узлах Qh\Q^2h путем решения уравнений

Lhv2h = fh на Qh \ Qyih,

v2h = 0 на rh \ r^^h, (11)

v22 = на Q^2h-

Отметим, что каждое уравнение (11) содержит только по одному неизвестному значению v^. В итоге получим новое приближение v^.

Повторяя шаги 1-4, получаем итерационный цикл, в ряде случаев сходящийся с множителем, не зависящим от шага h.

2. Фурье-анализ

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

е^ = ^ - ^

разлагается в ряд Фурье

емз = Е ¿(01,02)ехр(1(г01 + ¿02)). (12)

01,02

Здесь 1 — мнимая единица, а сумма берется по конечному набору параметров

п к п /

01 = ± —, к = 1, ..., П - 1 и 02 = ± —, / =1,...,/ - 1, п п

о которых нам достаточно знать, что

-п < 01 < п, -п < 02 < п. Вычисление невязки эквивалентно следующей операции:

тл = ГЛе£. (13)

Применительно к разложению (12) это дает равенство

1

01,02

= ^ (4 - 20С801 - 20С8 02^ А(01,02)ехр(1(г01 + ^02>).

(14)

Для упрощения записи введем коэффициент-символ оператора V", действующего на гармонику с параметрами 01, 02:

^(01,02) = А(4 - 2 сое01 - 2сое 02). (15)

п2

Вычисление правой части в (9) связано с символом оператора М^:

МЛ(01,02) = 2 + 4 СС8 01 + 4сс8 02. (16)

Оператор Г1^ умножает гармонику с параметрами 01, 02 в четных узлах на следующий множитель:

Ь^(01, 02) = (4 - 2 СС8(01 + 02) - 2со8(-01 + 02)). (17)

Как уже говорилось, оператор, обратный к Г1^, существует и невырожден, так что он производит деление гармоники на этот же множитель. В итоге гармоники ошибки е^ с параметрами (01,02) в четных узлах приобретают следующий вид:

Е (02) М(01, 02)^(01, 02^ А(01, 02) ехр(1(г01 + ¿02)). (18)

= ^^ ' _1_МЬ(0. 0. ) г^

Обозначим

WlA) = -Mh(0!,02)Lh(^l,02).

Добавление поправки в четных узлах приводит к равенству

eh,i,j = X] (l - D(0iA)) A(0iA)exp(i(i0i + j^))-

Введем функцию

/2 (0 0) в четных узлах,

j( 1, 2) | (cos 0i + cos 02)/2 в нечетных узлах.

Тогда на шаге 4 получаем в целом по сеточной области

= J] Ih(0i, 02Д1 - D(0i, 02) A(0i, 02) exp(i(i0i + j02)).

h

01,02

Поставим целью избавиться в /2 от зависимости i, j. Для этого представим /2 (0i, 02) в виде суммы двух значений других индексов Jh(0i, 02) и Jh(0i + п, 02 + п). А именно, гармонику с параметрами (0i, 02) с амплитудой A разложим в суммы гармоники с параметрами (0i, 02), (0i + п, 02 + п) и амплитудами Ai, A2. Для этого найдем такую функцию Jh(0i, 02), что

Ai = Jh(0i,02)A, A2 = Jh(0i + п, 02 + п)А

и

Ai exp (i(i0i + j02^ + A2 exp (i(i(0i + п) + j(02 + п))) = = Ij(0i,02)Aexp (i(i0i + j02)).

В четных узлах имеем Ai + A2 = A, а в нечетных узлах Ai — A2 = A(cos 0i + cos 02)/2. В итоге получаем систему

Jh(0i,02) + J2(0i + п,02 + п) = 1,

Jh(0i, 02) - Jh(0i + п, 02 + п) = (cos 0i + cos 02)/2. ( )

Ее решением является функция

Jh(0i,02) = 1 + 4cos 0i + 1 cos 02.

Итак, гармоника с параметрами (0i, 02) оказалась неинвариантной, она породила дополнительную гармонику с параметрами (0i + п, 02 + п). Аналогично гармоника с параметрами (0i + п, 02 + п) порождает дополнительную гармонику с параметрами (0i, 02) (в силу периодичности с периодом 2п гармоники с параметрами (0i + 2п, 02 + 2п) и (0i, 02) совпадают). Поэтому рассмотрим взаимодействующую пару гармоник, она оказывается инвариантной [4, 6]. Символ оператора перехода для этой пары является матрицей 2 х 2:

( Jh(0i,02)(1 - D(0i,02)) Jh(0i,02)(1 - D(0i + п, 02 + п)) \( Ai N =

I Jh(0i + п, 02 + п)(1 - D(0i, 02)) Jh(0i + п, 02 + п)(1 - D(0i + п, 02 + п)) i I A J

Jh(01,02КAl(1 - D(01,02)) + A2(1 - D(01 + П, 02 + п))) Jh(01 + п, 02 + пXA^ 1 - D(01,02)) + 1 - D(01 + П, 02 + п)))

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

S1 = exp(i(i01 + j02)) (A1 (1 - D(01, 02)) + A2(1 - D(01 + п, 02 + п))), а в нечетных узлах

cos 01 + cos 02

S = exp(i(i01 + j02))

2

(A1 (1 - D(01, 02)) + A2(1 - D(01 + п, 02 + п)))

Эти выражения можно оценить сверху:

5 = тах{|^1, |^2|}< А ||(1 - £(01,02))| + А1|(1 - £(01 + п,02 + п))| < (20)

< тах{|А11, |А21}(|(1 - £(01,02))| + |(1 - £(01 + п,02 + п))|-.

Отсюда видно, что чем меньше 11 - £ (01, 02) |, 11 - £ (01 + п, 02 + п) |, тем меньше ошибка. На рис. 1 представлен график функции |1 - £(01, 02)|. График функции |1 - £(01 + п, 02 + п)| совпадает с ним ввиду симметрии относительно двух осей координат.

Для 01 + 02 = 0 и 01 - 02 = 0 значение 1 - £(01,02) равно нулю. Значение функции 1 - £(01, 02) для параметров (0, ±п) и (±п, 0) не так мало, что обычно приводит к дополнительным сглаживающим итерациям на мелкой сетке для исключения гармоник ошибки с такими параметрами.

Попытаемся найти другой оператор Мк, такой, что выражение

1

(01,02;

Mh (01,02 )Lh (01,02 )

будет более близким к 1, чем £ (01, 02). В результате более тщательного подбора и тестирования рекомендуется более эффективный оператор сглаживания со следующим индексом:

М^ (01, 02) = Мн (01, 02) + ^(1 - СС8(01 + 02 ))(1 - СС8(-01 + 02)).

о

Рис. 1. График функции |1 - 0(01,02)|, максимум равен 0.5 в серединах сторон.

График функции

1 - В (0! ,02)

1 -

Ь^ (0! ,02 )

мк (0! ,02 )ЬЬ (0! ,02 )

изображен на рис. 2.

Этому индексу МЬ (0! ,02) соответствует следующий шаблон оператора МЬ:

32

х <

1

-2 4 -2

4 20 4

-2 4 -2

1

(21)

Замечание 2. В соответствии с^нечетным продолжением сеточных функций через границу (см. замечание 1) оператор Мк в приграничных узлах можно трансформировать так, что будут использованы лишь известные значения. Например, Д п = 0, а Дп+1 = - / п-!• Поэтому шаблон Мк в узле 2г п-! можно записать в виде

32

х

1 4 19 4 1 1 п - 1

-2 4 -2 п - 2

1 п - 3

/

г - 2 г- 1 г г + 1 г + 2

В углах он трансформируется еще сильнее, например:

1 4 18 1 п- 1

1

32

х

-2 4  п - 2 1 I п - 3 п — 3 п — 2 п — 1

1

1

1

Рис. 2. График функции |1 - Б(6\,02)|, максимум равен 0.12 в узлах (п?/2,пк/2), где ] + к — нечетно.

3. Вычислительный эксперимент

Рассмотрим задачу (3)-(4). Положим п = 32, Л, =1/32 и проведем вычислительный эксперимент с заданной сеточной функцией

= вт(пгЛ,г) 8т(п^в), г,^ = 0,1,...,п. (22)

Она является сеточной модой Фурье с частотами (пЛг, пЛ,в), которые мы будем брать для различных целых г, в = 1, 2, ...,п - 1. Правая часть /^ уравнения (3) вычисляется прямой подстановкой решения в оператор Для оценки точности решения будет использоваться дискретная Ь2-норма

\

31

,h )2

h2 ^ (uhj )2

В качестве начального приближения г^1 возьмем нулевую сеточную функцию. После этого реализуем процедуру Increase (h,vh, fh,vh) по формулам (8)-(11). В итоге мы получим новое приближение г^, для которого вычислим отношение

h = l~h-ТцГ1 • (23)

Значения этой величины являются коэффициентом подавления ошибки в результате применения одного V-цикла в двухсеточной реализации. Ее значения для различных r и s приведены в табл. 1. Таблица 2 содержит верхнюю оценку этой величины

|(1 - D(0 1, 02))| + |(1 - D(0 1 + п, 02 + п))|,

вычисленную по методике фурье-анализа (12) - (19).

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

Таблица!

r \ s 1 10 16 22 31

1 0 0.1439 0.2981 0.3632 0.3584

10 0.1439 0 0.0622 0.1638 0.3632

16 0.2981 0.0622 0 0.0622 0.2981

22 0.3632 0.1638 0.0622 0 0.1439

31 0.3584 0.3632 0.2981 0.1439 0

Таблица2

r\s 1 10 16 22 31

1 0 0.2162 0.4952 0.7744 0.9952

10 0.2162 0 0.1544 0.4718 0.7744

16 0.4952 0.1544 0 0.1544 0.4952

22 0.7744 0.4718 0.1544 0 0.2162

31 0.9952 0.7744 0.4952 0.2162 0

ТаблицаЗ

r\s 1 10 16 22 31

1 0 0.1116 0.1483 0.0798 0

10 0.1116 0 0.0173 0 0.0798

16 0.1483 0.0173 0 0.0173 0.1483

22 0.0798 0 0.0173 0 0.1116

31 0 0.0798 0.1483 0.1116 0

Таблица4

r\s 1 10 16 22 31

1 0 0.1676 0.2464 0.1702 0

10 0.1676 0 0.0428 0 0.1702

16 0.2464 0.0428 0 0.0428 0.2464

22 0.1702 0 0.0428 0 0.1676

31 0 0.1702 0.2464 0.1676 0

После этого вычисления были повторены с помощью процедуры Increase, в которой оператор Mh был заменен на оператор Mh в формуле (9). Таблица 3 содержит значения величины &Mh, вычисленной по аналогии с (23), а табл. 4 — теоретическую верхнюю оценку этой величины, определенную по методу фурье-анализа.

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

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

Для многосеточного метода рекуррентный учет ошибки в сеточной Ь2-норме подсчи-тывается путем оценки суммы [6]

q + q2 + ... + qm-1,

где m — число сеток. В нашем случае имеем

q + q2 + ... + qm-1 < q + q2 + ... = q/(1 - q).

Таким образом, множитель подавления ошибки в многосеточном V-цикле для этого метода не хуже 0.1764.

Подсчет арифметических операций показывает, что их количество пропорционально числу неизвестных на самой мелкой сетке с множителем 35, что меньше числа операций для пяти c половиной простых итераций на той же сетке.

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

[1] бахвалов н.с. О сходимости одного релаксационного метода при естественных ограничениях на эллиптический оператор // ЖВМ и МФ. 1966. Т. 6, № 5. С. 861-883.

[2] федоренко р.п. Релаксационный метод решения разностных эллиптических уравнений // ЖВМ и МФ. 1961. Т. 1, № 5. C. 922-927.

[3] федоренко р.п. О скорости сходимости одного итерационного процесса // ЖВМ и МФ. 1964. Т. 4, № 5. C. 559-564.

[4] Briggs w.l., Henson v.e., McCormick s.f. A Multigrid Tutorial. Second Edition. Philadelphia: SIAM, 1999.

[5] Diskin b. Multigrid solvers on decomposed domains. Thesis for the M.Sc. degree, Department of Applied Mathematics and Computer Science. The Weizmann Institute of Science, 1993.

[6] шайдуров в.в. Многосеточные методы конечных элементов. М.: Наука, 1989.

[7] Hackbusch w. Multi-Grid Methods and Applications. Berlin: Springer, 1985.

[8] ольшанский м.а. Лекции и упражнения по многосеточным методам. М.: Изд-во ЦПИ при механико-мат. фак. МГУ, 2003.

Поступила в редакцию 5 февраля 2004 г., в переработанном виде — 27 февраля 2004 г.