Об одном варианте метода неполного блочного разложения Текст научной статьи по специальности «Математика»

70

доказывает однозначную разрешимость соответствующих разностных уравнений.

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

1. Ильин В.П. Прямой анализ устойчивости метода прогонки // Актуальные проблемы вычислительной математики и математического программирования. Новосибирск: Наука, Сибирское отделение, 1985. С. 189 — 201

2. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978. 519 с.

3. Кащенко Н.М., Захаров В.Е. Численный метод интегрирования системы уравнений переноса ионосферной плазмы // Доклады международного математического семинара. Калининград: Издательство КГУ, 2002. С. 287 — 290

Об авторе

Н.М. Кащенко — канд. физ.-мат. наук, доц., КГу.

УДК 519.615.5

АА. Буздин, Е.А. Васильева

ОБ ОДНОМ ВАРИАНТЕ МЕТОДА НЕПОЛНОГО БЛОЧНОГО РАЗЛОЖЕНИЯ

Предлагается метод построения предобуславливателя для решения больших систем линейных уравнений, основанный на неполном блочном разложении блочных трехдиагональных матриц. Метод является обобщением методов, изложенных в работах [1; 2]. Его практическая реализация для модельной задачи сходна с реализацией метода модифицированной матричной прогонки [9]. Этот метод выгодно отличается от последнего своей эффективностью в случае, когда число блоков велико.

A preconditioner for large systems of linear equations based on the incomplete block-decomposition for block-tridiagonal matrices is presented. This method generalizes methods developed by Buzdin and Wittum [1; 2]. Practical realization of this method for the model problem is similar to the realization of the method of complete matrix factorization [9], but much more effective if the number of the blocks is large.

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

Вестник КГУ. 2005. Вып. 1—2. Сер. Информатика и телекоммуникации. С. 70 — 76.

практическая экономичность и широкие возможности конструирования адаптивных алгоритмов для разных классов задач. Метод неполной факторизации был предложен в 50-е годы известным специалистом по вычислительной гидродинамике Н.И. Булеевым [3] и в последующие годы неоднократно переоткрывался. Он получил дальнейшее развитие в работах О. Аксельсона [4], В.П. Ильина [5], Г. Виттума [6] и других, например [7; 8].

В статье предлагается метод построения предобуславливателя для решения больших систем линейных уравнений, основанный на неполном блочном разложении блочных трехдиагональных матриц. Этот метод является обобщением методов, изложенных в работах [1; 2], где для модельной задачи матрицы полного блочного разложения заменяются их линейными аппроксимациями. Мы рассматриваем метод, основывающийся на построении рациональных аппроксимаций для этих матриц.

71

1. Основная идея метода

Пусть блочная матрица А подлежащей решению системы уравне-

ний АУ = Б, Ає Ип

1, У, Бє Ипгп имеет виц:

А =

Г О,

Ьх

0

0

0

Ь 2

О Б 0 Ь

п-1

п-1

п-1

(1)

где Ц, Ьуе И™ m, Lj — ^, Dj — , Dj, Lj > 0, 01 > Ь1, Оп > Ьп_1, Dj > Ь—1 + Lj,

] — 2, ..., п-1. Полное блочное разложение для матрицы А имеет вид:

А — (Ь+Т) Т-1 (и+Т), (2)

где L — blocktridiag{L^> 0, 0}, и — blocktridiag{0, 0, ^}, ; — 1,..., п, Т — Ь1оск-diag{Tl, ..., Тп}, причем

Т1 — Ц Tj — Ц - Ь- ТД Lj■_l,) > 1. (3)

Неполное блочное разложение получается, если в (2-3) заменить Т_11или 1— ТД ^ _ их некоторой аппроксимацией.

В работах [1; 2] предложен метод аппроксимации, основную идею которого проиллюстрируем, рассмотрев задачу, в которой в матрице А все блоки 0j — С, а ^ — -I. Эту задачу мы в дальнейшем будем называть модельной. Такого рода матрицы возникают, например, при дискретизации уравнения Пуассона с помощью пятиточечного шаблона. В этом случае блоки Tj — _^(С), где ^ - некоторые рациональные функции. Приближая их линейными функциями от матрицы С, мы придем к касательному или двухчастотному разложению, в зависимости от способа построения приближения. Например, двухчастотное разложение можно определить следующим образом [2].

