Дж. Деммель - Вычислительная линейная алгебра, страница 81

Наша задача — показать, что собственные значения отделены от 1 расстоянием, не зависящим от г. Номера, упоминаемые в приводимой ниже записи, — это номера строк в алгоритме 6.16. (а) хбй = Я(Ь(г), х(1)) = Вт ~ хбй + Ь®/3 согласно строке 1) и уравнению (6.54), (Ь) т(0 = Тбй .х(0 Ь<0 согласно строке 2), Н(0 = 1п(МСЪ'(4. В(т®),0)) согласно строке 3) =1п((ТИ Ц) (4 В(тбй))) по нашему предположению о том, что задача на грубой сетке решается точно =1п((Т' ~~~ (4.Р,' 'тбй)) согласно уравнению (6.55) 6.9. Многосеточные методы 361 (с) = Р' ((Т)' ~)1 (4 ° Р' ~с<О)) согласно уравнению (6.56) (Н) хб) = хб) — дб) согласно строке 4) (е) х<О = $(Ь(з), х(1)) = В® х(О + ЬО'/3 согласно строке 5). Чтобы получить уравнение для пересчета вектора ошибки еб), вычтем из строк (а) и (е) тождество х = В~'~~ х + Ь)0/3, из строки (Ь) — тождество 0 = Т)О х — Ь)О и, наконец, из строки (а) тождество х = х.

В результате имеем: (а) е(О = В( ~ е(О г/зе (Ь) г'О = Т~О . е~О, (с) )160 = Р),([ТО Ц~ '(4.Р,.' 'т(О)), (а) е(О = е01 — а(О (е) еб) = В~ ) е(О г/з Подставляя каждое из этих уравнений в следующее, приходим к формуле, показывающей, как изменяется вектор ошибки в Ч-цикле: гй ~~ ) ВО) 1 — Р (Т) Ц1 . (4. Р ~Т(О) 1~~ ) ° ~® = М е10.

(6.57) Теперь нужно вычислить собственные значения матрицы М. Упростим внача- ле уравнение (6.57), используя соотношения Р г = 2. (Р,' ')т и ТО Ц = 4- Р,' ~Т)О)Рз ) —— 8 ° Р,' ~Т)О(Рз г)т (6.58) (см. вопрос 6.15) . Подставляя их в выражение для М в формуле (6.57), находим е=еГ (1 — (Р! ') - (Р.' т~~)Р'. ) $ .)Р.' т~~) )я)), или, опуская индексы, чтобы упростить запись, Вг/3 (~ Р (РТР ~ РТ /Вг/3' (6.59) Используем теперь то обстоятельство, что все матрицы, составляющие М (т. е.

Т, Вг/3 и Р), могут быть (почти) диагонализованы с помощью матриц Я = Я)О и ЯО Ц, состоящих из собственных векторов соответственно задач Т = ТО) и ТО '). Вспомним, что Я = Ят = Я ', Т = ЗЛЯ и Вг/3 = Я(Х вЂ” Л/3)Я =— ЕЛлЯ. Предоставляем читателю проверку соотношения Ф ПРУ(О = Лр, где Лр — почти диагональная матрица (см. вопрос 6.15), т.

е. (+1+ сов гз)/с/8, если /с = у, ~вне = ( — 1+ сов гг)/~/8, если й = 2* — У, (6.60) О, в противном случае. 362 Глава 6. Итерационные методы для линейных систем Это позволяет написать гмг = [гл,дзг) х 1 — (гР'гО-'>) <[г<'-оы)[гтрк)[Л"'го-0)~ ~ ~ х ~ » ~ | к[го-ИРг)[гтрк)) [гя„,г) =Л (1 — Л'[Л ЛЛ']-'Л Л) Л„. Поскольку Я = Я ', матрица ЯМУ подобна М, следовательно, имеет те же собственные значения. Кроме того, матрица ЯМУ почти диагональна, что означает: ненулевые элементы в ней могут находиться только на главной диагонали и «пердиагонапн» (т.

е. диагонали, идущей из левого нижнего угла матрицы в ее правый верхний угол). Это позволяет вычислить собственные значения матрицы М в явном виде. Теорема 6.11. Независимо от»', матрица М имеет собственные значения 1/9 и О. Поэтому многосеточный метод сходится с фиксированной скоростью, не зависящей от числа неизвестных. По поводу доказательства теоремы см. вопрос 6.15. Более общий анализ сходимости можно найти в [268].

Реализация данного алгоритма обсуждается в вопросе 6.16. Интернет- страница [91) содержит ссылки на обширную литературу по этой теме, программы и прочее. 6.10. Декомпозиция области Декомпозиция области при решении разреженных систем линейных уравнений является предметом текущих исследований. Обзоры последних результатов в этом направлении можно найти в [49, 116, 205) и, особенно, в [232]. Здесь мы приведем лишь простые примеры. Потребность в методах, отличных от описанных выше, проистекает из нерегулярности и размера реальных задач, а также из необходимости приспособить алгоритмы к параллельным компьютерам.

Самые быстрые методы из тех, что обсуждались до сих пор, а именно методы, основанные на блочной циклической редукции, алгоритме РРТ и многосеточном подходе, работают только или быстрее всего для очень регулярных задач типа нашей модельной задачи, т. е. уравнения Пуассона в прямоугольнике, дискретизованного посредством равномерной сетки. Однако с решением реальной задачи может быть связан не прямоугольник, а менее регулярная область, представляющая физический объект вроде крыла на рис. 2.12. Этот рисунок показывает также, что в тех частях области, где, как ожидается, решение является менее гладким, сетка может сгущаться по сравнению с частями, где решение гладкое. Кроме того, нам могут встретиться более сложные уравнения, чем уравнение Пуассона, или даже различные уравнения в разных подобластях.

