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

Он начинает работу с задачей на мелкой сетке, характеризуемой векторами Ьб) и хй). 2. Он улучшает текущее приближение, подавляя высокочастотную погрешность: х(') = Я(50),хй)). 3. Он вычисляет невязку тбй приближенного решения хб), 4. Он аппроксимирует невязку тбй для мелкой сетки с помощью соседней, более грубой сетки: В(тб)). 5. Задача на более грубой сетке решается рекурсивно при нулевом начальном приближении: МСУ(4 В(тб) ), О).

Коэффициент 4 обязан своим происхождением множителю л~ в правой части уравнения Пуассона; этот множитель изменяется ровно в 4 раза при переходе от мелкой сетки к грубой. 6. Алгоритм интерполирует решение, полученное на грубой сетке, к более мелкой сетке: с(0) = 1п(МСЪ'(4 В(тй)), О)). 7. Вычитает поправку, вычисленную с помощью грубой сетки, из приближенного решения на мелкой сетке: хб) = хбй — сг0. 8. Снова уточняет текущее приближение: хб) = Я(50), хб)), Дадим краткое обоснование алгоритма (детали будут восполнены позже). Сделаем (индуктивное) предположение, что 40) есть точное решение уравнения Тб) с(10 = тб) = Тбй хб) — Ьб).

Перегруппировывая его члены, находим Тй) . (хб) — бб)) = Ьб). Таким образом, хбй — сб«) есть искомое решение. Название «11-цикл» имеет следующее объяснение. Изобразим алгоритм схематически на плоскости переменных (номер сетки «, время), помещая на график точку при каждом рекурсивном обращении к МСЧ. Тогда получится нечто вроде рис. 6.11, где работа алгоритма начинается с обращения к МСУ(Ь®, х<е) ) (см. точку в левом верхнем углу). Затем происходят обращения к МСЧ для сетки 4, сетки 3, и так далее, пока не дойдем до самой грубой сетки 1, после чего сетки используются в обратном направлении с возвращением к сетке 5. Исходя только из того, что операторы Я, В и 1п алгоритма МС"11 заменяют текущие сеточные значения некоторыми взвешенными средними по со- 350 Глава б.

