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

Если 1 = з, то, суммируя геометрическую = О. = Фт, поэтому Ф = Ф*, и нужно лишь показать, — — и — „, ь~ ьа Б = Ее=о 4!ЙФьд = ~ э=о ы ы = Еь=о ы то эта сумма, очевидно, равна и. Если же 1 ф з, прогрессию и учитывая, что ыо = 1, находим П 338 Глава б. Итерационные методы дли линейных систем В алгоритме 6.13 показано, что для решения дискретного уравнения Пуассона Тетей'+ РТгг = й~г относительно T требуется возможность вычислять произведения с участием Дг х гУ-матрицы Я, где 2, х(у+ 1)(й+1) — яп Ж+ 1 г'т'+ 1 1 (Напомним, что в этом разделе строки и столбцы нумеруются от 0 до йг — 1.) Рассмотрим теперь ВРТ-матрицу Ф порядка 2йГ+ 2. Ее элемент Ц, й) равен / — 2хг1й'1 / — ягуй'1 куй, х)й Таким образом, Х х Ю-матрица Я с точностью до множителя — ~/ ~~, совпадает с мнимой частью подматрицы, образованной строками и столбцами матрицы Ф, имеющими номера с 2 до Ю + 1.

Поэтому, если мы найдем способ эффективного умножения на Ф с помощью РРТ, то сможем эффективно умножать и на т. (Для достижения наибольшей эффективности описываемый ниже алгоритм РРТ модифицируют таким образом, чтобы можно было умножать на Я непосредственно. Эта модификация называется быстрым синуспреобриоеаниелг. Но можно пользоваться и стандартным алгоритмом РРТ.) Итак, для быстрого вычисления произведения ЯР нужно применить операцию типа РРТ к каждому столбцу матрицы Р, а для произведения г л аналогичная операция должна быть применена к каждой строке Р. (В трехмерном случае символ И обозначал бы массив неизвестных размера Х х Ю х Ю и одна и та же операция применялась бы к каждому из Занге сечений, параллельных координатным осям.) 6.7.2.

Решение непрерывной модельной задачи посредством рядов Фурье Мы возвращаемся к нумерации строк и столбцов матриц числами от 1 до Л. В этом разделе будет показано, что алгоритм решения дискретной модельной задачи является естественным аналогом метода рядов Фурье для решения исходного дифференциального уравнения (6.1).

