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

В результате, имеем рь, Ать-г т Рь т р, Арь г (6.42) Недостаток этой формулы для рь состоит в том, что требуется еще одно скалярное произведение рь Ать 1 помимо двух, необходимых для вычислет ния иы Поэтому мы выведем другую формулу, не использующую новых скалярных произведений. Чтобы получить формулу для иь, мы вначале умножим обе части уравнения (6.39) слева на рь А, а затем используем то обстоятельство, что векторы рь н рь г являются А-сопряженными (лемма 6.8).

Это дает 325 6.6. Методы крыловского подпрострапства Выведем вначале иную формулу для ой. Умножим обе части уравнения (6.37) слева на т~т, снова используем соотношение тг,тй = О и разрешим полученное равенство относительно пй,что дает т „ ой= т (6.43) ттАрй ' Окончательную формулу для рй найдем, приравнивая два выражения (6.41) и (6.43) для пй й (заметим, что индекс здесь уменьшен на единипу), переупорядочивая и сравнивая с соотношением (6.42): рй,Атй т т рй Арй т тй т т тй гтй-2 (6.44) Объединяя рекурсии (6.37), (6.38) и (6.39), а также формулы (6.41) и (6.44), получим окончательную реализацию алгоритма сопряженных градиентов.

Алгоритм 6.11. Алгоритм сопряженных градиентов: й=О;хо=О;то=Ь;р1 =Ь; тереа1 й = й + 1 в=А рй пй = (т„тй 1)7(рй г) хй = хй 1 + пйрй тй — — тй 1 — пйг рй+, = (тй тй)7(тй,тй-1) Рй«ч = тй + Рй+1рй пока число Отйбг не станет достаточно малым 6.6.4. Анализ сходимости метода сопряженных градиентов Мы начнем с анализа сходимости алгоритма СО, основанного лишь на значении числа обусловленности матрицы А. Этот анализ покажет, что число итераций метода, необходимых для снижения ошибки в фиксированное число раз, большее единицы, пропорционально квадратному корню из числа обусловленности. Такой анализ наихудшего случая дает хорошую оценку скорости сходи- мости метода для нашей модельной задачи, т. е.

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

