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

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

В разд. 6.4 даны сводка итерационных методов нз этой главы и сравнение производительности (почти) всех этих методов при решении модельной задачи. 279 6.2. Интернет-ресурсы для итерационных методов Последовательность изложения методов в следующих пяти разделах примерно соответствует возрастанию их эффективности в применении к модельной задаче.

В равд. 6.5 описаны основные итерационные методы, а именно методы Якоби, Гаусса — Зейделя, последовательной верхней релаксации и их варианты. В разд. 6.6 рассмотрены методы крыловского подпространства, причем наибольшее внимание уделено методу сопряженных градиентов. Быстрое преобразование Фурье и его использование для решения модельной задачи составляют содержание разд. 6.7. Блочная циклическая редукция описана в разд. 6.8. Наконец, в равд. 6.9 обсуждается многосеточный метод, являющийся нашим самым эффективным алгоритмом для модельной задачи. В многосеточном методе на одно неизвестное приходится оптимальное количество работы 0(1). В равд. 6.10 изучается техника декомпозиции области, представляющая собой семейство приемов комбинирования простых методов, описанных в предыдущих разделах, с целью решения более сложных задач, чем наша модельная задача.

6.2. Интернет-ресурсы для итерационных методов Дпя уравнения Пуассона мы составим короткий список численных методов, которые явно превосходят все остальные из обсуждаемых нами. Однако для других линейных систем не всегда ясно, какой метод является наилучшим (вот почему мы говорим о столь многих!). Чтобы помочь пользователям выбрать среди множества существующих наилучший метод для их линейных систем, в ХЕТЫВ/1ешр1а1еа помещена справочная система. В этой директории содержатся небольшая книга [24] и программы большинства итерационных методов, обсуждаемых в данной главе.

Имеются версии книги на языках РоэсВсг1рс (см. ХЕТ11В/1ешр1асеэ/сешр1а1ее.ре) и Нурег1ех1 Маг1спр 1 апбпабе (НЕТЫВ/сешр1аСез/Тешр1асег.Ьсш1). Программы написаны на языках Мас1аЬ, фортран, С++. Слово 1етр1а1ег (шаблоны) указывает, что в описываемых книгах и программах детали, связанные с представлением матриц, отделены от алгоритмов самих по себе. Например, в методах крыловского подпространстеа (см. разд. 6.6) все, что требуется от матрицы А, это возможность умножить ее на произвольный вектор г.

Способ представления А определяет наилучший способ вычисления таких произведений, но никак иначе не влияет на организацию алгоритма. Иными словами, в шаблоне операция матрично-векторного умножения рассматривается как «черный ящик». Конкретная реализация такого черного ящика предоставляется пользователю. Аналогичный набор шаблонов разрабатывается в настоящее время для задач на собственные значения. Среди других недавно изданных учебников по итерационным методам назовем [15, 136, 214]. В самых трудоемких практических задачах, проистекающих из более сложных дифференциальных уравнений, чем наше модельное уравнение, линейная система Ах = 5 должна быть «предобусловлена», т.

е. заменена эквивалентной системой М 'Ах = М ~5, которую в том или ином отношении легче решить. Эта тема подробно обсуждается в равд. 6.6.5 и 6.10. Реализации (включая и параллельные) многих приемов предобусловливания помещены в Интернете в 280 Глава б. Итерационные методы для линейных систем вцде пакета РЕТБс (от РоггаЫе Ехгепз1Ые Тоо!)115 1ог Яс)епибс сошрпсапй) по адресу Ьир://ачачач.шсз.ап!.6оч/ресзс/регзс.Ыш1 [232]. 6.3. Уравнение Пуассона 6.3.1. Одномерное уравнение Пуассона ди(х) а1х иа иа-1 в=(а-.а)ь Ь аго(х) 5(х иа 51 иа х=(а-г.а)а Ь разность на Ь, приходим к центрироеанному Вычитая зти приближения и деля разностному приближению с)го(х) дхг 2иг — ог 1 — о+1 Ьг (6.2) Можно показать, что погрешность аппроксимации т; равна 0(Ь~ .

[[~ф[], ). Теперь уравнение (6.1) в точке х = хг можно записать как — ие 1 + 2оз — иеч.1 — — Ь /1 + Ь тп г г где 0 < 1 < )а'+1. Поскольку из граничных условий следует, что ио = и)ч 11 = О, мы имеем ))а уравнений с )ч' неизвестными и1, ..., ила: 2 — 1 0 — 1 0 — 1 2 /1 т1 (6.3) 1 Они называются условиями Дирнзяе. Возможны и другие типы граничных условий.

