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

Он численно неустойчив, так как матрицы А(") быстро растут: ((А(«)(~ 'бА(' ')((г 4г, поэтому при вычислении вектора Ь('+ значения векторов Ьгз, и Ь 1+, не будут из-за округлений отражены в результате. («) [«) 2. Если А — трехдиагональная матрица, то А(') имеет ленту ширины 2" + 1, следовательно, быстро заполняется. Поэтому умножение на А(") и решение систем с этой матрицей становятся дорогостоящими операциями.

Вот как можно исправить второй недостаток. Заметим, что А(') представляет собой многочлен р,(А) степени 2' от матрицы А, где ро(А) = А и р«+г(А) = (р«(А)) — 21. Лемма 6.10. Пусть 1 = 2 сов В. Тогда р (() = р,(2 сов В) = 2 соз(2"В). Доказательство. Утверждение леммы вытекает из простого тригонометрического тождества.

П Заметим, что р, (1) = 2 сог(2' агссоз( г )) = 2Тт- ( г ), где Тт. ( . ) есть многочлен Чебышева (см. разд. 6.5.6). Лемма 6.11. Имеет место представление р«(() = П. (1 — 11), где 1 2 сов(7Г 2 ) ° Доказагпельство. Нули многочленов Чебышева указаны в лемме 6.7.

П Таким образом, А(') = Пг „(А — 2 сов(.г1зг=„')1), поэтому решение системы А(')г = с эквивалентно решению 2" систем с трехдиагональными матрицами коэффициентов А — 2 сов(хогг:„')1. Решение каждой из них гауссовым исключением или алгоритмом Холесского обходится в 0()(г) операций. Чтобы добиться численной устойчивости, в алгоритм должны быть внесены дополнительные изменения. Окончательная версия алгоритма, принадлежащая Бунеману (Вппешап), описана в (47, 46]. Проанализируем сложность простого алгоритма; такова же будет сложность устойчивой версии. Умножение на трехдиагональную матрицу или решение трехдиагональной системы порядка Ю стоит 0(Ю) флопов. Поскольку 345 6.9.

Многосеточные методы А<") есть произведение 2" трехдиагональных матриц, умножение на А)') или решение системы с этой матрицей стоит 0(2"Ю) флопов. Поэтому стоимость внутреннего цикла на шаге 1) алгоритма составляет г,»т 0(2" Х) = 0(Л~) флопов, так как на нем пересчитывается Ю «.1 - г хт векторов Ь . Матриц (»-»1) ца А~" »") не вычисляется в явном виде.

Поскольку цикл на шаге 1) выполняется 1« — 1оягМ раз, полная стоимость шага 1) равна 0(Х 1оягЮ). Аналогично подсчитывается, что шаг 2) стоит 0(2"Ю) = 0(Юг) флопов, а шаг 3) — 0(Хг 1ояг Х) флопов, что дает для алгоритма в целом стоимость 0(А'г 1ояг Х) флопов.

Это объясняет элемент табл. 6.1, относящийся к блочной циклической редукции. Алгоритм обобщается на блочно-трехдиагональные матрицы, в которых все диагональные блоки равны симметричной матрице А, а все блоки на диагоналях, соседствующих с главной, равны симметричной матрице Р, перестановочной с А (т. е. ГА = АГ). См. по этому поводу вопрос 6.10. Такого рода структура часто встречается при решении линейных систем, возникающих при дискретизации дифференциальных уравнений, таких, как уравнение Пуассона. 6.9. Многосеточные методы Многосеточные методы были предложены в контексте решения дифференциальных уравнений с частными производными типа уравнения Пуассона, однако они пригодны для более широкого класса задач.

В отличие от обсуждавшихся до сих пор итерационных схем, скорость сходимости многосеточного метода, будучи независимой от размера )ч' сетки, не замедляется для болыпих систем. Как следствие, задачи с и неизвестными оказывается возможным решать за время 0(п), т. е. затрачивая фиксированную работу на каждое неизвестное. С точностью до (умеренной) константы, спрятанной в символе О( ), этот результат является оптимальным. Поясним, почему обсуждавшиеся ранее итерационные алгоритмы ие могут быть оптимальны для модельной задачи. Это утверждение, в действительности, верно в отношении любого итерационного алгоритма, вычисляющего приближение х +1 путем усреднения значений х и правой части 5 в соседних узлах сетки.

Под это описание подходят методы Якоби, Гаусса — Зейделя, КОЩЕЕ), ЯБОВ с чебышевским ускорением (для трех последних методов имеется в виду шахматное упорядочение) и всякий метод крыловского подпространства, основанный на матрично-векторных умножениях с участием матрицы Т)ч„гг (последнее верно потому, что умножение текущего приближения на матрицу Т)« „м равносильно усреднению по соседним узлам сетки). Предположим, что на сетке 31 х 31 задана правая часть 5 с единственным ненулевым элементом; она показана на левой верхней картинке рис.