72

А.А. Буздин, Е.А. Васильева

Определение 1. Пусть задана блочная трехдиагональная симметричная матрица А вида (1). Тогда ее двухчастотное разложение для точек (частот) X(1), X(2) имеет вид:

А — ^ + Т) Т -1^Т + Т), (4)

где Т —b1ockdiag{ Ту} и

~ ~ (1) //(я-(2)) _ //(А-(1)) (1)

Т1 — С, Т — /;(Х(1)) I^2^_(2Т-^(15-----(С-Я(1)1), у > 1. (5)

Параметры / (Я(г)), I — 1, 2 можно вычислить по формулам

Л(Я(г)) — Я(г), /у(Я(г)) — Я(г) -1 //^1)), у > 1.

Мы придем к определению касательного разложения матрицы А, если вместо секущей, проходящей через точки (Я(1), /у(Я(г))),

(Я(2), /у (Я(2))) кривой /у(Я), рассмотрим касательную в точке (X*, /(X*)) и сопоставим ей аналогично (5) линейную функцию матрицы С.

Легко проверить, что блоки Ту могут быть рассчитаны по следующим рекуррентным формулам:

Т1—с, т—с+у у т_1 - (ц(д+у) I, у > 1,

где ц(г) —V/уу > 1. Такое представление матриц Ту позволяет

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

Определение 2. Пусть А - симметричная блочная трехдиагональная матрица вида (1) и пусть заданы два тестовых вектора в1, е2 е Ит . Тогда для матрицы А можно определить двухчастотное разложение вида (4), где

а параметры

у — (Lj_ 1ег, ег)/(Туег, ег), г — 1,2.

Определение касательного разложения, отвечающее тестовому вектору е, получается формально из определения 2 при ц(2) ^ц(1) — (Lу_1е, е)/(Туе, е).

В работах [1; 2] доказано существование двухчастотного и касательного разложения для матриц вида (1). Результаты теоретического анализа и проведенных численных экспериментов показали, что особенно эффективным является использование этих разложений в качестве предобуславливателей в методе сопряженных градиентов. В этом случае доказано, что для модельной задачи скорость сходимости имеет порядок 1 - 0( И1 /3). Численные эксперименты показывают аналогичное поведение скорости сходимости и для более широкого класса задач.

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

Т-1 (С) — Г”аку(С -ук)1 )-1. (6)

В этом случае матрицы ТТу (С) можно эффективно обратить, выполняя

Ыр прогонок. Подобное разложение для блоков Ту(С) возникает в модифицированном методе матричной прогонки (см. напр. [9]), который требует для своей реализации 0(тп2) арифметических действий и поэтому является эффективным только тогда, когда число блоков в матрице А мало. Принципиальным отличием предлагаемого метода является то, что в сумме (6) число слагаемых Ы” мало и не зависит от номера блока.

Для построения приближения сопоставим функции от матрицы /у(С) функцию /(х), где хеИ и приблизим ее рациональной функцией

Ту (х) вида (6). Для того чтобы однозначно определить функции (х),

достаточно в некотором количестве точек Хк, к—1, ..., г задать значения функции и ее производных так, чтобы общее количество условий равнялось 2Ыр. В предлагаемом способе построения приближений значения функции / у (х) и ее производных в точках хк берутся равными значениям функции /(х) и ее производных в этих точках.

Значения функций /у (х) могут быть рассчитаны, например, по рекуррентным формулам

/1(х) — х, / — х-1//ул(х), у > 1.

Значения производных могут быть получены дифференцированием этих соотношений. Выбор точек, в которых вычисляются значения приближаемой функции и производных, обеспечивающих наилучшую скорость сходимости, представляет собой отдельную и до конца не решенную задачу. Отметим только, что если размер матрицы С невелик, то выбор в качестве точек хк всех точек спектра матрицы С приводит к точному решению задачи. В общем случае в качестве точек хк следует выбирать точки из отрезка \kmin, Xmax], где ^т^ Ятах _ соответственно минимальное и максимальное собственные значения матрицы С.

Для построения рациональных приближений может быть использован ряд методов, см., напр. [10; 11]. Разлагая полученную рациональную функцию на простейшие ~у(х) — Т1 а у /(х - уку) и сопоставляя ей функцию от матрицы С вида (6), получим искомую аппроксимацию.

2. Корректность метода и результаты численных экспериментов

Возникает вопрос о корректности предлагаемого метода, а именно, доказательстве того, что представление (6) будет существовать и что

74

А.А. Буздин, Е.А. Васильева

матрицы С - 2Хку1 из (6) будут невырожденными. Для модельной задачи ответ на этот вопрос может быть получен из следующих соображений. Как известно [9], в этом случае, используя полиномы Чебышева второго рода и;(х) степени 7, блоки Ту, 7=1,..., п можно записать в явном виде:

Ті = С, Ту (х) = и--1 (С/2) и (С/2), 7 > 2.

Так как нули полинома Ц(х) простые и расположены на промежутке (-1, 1), то для Ту(х) справедливо представление

Т-1 (х) = 21а м/(х - 2хм), (7)

где хк- корни полинома Чебышева Ц(х).

Из представления (7) следует, что функции 'Г-1 (х) являются

функциями1. Более того, так как корни знаменателя находятся на интервале (-2, 2), рассматриваемые функции будут принадлежать классу ЭТ[-2, 2]. Поэтому, по существу, вопрос о корректности метода сводится к вопросу о возможности приближения функции класса Ш[а, Ь] функцией того же класса. Ответ на него дает следующая теорема.

Теорема. Пусть дана функция /(х) є Ща, Ь], принимающая в точках хк значения /к, к=1,...,2Ыр. Тогда существует единственная функция g(x) того же класса, являющаяся правильной рациональной дробью вида g(x)=BNp_l(x)/DNp(x) и принимающей в точках хк, к=1,..., 2ЫР те же значения.