Мы начинаем с одномерной версии уравнения Пуассона 11ги(х) — г =/(х)а 0<х<1, (6.1) сах где /(х) — заданная функция, а и(х) — неизвестная функция, которую мы хотим вычислить. Эта функция должна еще удовлетворять граничным условиям' о(0) = и(1) = О. Мы заменим эту задачу разностной, стремясь к тому, чтобы найти приближенные значения решения в )аа + 2 равноотстоящих точках хо находящихся между 0 и 1: хг = 1Ь, где Ь = 1 и 0 < 1 < )ч'+ 1. Будем пользоваться сокращенными обозначениями иг = и(хе) и /1 = /(хг). Чтобы превратить дифференциальное уравнение (6.1) в систему линейных уравнений относительно неизвестных о1, ..., ии, аппраксимируем производные конечными разностями: 281 6.3.

Уравнение Пуассона 2.6 1.6 ол о О 2 4 6 6 10 12 14 16 16 26 Рис. 6.1. Собственные значения матрицы Тш. или Т„и =)121'+ аозт. (6.4) Чтобы решить эту систему, пренебрежем величиной т, которая мала по сравнению с у. Тогда получим Тнб = Ьз~. (6.5) (Ошибка и — 6 будет оценена позднее.) Матрица коэффициентов Тгг играет центральную роль во всем дальнейшем, поэтому мы исследуем ее несколько подробней. Сначала мы вычислим ее собственные значения и собственные векторы. С помощью тригонометрических тождеств легко проверяется следующая лемма (см. вопрос 6.1). Лемма 6.1. Собственнь.ми значениями магприцы Тгг являются числа Л = 2(1 — соз-~~д).

Собственные векторь1 г; имеют компоненгпв1 г (и) + ,~~„61п(Зйк/(Х+ 1)); при этом 2-норма вектора и. равна единице. Пусть о = [гг,...,гвг] — ортогональная магприца, составленная по столбцам иэ собствениыт векторов, и пусть Л = 111а6(Л1,..., Лиг). Тогда справедливо соОГПНОизение ТЬГ = ЯЛЯт. На рис. 6.1 изображены собственные значения матрицы Твг для 1ч' = 21. Наибольшим собственным значением является число Л1е = 2(1 — соз к ~ ) гг41 4. Наименьшее собственное значение1 есть Л1, где для мальгх 1 Л1=2 1 — соэ 2 1 — 1— 1 Обратим внимание, что, в отличие от гл, б, мы обозначаем через Лн наибольшее, в через Л1 — наименьшее собственное значение. 282 Глава 6. Итерационные методы для линейных систем Собственный вектор 5 1 Собственный вектор 1 -1 0 5 10 15 20 Собственный вектор 11 1 -1 0 5 10 15 20 Собственный вектор 2 1 0 5 1О 15 20 1 Собственный вектор 21 0 5 10 15 20 Собственный вектор 3 -1 0 5 10 15 20 Рис.

6.2. Собственные векторы матрицы Тм. -1 0 5 10 15 20 Таким образом, матрица Тм положительно определена и ее число обусловленности Лн/Лг составляет для больших Ю 4(М+ 1)г/в.г. Собственными векторами являются синусоиды, изображенные на рис. 6.2; наименьшая частота соответствует у = 1, а наибольшая — 1 = Ьг. Теперь мы знаем достаточно для того, чтобы оценить ошибку, т. е. разность между решением системы Тмй = Ьг/ и точным решением и дифференциального уравнения.

Вычитая (6.5) из (6.4), получаем о — о = ЬгТ„'т. Переходя к нормам, находим р~+ цг 6ео йо — ойг < Ь 6Тм бгйтйг Ь втаг — О(Ь бтпг) = О Ь е поэтому ошибка о — о стремится к нулю пропорционально Ьг при условии достаточной гладкости решения о (т. е. конечности величины )(2;-~(! ). Начиная с этого момента, мы не будем делать различия между о и его приближением 6, что упростит запись системы до Тчо = Ьг/.

Оказывается, что не только решение линейной системы Ь гТми = / аппроксимирует решение дифференциального уравнения (6.1), но и собственные значения и собственные векторы матрицы Ь гТм аппроксимируют собственные значения и собстлвеммые функции дифференциального уравнения. Говорят, что Л; есть собственное значение, а ге(х) — собственная функция уравнения (6.1), если ~(гх;(х) — — Лев;(х) и т;(0) — те(1) — О. 283 6.3.

Уравнение Пуассона 6.3.2. Двумерное уравнение Пуассона Теперь мы обратимся к двумерному уравнению Пуассона дев(х у) Отв(х у) (6.6) Охз Оут на единичном квадрате ((х,у): 0 ( х,у ( 1) с условием е = 0 на границе квадрата. Для дискретизации уравнения берется сетка с узлами (х;,у ), где х, = 16, у = 76 и 6 = + . Используются сокращенные обозначения е; е(16,76) и ~О = ~(16, у6); это иллюстрируется рисунком для Ю = 3.