6.9. Точное решение х показано на правой верхней картинке того же рисунка; отметим, что оно нигде не обращается в нуль и убывает по мере удаления от центра. Левая нижняя картинка рис. 6.9 показывает приближенное решение хлв, вычисленное посредством 5 шагов метода Якоби, исходя из нулевого начального вектора. Обратим внимание на то, что хдв равно нулю на расстоянии от центра, превышающем 5 узлов сетки. Причина в том, что усреднение по соседним узлам способно «передавать информацию» лишь на один узел за итерацию, а у начального 346 Глава б. Итерационные методы дли линейных систем Правая часть Точное решение 0.5 0.5 О 0 0 0 Наилучшее решение аа 5 шагов 1 5 шагов метода Якоби 1 0.5 0.5 0 0 0 0 Рис. 6.9. Ограничения, накладываемые усреднением по ближайшим узлам сетки. приближения была лишь одна ненулевая компонента в центре сетки.

Более общо, после 1« итераций отличны от нуля могут быть лишь значения в узлах сетки, удаленных от центра не более чем на й. На правой нижней картинке показано наилучшее решение хв,м а, которое можно получить за 5 шагов любым методом «ближайшего соседа». Это решение совпадает с х на расстоянии от центра, не превышающем 5 узлов сетки, и равно нулю на ббльшем расстоянии. Из рисунка видно, что погрешность и — ин„ца равна значению и на расстоянии шести узлов сетки от центра.

Это все еще большая погрешность. гРормализуя это рассуждение, можно показать, что какой бы алгоритм «ближайшего соседа» ни использовался, на сетке и х и потребуется по меньшей мере 0(1ойп) шагов, чтобы уменьшить погрешность в фиксированное число рнз, большее единицы. Чтобы получить лучший результат, чем 0(106 и) шагов (и, соответственно, сложность 0(п1окп)), нужно научиться «распространять информацию» более чем на один узел сетки за итерацию.

В многосеточных методах эта цель достигается путем обмена информацией с ближайшими узлами более грубых сеток; такие узлы могут оказаться гораздо дальше от данного, чем ближайшие узлы мелкой сетки. В многосеточных методах грубые сетки используются для реализации стратегии «разделяй-и-алас«науй» в двух родственных отношениях. Во-первых, начальное решение для сетки Ф х Ю строится посредством аппроксимирующей сетки (Ю/2) х (Ю/2), получаемой удержанием каждого второго узла исходной сетки Ю х ггг.

Более грубая сетка (Х/2) х (Х/2), в свою очередь, аппроксимируется сеткой (гт/4) х (Х/4), и так далее. Второй способ использования стратегии 347 б.9. Многосеточные методы «разделяй-и-властвуй» относится к частотпной области. Здесь погрешность нужно представлять себе как сумму собственных векторов, или синусоид с различными частотами. Тогда интуитивное описание будет таким: работа, выполняемая на конкретной сетке, уменьшает половину частотньгх компонент погрешности, для которых уменьшение не было достигнуто на более грубых сетках. В частности, работа, производимая для данной сетки (т.

е. усреднение приближенного решения в каждом узле сетки с использованием соседних узлов — некоторый вариант метода Якоби), делает приближенное решение глаже, что эквивалентно подавлению высокочастотных компонент погрешности. Ниже это описание будет развернуто в более подробное. ! Р!23: сетка узлов 5 х 5 сетка неизвестных 3 к 3 узлы, помеченные цифрой 3, ссставлвют слелуюшую более грубую сетку Р03: сеткаузаев3 х 3 сетка неизвестных 1 х ! Р<33: сетка узлов 9 х 9 сетка неизвестных 7 х 7 узлм, помеченные цифрой 2, ссставпюпг слелуюшую белее грубую сетку Рнс. 6.10.