Приведем основные этапы доказательства. Следуя общему принципу решения интерполяционных задач такого типа, можно получить описание совокупности всех функций класса Ш[а, Ь], принимающих в точке х1 значение /1 = /хг) > 0:

((х\ = /1 (х1 - а) + ю(х)(х - а)(Ь - х1) ^

(х - а) 1 + ы(х)(Ь - х) '

где ю(х) є ЭТ[а, Ь]. Далее, для того чтобы функция/х) в точках хк, к=2,..., 2ЫР принимала значения /к, необходимо и достаточно, чтобы выполнялись 2Np-1 условия на функцию ю(х), которые можно получить, выразив из (8) функцию ю(х) через /(х). Таким образом, мы получаем исходную задачу с числом условий на единицу меньше. Многократное применение этого приема позволяет построить совокупность функций класса Ш[а,Ь], удовлетворяющих всем 2ЫР условиям. Совокупность таких функций дается формулой

А^(х)ю(х) + BNp-і(х)

(х) =--------------------,

CNV (х)ю(х) + DNv (х)

где ю(х) є ЭТ[а, Ь], а Ап(х), Вп(х), Сп(х), Dn(x)- алгебраические полиномы степени п, которые можно построить рекуррентно. При ю(х) = 0 получим искомую дробь.

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

1 Определение и свойства функции класса Ш[а, Ь] можно найти в [12].

С целью определения эффективности метода были проведены численные расчеты для модельной задачи с матрицей С = tridiag{-1, 4, -1}, отвечающей с точностью до нормировки задаче Дирихле для уравнения Пуассона. В качестве расчетной области выбирался квадрат с равными шагами к в обоих направлениях. В этом случае все собственные значения матрицы С находятся на интервале (2, 6), поэтому точки аппроксимации хк, к = 1, ..., 2N удобно задавать по формуле

хк = 2 + 48т2(лХк/2п).

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

Таблица 1

Скорость сходимости и оптимальные параметры для разложения на четырех точках аппроксимации

к 1/32 1/64 1/128 1/256 1/512 1/1024

Х.1, Х.2, 1.0; 2.5; 1.2; 2.3; 1.2; 3.5; 1.3; 3.9; 1.4; 4.5; 1.5; 5.0;

Яд, Я4 5.5; 17.5 8.0; 27.0 11.0; 49.0 17.5; 75.0 26.0; 120.0 38.0;185.0

П(4)пит 0.0320 0.0787 0.1504 0.2035 0.2955 0.3914

П(4)сх 0.0146 0.0403 0.0723 0.1141 0.1664 0.2255

В табл. 2 приведены результаты расчетов для той же задачи с большим числом точек аппроксимации. Применение метода сопряженных градиентов приводило к уменьшению Цпит приблизительно в два раза.

Таблица 2

Скорость сходимости для разложений на шести и восьми точках аппроксимации

к 1/32 1/64 1/128 1/256

п(6) -Ч4 'пит 2.18е-3 8.27е-3 2.38е-2 4.95е-2

п(8) А|ч 'пит 1.72е-4 1.03е-3 5.35е-3 1.09е-2

Приведенные расчеты показывают высокую эффективность метода. С целью сравнения различных вариантов метода в табл. 3 представлены эффективные скорости сходимости п/ для различного числа точек аппроксимации. Так как для реализации одной итерации с 2Np точками аппроксимации требуется примерно в Np раз больше арифметических действий, чем для реализации одной итерации двухчастотного разложения, то в табл. 3 скорости сходимости всех методов пересчитаны в расчете на одну итерацию двухчастотного метода.

А.А. Буздин, Е.А. Васильева

Таблица 3

Эффективная скорость сходимости для разложений для различного числа точек аппроксимации

h 1/32 1/64 1/128 1/256

n(2V 0.402 0.567 0.700 0.799

n(4W 0.179 0.281 0.388 0.451

n(6)e/f 0.130 0.202 0.288 0.367

n(8)e/f 0.115 0.179 0.270 0.323

76

Из табл. 3 видно, что применение вместо двухчастотного разложения разложений с большим числом точек аппроксимации резко увеличивает эффективность метода. Различия между разложениями для числа точек 2Np > 4 не столь значительны. К сожалению, ограниченность объема статьи не позволила изложить обобщение этого метода на другие виды разложений и матрицы более общего вида. Отметим, что проведенные численные расчеты показывают его эффективность и в этом случае. Скорости сходимости для широкого класса задач близки к скорости сходимости для модельной задачи.

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

1. Buzdin A. Tangential decomposition // Computing (1998) 61:257276.

2. Buzdin A., Wittum G. TwoFrequency Decomposition // Preprint, Universität Heidelberg, 2000.

3. Марчук Г.И. Методы вычислительной математики. Учебное пособие для вузов. М.: Наука, 1977.

4. Axelson O., Polman B. A robust preconditioner based on algebraic substructuring and two-level grids // W. Hackbusch, ed., Robust multigrid methods, NNFM Bd.23, ViewegVerlag, Braunschweig, 1989.

5. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Наука. Физматлит, 1995.

6. Wittum G. Filternde Zerlegungen-Schnelle Löser für grosse Gleichungssysteme // Teubner Skripten zur Numerik, Band 1, Teubner-Verlag, Stuttgart, 1992.

7. Wagner C. Tangential frequency filtering decompositions for symmetric matrices // Numer. Math. (1997) 78: 143163.

8. Kettler R. Analysis and comparison of relaxation schemes in robust multigrid and preconditioned conjugate gradient methods // W. Hackbusch und U. Trottenberg, ed., Multigrid methods, Proceedings, Köln-Porz, 1981. Springer-Verlag, Berlin, 1982.

9. Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

10. Stoer J., Bulirsch R. Introduction to numerical analysis. Springer-Verlag, New York,1992.

11. Baker G.A. jr. Essentials of Pad'e Approximants. Academic Press, New York, 1975.

12. Крейн М.Г., Нудельман А.А. Проблема моментов Маркова и экстремальные задачи. М.: Наука, 1973.

Об авторах

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

Е. А. Васильева — ассистент, КГу.