Регулярна задача или нет, она может из-за своего размера не поместиться в памяти компьютера и тогда ее, возможно, придется решать «по частям». Или же мы можем захотеть разбить задачу на части, которые допускают одновременное решение на параллельной машине. 363 6.10. Декомпозиция области Теория декомпозиции области рассматривает все эти вопросы с тем, чтобы указать систематический способ конструирования «гибридных» алгоритмов из более простых методов, обсуждавшихся в предыдущих разделах. Эти более простые методы применяются к меньшим и более регулярным подзадачам исходной задачи, после чего из полученных частичных решений «монтируется» решение задачи в целом. Если вся задача не помещается в памяти компьютера, то подзадачи могут решаться по очереди, а на параллельном компьютере — одновременно.

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

В некоторых случаях (например, для задач, достаточно похожих на уравнение Пуассона) эта теория приводит к «оптимальным методам» (где выполняется работа 0(1) на одно неизвестное). Наше обсуждение будет разделено на две части, соответствующие методам без перекрытия и методам с ерекрытием. 6.10.1. Методы без перекрытии В литературе методы этого типа называют еще методами нодструктур, или методами дополнений Шура. Такие методы используются уже несколько десятилетий, особенно специалистами по строительной механике, для разбиения больших задач на части, могущие быть размещенными в памяти компьютера. Для простоты, проиллюстрируем данный тип методов с помощью обычного уравнения Пуассона с граничными условиями Дирихле, дискретизованного посредством 5-точечного шаблона.

Однако уравнение будем рассматривать не на квадрате, а на й-образной области. Эту область можно разбить на две: малый квадрат и большой с вдвое большей стороной, причем малый квадрат примыкает к нижней части правой стороны большого квадрата. Построим метод, использующий нашу способность быстро решать задачи на квадратах. Для грубой сетки, показанной на рисунке, каждый внутренний узел Ь- образной области помечен номером (стояшим слева сверху от узла). 364 Глава 6.

Итерационные методы для лннейных систем Заметим, что вначале пронумерованы внутренние узлы каждой из подобластей (номера с 1 по 4 и с 5 по 29) и только потом узлы на их общей границе (номера 30 и 31). В результате получается такая матрица: -1 -1 -1 4-1 -1 4 Аы 0 Азз Азз -4гз Азз Здесь Аы = Тг хз А22 — Тз ха и Азз — Тз х1 = Тз+ 21з, где матрица Тгз определена в (6.3), а матрица Тнхм — в (6.14). Одним из важнейших свойств этой матрицы является то, что А1 з = О, поскольку между внутренними узлами двух подобластей нет прямого сцепления.

Сцепление осуществляется лишь через занумерованные последними узлы общей границы подобластей (узлы 30 и 31). Таким образом, блок Аы соответствует сцеплению между границей и малым квадратом, а блок Азз — сцеплению между границей и большим квадратом. Чтобы понять, как извлечь выгоду из специальной структуры матрицы А при решении системы Ах = Ь, запишем блочное ЫП1-разложение этой матрицы 365 б.10.

Декомпозиция области в виде А= [ 0 1 0 0 Аы 0 Агг 0 1 0 0 1 0 ° 0 Агг Агз АтА 1 АтА 1 1 ~0 0 Я1 1 0 0 1 где матрица Я = Агз — А~~А~~'Агг — АггАгг Агг (6.61) 1 0 0 0 1 0 0 0 Я ~ А '= [ А~~~ 0 — АгггАгз 0 Агг' — Агг'Агз 0 0 о о~ 0 1 0 Ат А — ~ Ат Агг' 11 Таким образом, чтобы умножить вектор на А г, нам нужно умножать на блочные элементы этой факторизованной формы, а именно на блоки Агг и Агз (и транспонированные к ним), А,~г, Аг ' и Я '. Умножение на Аш и Агг стоит дешево, так как эти блоки очень разреженьь Умножение на А~,~ и Аг ' тоже дешево: ведь подобласти были выбраны так, чтобы были применймы алгоритм РГТ, блочная циклическая редукция, многосеточный метод и некоторые другие быстрые методы из обсуждавшихся выше.

Остается объяснить, как происходит умножение на матрицу Я ~. Поскольку на границе гораздо меньше узлов, чем в подобластях, порядок блоков Агг и Я много меньше порядка блоков Аы и Агг, этот эффект еще более выражен при сгущении сетки. Как и А, блок Я симметричен и положительно определен, но является (в данном случае) плотным. Чтобы вычислить его в явном виде, потребовалось бы решить по одной задаче для каждой подобласти на каждый из общих граничных узлов (что соответствует произведениям Агг'Агг и Агг~Агг в формуле (6.61)). Это, разумеется, можно сделать, после чего можно было бы факторизовать Я посредством плотного алгоритма Холесского, а затем решить систему с этой матрицей. Но такой способ действий был бы чересчур дорогостоящим, куда дороже, чем простое умножение вектора на матрицу Я, для чего, согласно (6.61), понадобится лишь одно решение каждой задачи на подобласти.