Это будет сделано для одномерной модельной задачи. Напомним, что уравнение Пуассона — Ц = /(х) рассматривается на отрезке [О, Ц и имеет краевые условия о(0) = о(1) = О. Чтобы решить эту задачу, разложим о(х) в ряд Фурье: о(х) = 2 ~~ сг эшЦ~гх). (Отсутствие косинусов в этом разложении объясняется краевым условием о(1) = 0.) Подставляя о(х) в уравнение Пуассона, находим гту(г' х ) вш(ргх) = /(х). Умножим обе части на яп(йхх), проинтегрируем в пределах от 0 до 1 и используем тот факт, что Д яп(ухх) э!п(йях) г(х равен нулю при у ~ й и 1/2, 1 340 Глава 6.

Итерационные методы для линейных систем адхь и Ь(х) = ~ ~: Ььхь степени Х вЂ” 1. Тогда их произведением является многочлен с(х) = а(х) Ь(х) = 2 „„сдхь, где коэффициенты со, сгдг д определяются с помощью дискретной свертки. Одной из задач, решаемых посредством преобразований Фурье, рядов Фурье или ПУТ, является превращение свертки в произведение. В случае преобразования Фуоье имеем У(1 * д) = У(1) У(д), т. е. образ Фурье свертки есть произведение образов Фурье.

Для рядов Фурье справедливо соотношение сз(1 ед) = с,(З) с,(д), т. е. коэффициенты Фурье свертки суть произведения коэффициентов Фурье. То же самое верно для дискретной свертки. Теорема 6.10. Пусть а=(ар,...,аи д,О,...,0)~и 6=(Ьо,...,Ьгд — д,О,...,0] векторы размерности 2ддд и пусть с = а е Ь = [со,..., сгы д)т. Тогда (Фс)д = (Фа) д (ФЬ) ь . Доказательство. Если а' = Фа, то ад — — 2 "о а до"з, т. е.

ад, равно значению э=о многочленаа(х) = ~С1:о а хд в точке х = до". Точно также, из Ь' = ФЬ следует, дч-д что Ьь — — ~,:о 6 ыдз = Ь(од~), а из с' = Фс следует, что с'„= 2 ~~~о ' с до"д = с(ы"). Поэтому а1 Ь~ь — — а(до~) Ъ(до~) = с(ы~) = сд„ что и требовалось. Иными словами, РРТ эквивалентно вычислению значений многочлена в точках ыо,..., ьдддд '. Обратно,ПЬЕТ эквивалентно интерполированию много- члена, т.

е, вычислению его коэффициентов по значениям многочлена в точках ЮО дои-д 6.7.4. Вычисление быстрого преобразования Фурье Мы выведем алгоритм РГТ, пользуясь только что полученной интерпретацией РРТ как операции вычисления значений многочлена. Наша цель — вычислить значения многочленаа(х) =;> „' аьхь в точках х = ыз для 0 < у < Ю вЂ” 1. Для простоты, предположим, что Х = 2 .

Теперь имеем а(х) = ао + ад х + адхг + .. + аы дхддд ' =(ао+агх +адх~+ )+х(ад+агх +авх +.-.) — аеееь(х ) + х ' аоедед(х ). Таким образом, нужно вычислить значения двух многочленов а,„,„и а,гв степени г — 1 в точках (ыд), 0 < д < Ю вЂ” 1. В действительности, можно ограничиться — точками до д, отвечающими значениям 0 < д < — — 1, так как ьд г' ду ,„гз — ыг0ь е д Итак, вычисление многочлена степени 1де — 1 = 2 — 1 во всех Х корнях Л-й степени из единицы равносильно вычислению двух многочленов степени — — 1 дде 2 во всех — корнях степени — из единицы с последующим комбинированием дч ддд г г результатов, требующим Х умножений и сложений.

Эту связь можно использовать рекурсивно. 341 6.7. Быстрое преобразование Фурье Алгоритм 6.14. РРТ (рекурсивная версия): /ипс1зоп РРТ(а,Ю) ЧА =1 возвратить а е(зе а',„,„= РРТ(а,„,„, Х/2) а',зз — — РРТ(а зз, 1х'/2) — 2з(/Ю [ о коз-~) возвратить а' = [а',„,„+ и. * а'зз,а,'„,„— ш. з а'зз[ жй/ Здесь символ .ь обозначает покомпонентиное умножение массивов (как в МаИаб'е) и используется соотношение ыгь~~~ = — ое'. Обозначим сложность этого алгоритма через С(Х).

Мы видим, что функция С(Х) удовлетворяет рекуррентному соотношению С(Ю) = 2С(Л/2)+ЗЮ/2 (при этом предполагается, что степени ы вычислены и записаны в таблицы заранее). Решение этой рекурсии дает выкладка С(У) = 2С вЂ” + — = 4С вЂ” + 2 — = 8С вЂ” +3 ЗЮ = 108з ь7 2 Применение РРТ к каждому столбцу (или каждой строке) матрицы порядка Х обходитсЯ, таким обРазом, в 1ойз Х вЂ” опеРаций. Этот анализ сложности зло объясняет элемент табл. 6.1, относящийся к РРТ.

В практических реализациях РРТ для достижения наибольшей эффективности используются простые вложенные циклы, а не рекурсии. Кроме того, в этих реализациях компоненты результата иногда выдаются в двоичноинверсном порядке. Это означает, что в вьгходном векторе у = Фх вместо обычного порядка компонент уо, уз, ..., ум з принята последовательность, отвечающая обращению порядка битов в двоичном представлении индексов. Например, при Ю = 8 индексы изменяются от 0 = 000з до 7 = 111з. В таблице показаны обычный порядок индексов и двоична-инверсный порядок.

В обратном РРТ происходит возврат от модифицированного порядка к исходному. Следовательно, эти алгоритмы могут быть применены к решению 342 Глава б. Итерационные методы для линейных систем модельной задачи; нужно лишь, чтобы деление на собственные значения про- изводилось с учетом двоично-инверсного порядка индексов. [Заметим, что в Маз1аЪ'е результаты всегда выдаются в обычном порядке возрастания ин- дексов.) 6.8. Блочная циклическая редукция А — 1 — 1 — 1 А и предположим, что порядок М матрицы А = Т22+ 21м нечетен. Отметим еще, что хз и Ь, суть Ю-векторы.

Рассматривая три подряд идущие блочные строки, применим блочный метод Гаусса: + [ — хз — г +Аз [ ) ). =ь = Ь. =ь,+, + Ах — х х' 3 + Ах. Х' 4 — Х 4.1 + АХ,4.1 Это пРивеДет к исключению вектоРов хз 1 и хть1. — х, 2 + (А — 21) хз — хз~ьг = 61 1 + АЬ, + Ьутз. Если проделать это для каждой группы из трех подряд идущих блочных строк, то получим две системы уравнений.

Одна из них соответствует векторам х. с четными индексами 1:  — 1 — 1  — 1 61+ Аьг+ ьз Ьз + АЬ4 + Ьз Х2 Х4 (6.50) хм Ьз1-2 + Аьз1-1+ Ьгч — 1 — 1 В Здесь В = Аг — 21. Другая система соответствует векторам х с нечетными индексами: 61 + Х2 Ьз + х2 + х4 Х1 хз 61ч + Х24-1 Блочная циклическая редукция — зто еще один быстрый метод (т.

е. метод со сложностью 0(11'2!окг 1т )) для решения модельной задачи, однако область его применимости несколько шире, чем метода, основанного на РРТ. Наиболее быстрое решение модельной задачи на векторных компьютерах часто достигается комбинированием блочной циклической редукции и РРТ. Вначале будет описана простая, но численно неустойчивая версия алгоритма; затем мы кратко обсудим, как стабилизировать ее. Запишем модельную задачу в виде 343 6.8. Блочиая ииклическая редукция Она может быть решена после того, как из уравнения (6.50) определены хз с четными индексами. Заметим, что уравнение (6.50) имеет ту же форму, что и исходная задача, поэтому процесс может быть повторен рекурсивно.

Так, на следующем шаге мы получили бы системы С вЂ” 1 — 1 С вЂ” 1 [" ] где С = Б2 — 21, — 1 — 1 С Алгоритм 6.15. Блочная циклическая редукция: 1) Приведение1 1ог т = 0 со й — 1 А(г+1) (А(г) ) 2 21 1ог 2 = 1 со Х .~1 + АБО 5~ ~ + Ь~~~ 1 — 21-1 21 21+1 епд 1ог еаза 1ог Комллентарии1 после г шагов задача приводится к виду: А00 — 1 — 1 — 1 А1 "1 2) Система А1Я1х1Я) = 51а1 реи1ается другим методом. 3) Обратная подстановка: уатт = й — 1,...,0 Будем продолжать процесс до тех пор, пока не останется только одно (блочное) уравнение.

Оно будет решено иным методом. Мы формализуем этот алгоритм следующим образом: предположим, что Х = Жв —— 2"+1 — 1, Ф„= 2"+' " — 1. Положим А1в) = А и Ь~ ~ = Ь для 2 = 1,..., Ю. 344 Глава б. Итерационные методы для линейных систем ~ог 1 = 1 (о Ю«+г х21 хд епа 1ог гог г = 1 го )(г«з(ер В решить систему А(")х(') = Ь(") + х("~д + хт'+д («) относительно х( ) (полагаем хо = хы = 0) («) (1 ) епд 1ог епд 1ог Вектор х = х(о) есть искомый результата. У этого простого метода есть два недостатка: 1.