Более подробно о деталях реализации, включая то, как определить, что «ОтйОг достаточно мало», можно узнать по адресу ХКТЫВ/сешр1аСег/Тетр1асез.'пСш1. 326 Глава б. Итерационные методы для лииейиых систем После вывода оценки, основанной на числе обусловленности, мы укажем, ког- да можно рассчитывать на более быструю сходимость. В качестве начального приближенного решения возьмем хо = О. Вспомним, что ха минимизирует А '-норму невязки г» = Ь вЂ” Аха по всем возможным выборам ха Е Ка(А, Ь). Это означает, что ха минимизирует величину [[Ь А«[[гА 1(«) (Ь А«)тА-1(Ь А«) (х «)тА(х «) по всем векторам «Е Ка = арап[Ь,АЬ,А«Ь,...,А" 1Ь).

Всякий вектор «Е Кь(А,Ь) может быть записан в виде « = 2 ~ о с«1А)Ь = ра 1(А)Ь = р(, 1(А)Ах, где рь-1(с) = 2 . о с«,бт — многочлен степени й — 1. Поэтому 1(«) = [(1 — р» 1(А)А)х) А[(1 — рь 1(А)А)х[ = (д»(А)х) тА(дь(А)х) = хада(А)Ада(А)х, где дь(с) — = 1 — ра 1(с) с — многочлен степени Ь, такой, что да(О) = 1. Заметим, что (да(А))г = дь(А), так как А = Ат. Обозначим через Яа множество всех многочленов степени й, принимающих значение 1 в нуле. Тогда Дха) = ппп 1(«) = ппп х да(А)Ада(А)х. (6.45) «Е(сь Ф»Е(«» Чтобы упростить это выражение, введем спектральное разложение А = („"«ЛЯт и положим (ь(~х = у. Получим 1(хь) = пбп 1(«) = ппп хт(д«ЯЛЯт))ЯЛЯт)(даЩЛЯ ))х ппп х~фда(Л)Я~)ЯЛЯ~)Яда(Л)()~)х ппп ух да (Л) Лда (Л) у вьем» ппп у 61ая(да(Л1)Л«да(Л»)) ° у Е Е(«ь ппп ~~ угЛ;(да(Л;)) е ей».

1=1 (а(» а') «. «1» Е»ЕЙ» Ь,Л«ЕЛ(А) ппп ( шах (да(Л;))г 1(хо), Е»ЕЯ» »,Л,ЕЛ(А] так как из хо = 0 следует, что 1(хо) = х Ах =у Лу = 2',,".ыьугЛ;. Поэтому [[ ~[[' - У(*~) ( гпш игах (да(Л1)), или [[та[[А-« < ппп шах [да(Л»)[. [[го[[А-«чьейл.ел(А) 327 6.6. Методы крыловского подлрострвнствв Итак, вопрос о скорости сходимости алгоритма СО сведен к вопросу о многочленах: насколько мал может быть многочлен и?ь1Я) степени й, удовлетворяющий условию и?в[0) = 1, когда С пробегает множество собственных значений матрицы А? Поскольку А положительно определена, ее собственные значения принадлежат отрезку [Л ы, Л ], где 0 < Л;и < Л,„.

Чтобы получить простую верхнюю оценку, будем вместо де[~) искать многочлен де[С), как можно более малый на всем отрезке [Л м, Л ах] и принимаюший значение 1 в нуле. Многочлен и?э[С) с этими свойствами легко конструируется из чебышевских многочленов Тл[г), обсуждавшихся в рэзд. 6.5.6. Вспомним, что ]Те[~)[ < 1 при [С[ < 1 и Тл[г) быстро возрастает по модулю при ]Я ) 1 [см. рис. 6.6). Положим теперь Липах + Лпиы 21'~ ~Т (Липах+ Лппп [ Лпиах Лпиип / ии Липах Лппп / Легко видеть, что $,[0) = 1 и если С 6 [Л;и, Л ], то Лпиах + Лип|а 2~ < 1.

Липах Лппп Поэтому < ппп шах [и?ь[Л )[ []ел [[А [[ГО[]А и — даэаЛиЕЛ(А) < 1 1 1 Т [л +л ) Т„[Я) Те[1+ „~,) [6.46) где к = Липах/Ливио есть число обусловленности матрицы А. Если число обусловленности к близкб к 1, то величина 1+ 2([к — 1) велика, а число 1 (Те[1+ „э ) мало; поэтому сходимость будет быстрой. Если к велико, то сходимость замедляется, причем А '-норма невязки гь стремится к нулю со скоростью 1 2 < Пример 6.14. Для модельной задачи на сетке Ю х Ф имеем к = 0[иигз), поэтому после к шагов метода СС модуль невязки приобретает множитель, приблизительно равный [1 — 0[И "))", т.е.

такой же, как в методе ЗОВ с оптимальным значением параметра верхней релаксации пи. Другими словами, для сходимости метода требуются 0[иЛи) = 0[я"гэ) итераций. Поскольку стоимость каждой итерации равна 0[п), общая стоимость алгоритма составляет 0[пег~). Это объясняет элемент таблицы 6.1, отвечающий методу 00. О Проведенный анализ, опиравшийся на число обусловленности, объясняет не все важные особенности в характере сходимости метода СО.

Как показывает следующий пример, имеет значение все распределение собственных чисел, а не только отношение наиболыпего из них к наименьшему. Пример 6.15. Рассмотрим график относительной невязки [[ге [[эДго[[э на шаге?с алгоритма СС для восьми различных линейных систем; этот график показан на рис 6.7. Относительная невязка [[ге[[э/[[го[[э измеряет скорость сходимости. Наша реализация метода прекращает работу, когда это отношение 328 Глава 6. Итерационные методы для линейных систем Относительные невязки как функции числа СО-итераций 10 1О 1О 10 1О 1О 10 м4 0 20 40 60 60 100 120 140 160 160 200 Рнс.

6.7. График относительных невязок, вычисляемых в алгоритме СС. становится меньше, чем 10 'а, либо когда уже проделаны й = 200 шагов. Выход происходит по первому выполненному из этих условий. Все восемь линейных систем имеют одну и ту же размерность и = 10 и одно и то же число обусловленности и 4134; тем не менее характер сходимости метода для них радикально различен.

Самая верхняя линия (состоящая из точек и тире) соответствует величине ЦГа(1+ а, ); согласно неравенству (6.46), эта величина является верхней границей для йга'йл- Дго~йл- . Оказывается, ЧтО ГРафИКИ ОтНОШЕНИй йтайа/ОГОйа И йтайл — з/йтейх-а ПРИМЕРНО ОДИНаКОВЫ, поэтому на рисунке представлены графики первого отношения, допускающего более простую интерпретацию. Сплошной линией изображено отношение йгайаДгайа для уравнения Пуассона на сетке 100 х 100 со случайной правой частью Ь.

Мы видим, что верхняя граница отражает общий характер сходимости. Семь штриховых линий— зто графики отношения йга'йа/йгайа для семи диагональных линейных систем Й,х = Ь, пронумерованных от В1 (самый левый график) до Пу (самый правый график). Все матрипы Р1 имеют ту же размерность и то же число обусловленности, что и уравнение Пуассона, поэтому, чтобы понять различия в характере сходимости, нужно присмотреться к этим матрицам внимательнее.

Матрица Б1 построена таким образом, чтобы тп1 ее наименьших и уп1 наибольших собственных значений совпадали с соответствующими собственными значениями уравнения Пуассона, а остальные п — 21п1 собственных значений были равны среднему геометрическому из наибольшего и наименьшего собственных чисел. Другими словами, 111 имеет лишь 111 = 2пт1+ 1 различных 329 б.б. Методы крыловского подпрострвкствв собственных значений. Обозначим через к; число Сб-итераций, требуемых для того, чтобы приближенное решение системы .0;х = о удовлетворило условию !!гь!!т/!!то!!г ( 10 "е.

Данные о сходимости представлены в таблице. Мы видим, что число к; шагов, требуемых для сходимости, растет вместе с числом 4 различных собственных значений. Матрица йт имеет тот же спектр, что и уравнение Пуассона, и сходимость для нее примерно столь же медленная. Мы утверждаем, что, в отсутствие округлений, методу С6 для сходнмости потребовалось бы ровно Ц = 4 шагов. Объяснение состоит в том, что можно найти многочлен дз,(с) степени 4, обращающийся в нуль на собственных значениях аз матрицы А и равный 1 в нуле, а именно П,"=,(-1 -4) П,=' (1) Согласно уравнению (6А5), результат о, шагов метода С6 минимизирует величину !!тз, !!л, = Дтз,) по всем многочленам степени 4, принимающим значение 1 в нуле.

Поскольку щ,. — один из таких многочленов и де,(А) = О, должно быть !!гл,.!!~~, = О, или тл, = О. О Один из уроков примера 6.15 заключается в том, что если число собственных значений на периферии спектра невелико (или же спектр сконденсирован на краях), то метод С6 сходится гораздо быстрее, чем показывает анализ, основанный лишь на числе обусловленности. Другим уроком является то, что поведение метода в арифметике с плавающей точкой значительно отличается от его поведения в точной арифметике.