Последовательность сеток, используемых в двумерном многосеточном методе. 6.9.1. Очерк многосеточного метода на примере двумерного уравнения Пуассона Мы дадим вначале общую формулировку нашего алгоритма, а затем объясним детали. Как и в случае блочной циклической редукции (рнзд. 6.8), более удобно рассматривать сетку неизвестных размера (2" — 1) х (2" — 1), а не сетку 2" х 2в, предпочтительную для г ГТ (разд. 6.7). Для лучшего понимания и более удобной реализации разумно добавить еще граничные узлы, где решение принимает известное значение О, тогда получится сетка размера (2" + 1) х (2" + 1), показанная на рис. 6.10 и 6.13. Положим еще гтгв = 2" — 1.

Обозначим через Рбй задачу решения дискретного уравнения Пуассона на сетке размера (2'+ 1) х (2!+ 1) с (2г — 1)з неизвестными, иначе говоря, на сетке размера (гггг + 2) х (уггг + 2) с гч'Р неизвестными. Задача Рбй характеризуется правой частью Ьбй и (неявно) размером сетки 2' — 1 и матрицей коэффициентов Т<г! = Тлгг „лгг. Приближенное решение задачи РО> обозначим через т(г!. Таким образом, ЬО! и юб! суть массивы размера (2г — 1) х (2« — 1), составленные из значений в узлах сетки. (Нулевые граничные значения присутствуют неявно.) Будет построена последовательность связанных задач Р!г), РО '3, РО а>, ..., Р(г) на 348 Глава 6. Итерационные методы для линейных систем все более грубых сетках; при этом решение задачи Р(' ') будет хорошим приближением погрешности решения задачи Р(г). Для того чтобы объяснить, как работает многосеточный метод, нам потребуется ввести операторы, которые либо улучшают приближенное решение для данной сетки, либо преобразуют данную задачу в родственную задачу для другой сетки: ° Оператор сглазсивания Я, будучи применен к задаче Р(') с приближенным решением х(1), вычисляет уточненное приближение х('): улучшенное хбй = Я(Ь(*), хбй).

(6.51) Улучшение состоит в подавлении «высокочастотных компонент» погрешности. Ниже будет объяснено, что это значит. Уточнение решения достигается усреднением значения в каждой точке сетки по его ближайшим соседям; оно представляет собой некоторый вариант метода Якоби. ° Оператор сужения В отображает правую часть Ь(') задачи Рбб в ее приближение Ь(' ') на более грубой сетке: Ь(1-1) — В(Ь(г) ) (6.62) Реализация этого оператора также сводится к взвешенному усреднению по ближайшим соседним узлам сетки.

° Оператор интерполяции 1п преобразует приближенное решение х(' «) задачи Р(' ') в приближенное решение х(') задачи Р('), поставленной на соседней, более мелкой сетке: (') — 1п( (' — ')) И здесь реализация оператора требует лишь взвешенного усреднения по ближайшим соседним узлам сетки. Итак, все три оператора выполняются путем замены текущих значений в кззкдом узле сетки некоторыми взвешенными средними по ближайшим соседним узлам.

Поэтому каждый оператор производит работу О(1) на одно неизвестное или работу О(п) для всех п неизвестных. Это объясняет низкую стоимость алгоритма в целом. Многосеточный Ч-цикл Сказанного выше достаточно для того, чтобы сформулировать основной алгоритм, так называемый многосеточный Ч-цикл (МСЧ). Алгоритм 6.16. МСЧ (строки пронумерованы для последующих ссылок): зипсйоп МСЧ(Ь('~, х(')) ... заменить приближенное . решение х(') задачи Рбй утаочненнт«м приближением »1» = 1 случай одного неизвестного вычислить точное решение х(1) задачи Р(1) возвратить х(~) е)зе 1) х(') = Я(Ь('),х(')) 349 6.9. Многосеточиые методы приближение вычислить невязку использовать рекурсию по более грубым сеткам уточнить решение на мелкой сетке 2) гб) = Тб) ° хб) — Ьб) 3) сК0) = 1п(МСЪ'(4 В(т(0),0)) 4) .0) .0) 10) 5) хб) = Я(50),хб)) снова уточнитпь приближенное решение возвратить хб) епй1 Если прибегнуть к словесному описанию, то алгоритм делает следующее: 1.