Итерационные методы для линейных систем Время Рис. 6.11. МСЧ. 0(4') = 0(4") = 0(п). а=1 Полный многосеточиый метод Неш окончательный алгоритм, называемый полным мпозосеточным методом (РМС), использует описанный выше алгоритм МСЧ как строительный блок. Алгоритм 6.17. РМОс 1ипсйоп РМС(МЯ~, х60) вычисляет приближенное решение х~я) ... хорошей точности для задачи Р® решить задачу Р<Ц точно, определив хйй 1огз= 2 сок х<0 = МСЧ(Ь<0,1п(хй Ц)) епс( 1ог В словесном описании алгоритм делает следующее: 1. Решает самую простую задачу Р00 точно. седним узлам, мы уже можем провести анализ сложности алгоритма в смысле символа 0( ).

Поскольку каждый из операторов производит фиксированную работу в каждом узле сетки, общее количество выполняемой им работы пропорционально числу узлов. Таким образом, если точка на букве Ч в нашем Ч-цикле соответствует сетке с номером г, то с этой точкой связана работа 0((2' — 1)~) = 0(4') операций. Если самой мелкой сетке соответствуют уровень я и и = 0(4Я) неизвестных, то полная стоимость цикла выразится геометрической прогрессией 351 6.9. Многосеточные методы Время Рис. 6.12. РМС-цикл. 2. При имеющемся решении хр Я1 более грубой задачи РВ '1 отображает хр '1 в начальное приближение х61 задачи РВ1 на более мелкой сетке: Хп(хр '1).

3. Решает задачу на более мелкой сетке с указанным начальным приближением, пользуясь алгоритмом МСЧ: МСЧ(Ь61, 1п(хр Я1)). Теперь можно провести анализ сложности в смысле символа О( ) для алгоритма РМС в целом. Рис. 6.12 дает изображение алгоритма на плоскости переменных (номер сетки 1, время). Каждому обращению к МСЧ во внутреннем цикле алгоритма РМС соответствует один символ «Ч» на рисунке. Как и прежде, стоимость такого символа, начинающегося с уровня в', составляет 0(4').

Поэтому общая стоимость алгоритма снова выражается геометрической прогрессией Этот результат оптимален, поскольку на каждое из п переменных приходится одна и та же, не зависящая от п работа. Он объясняет элемент табл. 6.1, соответствующий многосеточному методу. Ма11аЬ-реализацию многосеточного метода (как для одномерной, так и для двумерной задачи) можно найти в НОМЕРАСЕ/Ма$1аЬ/МС ВЕАПМЕ.ЬФш1. 6.9.2. Подробное описание многосеточного метода дпя одномерного уравнения Пуассона Теперь мы детально разберем конструкцию различных операторов Я, Л и 1п, составляющих многосеточный алгоритм, и дадим доказательство его сходимости. Мы сделаем это для одномерного уравнения Пуассона, поскольку здесь будут схвачены все особенности поведения алгоритма, но записи будут проще.

В частности, вместо последовательности вложенных двумерных задач теперь можно рассматривать показанные на рис. 6.13 аналогичные одномерные задачи. 362 Глава б. Итерационные методы для линейных систем 2 2 2 2 2 1 зйз>: одномерная сетка из 9 узлов 7 нензвестнаа узлы, помеченные цифрой 2, сосавляюг следуююую более грубую сетку Рх~у: одномерны сетка из 5 узлов 3 неизвестных узлы, помеченные цифрой 1, составляют следующую более грубую сетку уи гУ; одномерная сетка из 3 узнав одно неизвестное Рнс.

6.13. Последовательность сеток, используемых в одномерном многосеточном методе. Как и выше, будем обозначать через Р(6 задачу, решаемую на сетке з, а именно Тб> хбй = Ь<'7, где Тбй = Тгз„а Хг = 2' — 1, как и раньше. Начнем с описания оператора сглаживания Я, который представляет собой вариант взвешенного метода Якоби. е =Ве =В ео = (Я(1 — юЛ/2)ют) ео = Я(1 — юЛ/2) Е~еп, откуда Я~е = (1 — юЛ/2) ю ~еп, или (л~е )у = (1 — юЛ/2) (Хтеп)у. Поскольку е = ю(Ете ) есть сумма столбцов матрицы Я, где числа (юге ) выступают как веса, иначе говоря, е есть сумма синусоид различной частоты (см. рис. 6.2), то число (юге,„) называется у-й частотной компонентой погрешности е .

Собственное значение Лу(В ) = 1 — юЛ,/2 определяет, с какой скоростью стремится к нулю эта компонента. На рис. 6.14 показан график Лз(В„) для М = 99 и различных значений веса ю. При ю = — и 3 ) —, т. е. для правой половины частот Л, имеем з зч ~Л,(В )~ < -'.

Это означает, что, независимо от Лг, соответствующие компоненты погрешности (лте )у умножаются на каждой итерации на -' или меньшее число. Как показывает рис. 6.16, низкочастотные компоненты погрешности не , Оператор сглаживания в одномерном случае В этом разделе для упрощения записей верхние индексы в символах ТП>, хбй и Ьбй будут опущены. Пусть Т = ЯЛЕт — спектральное разложение матрицы Т, описанное в лемме 6.1. Стандартный метод Якоби для решения системы Тх = Ь задается формулой х ч.г = Вх + с, где В = 1 — Т/2 н с = Ь/2. Рассмотрим взвешенный метод Якоби х эз = В х +с, где В = 1 — юТ/2 н с = юЬ/2; значение ю = 1 соответствует обычному методу Якоби. Заметим, что спектральное разложение матрицы В имеет вид В = Я(1 — юЛ/2)хт. Как обычно, собственные значения этой матрицы определяют скорость сходимости взвешенного метода Якоби.

Пусть е = х — х есть погрешность на т-й итерации метода, тогда 353 6.9. Миогосеточные методы ав ав ол ав о.« и о ю ю во «о во ю во во ю кю Рис. 6.14. График спектра матрицы Я для 1«' = 99 и «о = 1 (метод Якоби), «о = 1/2 и «о = 2/3. слишком убывают. Итак, взвешенный метод Якоби для ю = в эффективен в уменьшении высокочастотной погрешности. Поэтому наш оператор сглаживания для уравнения (6.51) соответствует одному шагу взвешенного метода Якоби для «о = ф (6.54) 'т'(Ь,х) = Вв/в ° х+Ь/3. Будем писать В /в вместо Вт/в, если нужно указать сетку «, на которой () действует Вг/в.

На рис. 6.15 показан результат выполнения двух шагов типа Я для г = 6; эта задача имеет 2' — 1 = 63 неизвестных. Рисунок состоит нз трех групп графиков. Первая группа показывает начальное приближение х,„и погрешность е,„ после соответствующего примененця оператора о'. Точное решение есть синусоида, образованная точками на левом графике каждой группы. На том же графике приближенное решение показано сплошной линией. На средних графиках изображена погрешность, причем под графиками приведено значение ее 2-нормы.

На правых графиках показаны частотные компоненты погрешности, составляющие вектор Я~е . Можно видеть на этих графиках, как применение оператора о' подавляет правую половину частотных компонент. Это же явление обнаруживается в средних и левых графиках тем, что приближенное решение становится глаже. Причина в том, что высокочастотная погрешность проявляет себя как «биения», а низкочастотная погрешность выглядит как «гладкая» функция. Поначалу норма погрешности убывает быстро — с 1.65 до 1.055. Однако затем убывание замедляется, поскольку высокочастотные компоненты погрешности уже в значительной мере подавлены.