МАТЕМАТИКА
УДК 519.615.5
А. А. Буздин, Е. А. Васильева
ОПТИМАЛЬНЫЕ ПАРАМЕТРЫ ДЛЯ ПОСЛЕДОВАТЕЛЬНОСТЕЙ КАСАТЕЛЬНЫХ И ДВУХЧАСТОТНЫХ РАЗЛОЖЕНИЙ
Рассматривается метод построения предобусловливателя для решения больших систем линейных уравнений, основанный на неполном блочном разложении блочных трехдиагональных матриц. Метод является обобщением методов, изложенных в [1; 2]. На основе анализа ряда модельных задач был предложен способ выбора квазиоптимальных параметров для последовательностей касательных и двухчастотных разложений, обеспечивающий скорость сходимости лучше, чем метод переменных направлений в коммутативном случае.
A preconditioner for large systems of linear equations based on an incomplete block-decomposition for a block-tridiagonal matrices is considered. This method generalizes methods are developed by Buzdin and Wittum [1; 2]. A method of choosing quasioptimal parameters for consequences of tangential and two-frequencies decompositions based on analysis of model problems is developed. It's convergence rate is better than ADI's convergence rate in commutational case.
Ключевые слова: системы линейных уравнений, неполные блочные разложения, касательное
разложение, двухчастотное разложение, оптимальные параметры.
Key words: systems of linear equations, incomplete block-decomposition, tangential decomposition, two-frequencies decomposition, optimal parameters.
Введение
Пусть блочная матрица К подлежащей решению системы уравнений Ки = F, К е Rnmxnm, и, F е Rnm имеет вид K = blocktridiag{ Lj-i, D;, Lj}.
Полное блочное разложение для матрицы К имеет вид: К = (L + T) T -1 (U + T), где L и U —
соответственно нижняя и верхняя блочно-треугольные матрицы, а T — блочнодиагональная матрица
T = blockdiag{Ti, ..., Tm}, причем
Ti = Di, Tj = Dj - Lj-i ТД Lj-i, j > 1. (1)
Если заменить Tj их некоторой аппроксимацией Tj, то получим неполное блочное разложение
матрицы K. Для модельной задачи, когда в матрице К все блоки Dj = C, а Lj = - I, блоки Tj = fj(C), где fj — некоторые рациональные функции. Приближая их линейными функциями от матрицы C, придем, в зависимости от способа построения приближения, к касательному или двухчастотному разложению. Двухчастотное разложение матрицы K получим, если рассмотрим секущую, проходящую через точки (а(1), fj (аЯ)), (а(2), fj (а(2))) кривой f(k) и сопоставим ей
линейную функцию матрицы C. К определению касательного разложения матрицы К придем, если вместо секущей рассмотрим касательную в точке (A*, fj(A*)). Нетрудно проверить (см. [1; 2]), что блоки Tj могут быть рассчитаны по следующим рекуррентным формулам:
~i = с, т = с+j j т-1 - (j + j) i, j > i,
где j = 1/fj (A(0); j 1. Такое представление матриц Tj позволяет обобщить определение
двухчастотного разложения на матрицы вида K = blocktridiag{Lj-i, Dj, Lj}.
Определение 1. Пусть K = blocktridiag{Lj-i, Dj, Lj} - симметричная блочная трехдиагональная
матрица, ег, ei е Rm - два тестовых вектора. Тогда для матрицы К можно определить двухчастотное разложение
M = (L + ~) T-1(LT + ~),
где T1 = D1, ~ = Dj + j j ~-1 - (j + j) Lj_1, j > 1,
j = (Lj-1ei, ei)/(~jei/ e/)z l - 1,2
Определение касательного разложения, отвечающее тестовому вектору е, получается из определения 1 при ц(2) ^ ц(1) = (L j-1e, e)/(Tj.e, e).
Результатах теоретического анализа и проведенных численных экспериментов показали, что при использовании этих разложений в качестве предобусловливателей в методе сопряженных градиентов для модельной задачи скорость сходимости имеет порядок 1 - O( к173). Численные эксперименты показывают аналогичное поведение скорости сходимости и для более широкого класса задач.
Дальнейшее увеличение скорости сходимости достигается применением набора касательных или двухчастотных разложений (см.: [1—4]). В [3 и 4] получены качественные оценки скорости сходимости такого метода. В [1] доказано, что для модельной задачи этот способ обеспечивает скорость сходимости не худшую, чем метод переменных направлений в коммутативном случае. В отличие от этих работ, здесь основное внимание уделено вопросу о выборе оптимальных параметров последовательностей касательных и двухчастотных разложений.
Задача о выборе оптимальных параметров
Решению задачи о выборе оптимальных параметров для последовательностей касательных и двухчастотных разложений посвящены работы [1 — 7] и некоторые другие.
Пусть Б; — итерационный оператор, отвечающий касательному или двухчастотному
разложению, Бп = ^ — итерационный оператор последовательности касательных или
двухчастотных разложений. Чтобы найти оптимальные параметры последовательностей
касательных разложений Х(г) (I = 1, ..., к), найдем такие их значения, чтобы
IIS
nil к
ntlSi
(2)
Аналогично для последовательностей двухчастотных разложений требуется определить 2к параметров Х(21 -1), х(2г)(I = 1, ..., к) с такими же свойствами. Используя результаты работы [2], можно показать, что для нормы итерационного оператора для последовательности касательных и двухчастотных разложений в случае модельной задачи справедливы соответственно оценки
Ц^пЦк ^Skoc(v, v(1), ..., v(fc))= max
ve[vmn.v
min ' max J l =
k (
Л
] l=1
v - V
(l)
f» (tf0)
v + v
(l)
IsA^v, v(1),v(M)) = max
VG[v • , V 1
L vmm' maxJ
(v-v^Xv-v™)
.(21)4
П77______________ ______
l-1 v+^/V2^
(2l)
(З)
соответственно
где V - X - 2; V max -Xmax - 2; Vmin -Xmin - 2; fx (X)-X/2 ^X2 /4 - 1; а Xmin, Xma минимальное и максимальное собственные значения матрицы C. Поэтому для модельной задачи задачу поиска оптимальных параметров последовательностей касательных и двухчастотных разложений можно свести к минимизации функций Sкас(v, v(1), ..., v(k)) и Бдв(v, v(1), ..., v(2k)).
Так как для того чтобы при любом числе блоков матрица системы K была положительно определена, нужно, чтобы хттЗ:2, а при Х^2 функция ^(X^l, то мы можем получить
приближенное решение задачи, если в (3) и (4) заменим (X(1)) на единицу. Тогда задача минимизации существенно упрощается. В этом случае для последовательностей касательных разложений минимизируемая функция примет вид
Sа (V, V(1).......V(k)) = П
l=1
vv+v(l),
К
2
2
2
Минимизационная задача с функцией вида (5) возникает при исследовании метода переменных направлений. Ее решение подробно представлено в [8]. В [4] показано, что скорость сходимости не будет зависеть от шага сетки, если число разложений растет так же, как логарифм числа неизвестных в блоке матрицы системы K. Подобная зависимость числа разложений от размера блоков, как было показано в [5], наблюдается и для последовательностей двухчастотных разложений.
В данной работе под минимизационной задачей понимается задача минимизации функций Sкас(v, v(1),..., v(k)) и Sde(v, V(1), ..., v(2k)). В отличие от вышеупомянутых работ, исследования были проведены без замены функции /»(X) на 1. Аналитическое решение этих задач вряд ли возможно и для ее решения использовались численные алгоритмы.
Для параметров v(l), являющихся точными или приближенными решениями этих задач, мы сохраним название оптимальных, несмотря на то, что правильнее их было бы называть квазиоптимальными.
Численный алгоритм поиска оптимальных параметров
Фактически решение сформулированных минимизационных задач сводится к построению функций, наименее уклоняющихся от нуля, причем интересующие нас параметры v(l) являются нулями этих функций. Другими словами, для выбора оптимальных параметров касательных разложений требуется среди множества функций Sкас(v, v(1),..., v(k)), однозначно определяемых своими параметрами v(l), выбрать такую, чтобы
РиС = max \Skuc(v/ у(1)/•••/ v^l = min max \Skuc(v/ v(1)/•••/ v(k>>|.
ve[vmln/vmax1 1 (v(1>,...,v(k>) ve[vmln/vmax] 1
Функция SKac (v, v *1),v *k)) характеризуется свойством: число последовательных точек отрезка [vmin, Vmax], в которых она принимает с чередующимися знаками экстремальное значение р® , не меньше k + І.
Аналогично можно сформулировать задачу выбора оптимальных параметров для последовательностей двухчастотных разложений и функцию Sde(v, v(1),..., v(2k)), для которой число последовательных точек [vmin, Vmax], в которых она принимает с чередующимися знаками экстремальное значение р®, не меньше 2k + І.
Сформулируем алгоритм, который в дальнейшем мы будем называть «алгоритмом выравнивания экстремумов» (в дальнейшем, для краткости, алгоритм ВЭ), позволяющий найти k параметров v(l) для последовательностей касательных и двухчастотных разложений для различных видов функции S(v, v(1),..., v(k)), имеющих k корней. С его помощью можно определить точки v(l) є [vmin, vmax] такие, чтобы выполнялось условие
|S(vmrn' v<1) / ./ VW | = max |S(V, V(1), V(k))| = |S(Vmax, v(1),v^, (6)
I I ve[v(l), v(l+1)^ 1 1 1
где l = І, ..., k - І.
Работа алгоритма основана на следующем свойстве функций S(v, v(1), ..., v(k)): если мы
двигаем параметр v(lo) влево (уменьшаем) при фиксированных остальных значениях параметров v(l), то ее значение в точке v(l°-1) и значения во всех точках v(l), лежащих левее ее, уменьшаются, а в точках, лежащих правее, — увеличиваются. Обратная картина наблюдается при увеличении значения v(l°) (движении вправо). Поэтому если | El | > j El+i |, то для выравнивания экстремумов параметр v(l) сдвигается влево, в противном случае — вправо. Многие численные исследования показали надежность предложенного алгоритма и его высокую работоспособность.
Алгоритм выравнивания экстремумов
Для нахождения к параметров У(г) е [у^д, Утах] функции Б(у, у(1), ..., У(к)) таких, чтобы было истинным условие (6), последовательно выполняются следующие действия.
1. Задаются произвольным образом к параметров У(г).
2. Ищется к - 1 локальный экстремум функции Б(у, ..., V®) Е, = тах р(у, у(1), V®),
,„-Г„(/) „(1+1)1І I
Vф('), /+1)]1 0 :
(I = 1, к - 1) и вычисляются значения этой функции на концах отрезка Е0 = |5(ут1п, ..., у®|,
Ек = |^(утах, у(1),... , у(к)|.
3. Среди значений Ег (I = 0, ..., к + 1) ищутся два соседних Е^ и Ег +1, разность между абсолютными величинами которых максимальна:
Е, - Е
Мо
10 + 1
=таХ |е|- Iе/+1 1=^-
(7)
Х/о +1)
Определяются величины I = у(г° ) и Г = У(
4. Параметр у(го) двигается влево или вправо, в зависимости от знака разности |е^ | - |е^ + ^ до
тех пор, пока величина ДЕ не уменьшится в заданное число раз.
5. Последовательно повторяются пункты 3 — 4 до тех пор, пока с заданной точностью не выполнится условие (6).
Отметим, что поиск экстремумов можно проводить, например, по методу «золотого сечения»
(см., напр., [9]), а движение параметра
Л)
методом деления отрезка пополам. При этом удобнее
использовать значения не самих параметров V параметры ^г) следующим образом:
(/)
ДО -
= 48ІП
2 ЛЮ
,(')
2(п +1)
а параметров ю(/) , однозначно определяющих
(8)
Работа алгоритма ВЭ тестировалась на задаче (5), при к = 2р имеющей простое аналитическое решение [10]. Во всех случаях наблюдалось полное совпадение полученных значений с теоретическими.
Проиллюстрируем работу алгоритма ВЭ на примере функции 5кас(V, ..., Vк)) из задачи (5)
при помощи рисунка. На нем отображен процесс «выравнивания» экстремумов, описанный в алгоритме ВЭ. На графике представлена зависимость значения максимальной разности между локальными экстремумами ДЕ от номера шага алгоритма ] для последовательности девяти касательных разложений с шагом сетки к = 1/512. Под «шагом алгоритма» понималась совокупность действий, описанная в пунктах 3—4 алгоритма ВЭ, приводящая к уменьшению максимальной из разностей соседних экстремумов ДЕ в 10 раз. Параметры ю(г) (I = 1, ..., к) задавались как ю(г) = I + 1, а ю(0) и ю(к+1) полагались равными 1 и п + 1 соответственно. Критерием выхода являлось выполнение условия (6) с точностью 10-8 .
Рис. Иллюстрация алгоритма «выравнивания экстремумов»
Этот алгоритм использовался для получения значений теоретических оценок скорости сходимости и параметров разложений во всех рассмотренных в дальнейшем задачах.
Результаты исследования и их анализ
Для того чтобы иметь возможность сравнивать эффективность использования последовательностей различного числа касательных или двухчастотных разложений, скорости сходимости всех последовательностей следует пересчитать в расчете на одну итерацию соответствующего разложения, так как для реализации последовательности к разложений требуется примерно в к раз больше арифметических действий, чем для реализации одного разложения. Введем понятие «эффективной скорости сходимости», которое мы будем в дальнейшем широко использовать.
Определение 2. Эффективная скорость сходимости составной итерации % из к касательных или
двухчастотных разложений - это ф|8П||К .
В работе [1] было показано, что скорость рассматриваемого метода не хуже, чем в методе переменных направлений при оптимальном выборе параметров. Использование алгоритма ВЭ позволяет дать количественную оценку эффективности метода.
В первой строке таблицы 1 через %П'н обозначены значения теоретической оценки эффективной скорости сходимости для метода переменных направлений для различного числа касательных разложений к = 1, ..., 8 при шаге сетки к = 1/512, а во второй через %(2)^3) обозначены значения оценки эффективной скорости сходимости для той же задачи, полученные при использовании параметров решения минимизационной задачи (2) — (3).
Таблица 1
Оценки эффективной скорости сходимости для различного числа касательных разложений для уравнения Пуассона
к 1 2 4 8 16
пГ 0,988 0,855 0,659 0,555 0,509
(2)-(3) Л1 0,917 0,747 0,573 0,512 0,449
В таблице 2 для последовательности восьми касательных разложений представлены теоретические и численные значения скорости сходимости для задачи Дирихле для уравнения Пуассона с однородными граничными условиями при различном размере задачи. Теоретические значения оценки для нормы итерационного оператора ^теор получены численным решением минимизационной задачи (2) — (3). Для сравнения во второй и третьей строках приведены значения скорости сходимости Лчисл и эффективной скорости сходимости Л1тасл, полученные для
реальной двумерной задачи при выборе таких же параметров ю(г), что и при получении оценки
Лтеор*
Таблица 2
Теоретические и численные значения скорости сходимости последовательности восьми касательных разложений для уравнения Пуассона
Н 1/16 1/32 1/64 1/128 1/256 1/512 1/1024
Лтеор 6,00е-8 3,18е-6 4,43е-5 2,80е-4 1,08е-3 3,01е-3 6,70е-3
Лчисл 2,07е-9 7,87е-7 2,14е-5 1,74е-4 7,31е-4 2,25е-3 5,35е-3
Л1числ 0,08 0,17 0,26 0,34 0,41 0,47 0,52
) 3,45е-4 3,24е-4 6,98е-4 1,04е-3 1,40е-3 2,00е-3 2,37е-3
Лі(ю«) 0,13 0,20 0,30 0,37 0,43 0,49 0,54
Из таблицы 2 видно, что значения скорости сходимости %теор и %числ близки, причем они становятся ближе с ростом размера задачи.
Вызывает интерес следующий вопрос: существует ли более простой способ задания параметров ю(г) последовательностей касательных разложений для уравнения Пуассона, обеспечивающий высокую скорость сходимости. Анализы значений оптимальных параметров ю(г) последовательностей к касательных разложений, численно рассчитаных для задачи, где число
і пк
разложении к и размер задачи п при этом связаны соотношением п = 2 , показывает, что они с высокой степенью точности (за исключением, быть может, последних) составляют геометрическую прогрессию со знаменателем ~ 2. Таким образом, можно предположить, что, выбрав для задачи размером п = 2к оптимальные параметры следующим простым способом
ю(г) = 2г_1 (I = 1, ..., к), (9)
можно получить высокую скорость сходимости.
В таблице 2 в последних двух строках представлены значения скорости сходимости п(ю«) и эффективной скорости сходимости п1(ю«) последовательностей к касательных разложений для
задач с шагом к = 1/ 2к при приближенном выборе параметров (9), полученные за 30 итераций. При таком выборе параметров вплоть до сетки 1024 х 1024 значения эффективной скорости сходимости не превосходят 0,54.
В таблице 3 представлены теоретические и экспериментальные значения оценки для нормы итерационного оператора для последовательности четырех двухчастотных разложений. В ней через Лтеор обозначены значения теоретической оценки скорости сходимости, полученные численным решением минимизационной задачи (2), (4). Через Пчисл и ПітаСл обозначены экспериментальные значения скорости сходимости и эффективной скорости сходимости, полученные при тех же значениях параметров ю(г).
Таблица 3
Оптимальные параметры и численные значения скорости сходимости последовательности четырех двухчастотных разложений
Н 1/16 1/32 1/64 1/128 1/256 1/512 1/1024
Лтеор 1,22е-4 1,17е-3 5,15е-3 1,44е-2 3,06е-2 5,43е-2 8,50е-2
Лчисл 1,13е-4 1,14е-3 5,09е-3 1,43е-2 2,47е-2 4,68е-2 7,46е-2
Л1числ 0,01 0,17 0,25 0,33 0,40 0,46 0,52
п(ю*) 3,19е-5 1,48е-4 4,96е-4 9,58е-4 1,32е-3 1,53е-3 1,65е-3
Л1(ю*) 0,073 0,17 0,28 0,37 0,43 0,48 0,53
Сравнение результатов показывает, что значения теоретической оценки с высокой степенью точности совпадают со значениями скорости сходимости, полученными экспериментально.
Предложим простой способ выбора параметров разложений последовательностей двухчастотных разложений, обеспечивающий высокую скорость сходимости. Анализ численных результатов показывает, что значения оптимальных параметров ю(г) (I = 1, ..., 2к) для
последовательностей к двухчастотных разложений для задач с шагом п = 2к можно приближенно выбирать следующим образом:
ю(2'-1) = 21-1, ю(2г) = [1,5о/2'-1)] (I = 1, ..., к), (10)
где [ •] обозначает округление числа до ближайшего целого.
Этот вывод подтверждается результатами, представленными в последних двух строках таблицы 3, в которых через %(ю~) и %1(ю~) обозначены значения скорости сходимости, полученные экспериментально при выборе параметров способом из (10).
Анализ численных значений эффективной скорости сходимости для последовательностей двухчастотных разложений показывает, что она несколько выше, чем у последовательностей касательных разложений. Однако на реализацию одного двухчастотного разложения затрачивается несколько большее число действий, чем на реализацию касательного. Поэтому однозначного предпочтения ни одному из видов последовательностей отдать нельзя. Отметим, как и в [2], более слабую зависимость скорости сходимости от параметров у последовательностей двухчастотных разложений при применении метода сопряженных градиентов. При небольшом числе разложений (к = 1, ..., 4) эффективная скорость сходимости резко уменьшается, а затем
практически перестает изменяться. Поэтому для задач среднего размера (h и 1/200) можно рекомендовать использование небольшого числа разложений (два-четыре) в сочетании с методом сопряженных градиентов.
В работе [6] вопрос приближенного выбора оптимальных параметров последовательностей касательных и двухчастотных разложений рассматривался на примере анизотропного уравнения диффузии єиxx + Uyy = f в зависимости от значения є, смешанных краевых задач для уравнения Пуассона и стационарного уравнения диффузии с разрывными коэффициентами. Сформулируем некоторые выводы, которые можно сделать, анализируя полученные результаты.
1) Выбор параметров (9) и (10) эффективен для решения смешанных задач и задач с разрывными коэффициентами. Для увеличения скорости сходимости последовательностей разложений для анизотропного уравнения диффузии значения параметров ю(г) следует выбирать при є >> 1 меньшими, а при є << 1 большими, чем значения параметров в (9) и (10). При є и 1 их можно оставить без изменения.
2) Для задач с разрывными коэффициентами как при точном, так и при приближенном выборе параметров последовательностей касательных и двухчастотных разложений (9) и (10), скорость сходимости практически не зависит от величины скачка.
В заключение приведем результаты численных исследований для уравнений с переменными коэффициентами и покажем эффективность метода последовательностей касательных и двухчастотных разложений. Касательные разложения строились при помощи тестовых векторов e(l) = (эт(лю(г) zh))i=1,...,n, ю(г) задавались по закономерностям (9) и (10). Число разложений k и шаг сетки связаны равенством h = 1/ 2k .
Пусть коэффициент ф^, y) уравнения V ^(x, y) Vu) = F имеет вид
ф^, y) = 1 + ql(x(1 - x) + y(1 - y)) при ql > 0.
Численные исследования показали, что эффективная скорость сходимости последовательностей как касательных, так и двухчастотных разложений при таком выборе параметров ю(г) от коэффициента ql практически не зависит и не превышает при всех значениях q1 є [0, 103 ] вплоть до h = 1/1024 значения 0,54, даже без использования метода сопряженных градиентов. Это не является характерным для задач с быстро меняющимися коэффициентами.
Пусть теперь ф^, y) = 1 - exp(- xy). Особенностью этой задачи является то, что ее дифференциальный оператор будет вырождающимся эллиптическим оператором. Решение задач подобного рода связано с большими трудностями. Значения эффективной скорости сходимости для последовательностей касательных разложений для этой задали при всех 1/1024 не превосходят 0,56, а для двухчастотных — 0,69. Использование метода сопряженных градиентов позволяет уменьшить эти значения соответственно до 0,48 и 0,57. Это поведение эффективной скорости сходимости является характерным и для других задач с переменными коэффициентами: если не применять метод сопряженных градиентов, то последовательность касательных разложений дает, как правило, лучшую скорость сходимости, чем двухчастотных, а применение метода сопряженных позволяет уменьшить эту разницу.
Рассмотрим ф^, y) = 1 + q2 sin(14^x) sin(14^y). Несмотря на быстрое колебание значений коэффициента ф^, y) при q2 < 0,5, скорость сходимости последовательностей касательных и двухчастотных разложений примерно такая же, как и для уравнения Пуассона. При больших значениях q2 скорость сходимости несколько замедляется, не превосходя при всех значениях q2 < 0,99 вплоть до h = 1/1024, значения 0,81.
Таким образом, многочисленные численные эксперименты показали, что предложенный способ выбора оптимальных параметров для последовательностей касательных и двухчастотных разложений позволяет достичь высокой скорости сходимости не только для модельных задач, но и для задач с переменными коэффициентами.
Работа выполнена при поддержке Российского фонда фундаментальных исследований, проект № QS-Q1-QQ431.
Список литературы
1. Buzdin К. Tangential decomposition // Computing (1998) 61:257-276.
2. Buzdin К., Wittum G. Two-Frequency Decomposition // Numerische Mathematik. 2004. Vol. 97. P. 269-295.
3. Wagner C. Tangential frequency filtering decompositions for symmetric matrices // Numer. Math. (1997) 78: 143-163.
4. Wittum G. Filternde Zerlegungen-Schnelle Loser fur grosse Gleichungssysteme / / Teubner Skripten zur Numerik, Band 1, Teubner-Verlag, Stuttgart, 1992.
5. Буздин А. А., Васильева Е. А., Латышев К. С. О последовательностях двухчастотных разложений // Вестник Калининградского государственного университета. 2006. Вып. 10. C. 64 — 69.
6. Васильева Е. А. Неполные блочные разложения, основанные на аппроксимантах Паде: дис. ... канд. физ.-мат. наук. Калининград, 2008.
7. Gander M. J.; Nataf F. A preconditioner based on the analytic factorisation of the elliptic operator. Numerical Linear Algebra with Applications. 2000. V. 7, № 7—8. P. 543-567.
8. Самарский А. А, Николаев Е. С. Методы решения сеточных уравнений. М., 1978.
9. Федоренко Р. П. Приближенное решение задач оптимального управления. М., 1978.
10. Hackbusch W. Iterative solution of large sparse systems of equations. New York, Springer-Verlag, 1993.
Об авторах
А. А. Буздин — канд. физ.-мат. наук, доц., РГУ им. И. Канта.
Е. А. Васильева — канд. физ.-мат. наук, доц., КГТу.
Authors
Dr. A. A. Buzdin — assistant professor, IKSUR. Dr. E. A. Vasilieva — assistant professor, KSTU.