Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 12
Текст из файла (страница 12)
Метод СЕРР— зто наиболее распространенный способ практической реализации гауссова исключения. В следующем разделе мы обсудим его численную устойчивость. Другой, более дорогостоящий способ выбора Рь и Рг указан в следствии 2.2. Этот способ почти никогда не используется на практике, хотя и существуют немногочисленные примеры, когда он дает решения хорошей точности, между тем как метод СЕРР «проваливается» (см.
вопрос 2.14). Новый способ также кратко обсуждается в следующем разделе. Следствие 2.2. Можно выбрать Р,' и Р' таким образом, чпьобьь ам был наибольшим по абсолютной величине элементом всей мапьрицы. Более общо, при вычислении ь-го столбца матрицы Е на ь'-м шаге гауссова исключения строки и столбцы с номерами от ь до п переупорядочиваются так, чтобы наибольший по абсолютной величине элемен«п образованной ими подматрицы находился в позиции (ь,ь).
Такая организация процесса называется «гауссовым исключением с полным выбором главного элеменьпа» или, для краткости, методом СЕСР. Следующий алгоритм реализует теорему 2.5, выполняя перестановки, вычисляя первый столбец матрицы Е и первую строку матрицы У и перестраивая Агг в матрицу Агг = Агг — Еьп Уьг. Мы запишем алгоритм вначале в привычных обозначениях языков программирования, а затем в обозначениях Мас1аЬ'а. Алгоритм 2.2. ЕУ-разлоэьсение с выбором главного элемента: 1ог ь = 1 со и — 1 использовать перестановки, чпьобы обеспечипъь условие аьь ф 0 (перестановки применяю«ноя и к матрицам Е и сь) 1'* так, в методе СЕРР меняются местами строки 1' и ь матриц А и Е, где ~1ауь~ — наибольший элемент, в ~А(ь': и, ь) ); в методе СЕСР меняются местами строки 1' и ь матриц 2.3.
Гвуссово исключение А и Хч а танисе сп1олбцы й и г' матриц А и У, где ~а ь~ — наибольший элемент в ~А(з: п, 1: п)! '/ /* вычислить столбец 1 матрицы Х (матрицу Хг1 в (2.10)) '/ /от у = г + 1 1о п 1ч = аз;/аа епд /ог /* вычислитаь стпроку г матприцы У (матрицу Хттг в (2.10)) '/ /от у = г 8о и иб =а,у епд /от /* обновить Агг ( ппобъг получить матрицу Агг = Ага — Хггбю в (2.10)) '/ /от у' = г' + 1 Фо п /от (г = г' + 1 1о п ав = а ь — 1,*им епд /от епд /от епд /от Заметим, что после того как 1-й столбец в А использован для вычисления 1-го столбца матрицы Ь, он нигде не используется в дальнейшем. Точно так же, 1-я строка в А не используется после того, как вычислена 1-я строка матрицы (Х.
Это позволяет нам записывать Х и (Х по мере того, как они вычисляются, на место матрицы А, так что для хранения этих матриц не требуется дополнительной памяти. При этом Х, занимает нижний треугольник в А (без главной диагонали, поскольку диагональные единицы в Х, не хранятся явно), а У занимает верхний треугольник. В результате приходим к следующему упрощенному алгоритму. Алгоритм 2.3.
Х 1)разложение с выбором главного элемента и записью Х, и У на место А: /отг'=1 1ои — 1 использоватпь перестановки (см. детали в алгоритме 2.2) /от у' = г' + 1 1о п аз = ач/аи спи /от /от у = г' + 1 1о и /от Й = г + 1 1о п аз =аль — а;вать епд /ог спи /от епд /от Использование обозначений Ма11аЬ'а приводит к такому алгоритму: Алгоритм 2.4. Х 1Х-разложение с выбором главного элемента и записью Х и У на место А: /ог1=1 1ои — 1 использоватаь перестановки ('см.
детали в алгоритме 2.2) 52 Глава 2. Решение линейных уравнений А(1 + 1 с и, с) = А(1+ 1: и, г)/А(г, с) А(1+ 1: и,'+ 1: и) = А(1+ 1: и, с + 1: п) — А(1 + 1: п, 1) * А(г, 1+ 1: и) еаза /ог В последней строке алгоритма А(1+1: п,1) «А(«,1+1: п) есть произведение матрицы размерности (и — 1) х 1 (матрицы Ьм) и матрицы размерности 1 х (и — 1) (матрицы Есш); оно представляет собой квадратную матрицу порядка и — 1. Теперь мы выведем данный алгоритм заново, отправляясь от, видимо, самого известного описания гауссова исключения: «взять по очереди каждую строку и вычитать ее кратные из последующих строк с тем, чтобы аннулировать поддиагональные элементык Прямой перевод этого описания дает следующий алгоритм: Еог1=1 Со п — 1 /' для каждой строки 1*/ Сот у' = 1+ 1 Со п /* вычесть кратное строки 1 из строки у */ Еог Ее= с Сои /*... для столбцов с номерами отздо и...'/ а,а = ась — -,л-ась /"...
аннулировать поддиагональные элементы столбца 1 '/ епс( Еог епсЕ Еог епс1 Еог Произведем теперь некоторые усовершенствования в этом алгоритме, изменяя его так, что в конечном счете он окажется идентичен алгоритму 2.3 (кроме выбора главного элемента, который мы опускаем). Обратим, прежде всего, внимание на то, что нет нужды вычислять нулевые элементы под главной диагональю, заранее зная, что они нулевые. Это сокращает цикл по Ес и приводит к такому алгоритму: Еог 1 = 1 Со п — 1 Сог у' = 1 + 1 Со и Еог й = г + 1 Со п а с а ь = а ь — -нс-ась ан епс( Еог епс1 Еог епсЕ Еог Следующее усовершенствование состоит в том, чтобы вычислять отношение а,/ан вне внутреннего цикла, поскольку оно не зависит от индекса Й этого цикла: Еог 1 = 1 Со и — 1 Еогу =1+1 Соп лч а« епсЕ Еог Еог у' = Е + 1 Со и Еог Ес = 1 + 1 Со и а ь = пи — (сань епсЕ Еог епсЕ Сог епс1 Еог 63 2А.
Анализ ошибок Наконец, мы помещаем множители 1,.; на места поддиагональных элементов а;, которые перед тем были аннулированы; последние нигде уже больше не используются. В результате получаем алгоритм 2.3 (с точностью до выбора главного элемента). Число операций в ВП-разложении можно подсчитать, заменяя циклы суммированиями в тех же пределах, а содержимое внутренних циклов — числом операций в них: й(ь +Е Е и — 1 ((ив 1)+2(н-1)') = йпг+0(па) 3 2.4. Анализ ошибок Напомним нашу двухшаговую парадигму для получения оценок ошибки в при- ближенном решении системы Ах = 6: 1. Проанализировать ошибки округлений с тем, чтобы показать: результат х решения системы Ах = 6 является точным решением возмущенной системы (А+ бА)х = Ь+ бЬ с малыми возмущениями бА и бЬ. Это пример обратного анализа ошибок, причем бА и бЬ называются обратямми ошибками.
2. Применяя теорию возмущений из равд. 2.2, получить оценку ошибки; можно использовать, например, оценку (2.3) или (2.5). В этом разделе мы преследуем две цели. Первая состоит в том, чтобы показать, как нужно реализовать гауссово исключение, чтобы обеспечить малость обратных ошибок бА и бЬ. Мы хотели бы в частности чтобы отношения 66АЬ \ ЬАЬ и были величинами порядка 0(г). Мы не можем надеяться на ббльшую ~16ЬП ~~ь|~ малость, поскольку простое округление наибольших элементов в А (или Ь), превращающее их в числа с плавающей точкой, уже может повлечь за собой соотношение )(Аа" > г (или;Гь((г > г).
Выясняется, что бА и бЬ вовсе не обяза- 66АЬ ИЬЬИ ны быть малыми, если не позаботиться о разумном выборе главных элементов. Этот вопрос обсуждается в следующем разделе. Прямая подстановка и обратный ход, выполняемые для 6 и У с тем, чтобы завершить решение системы Ах = 6, требуют 0(нг) операций, поэтому полная стоимость решения системы посредством гауссова исключения составляет -нг + 0(пг) операций. Здесь использована формула 2',.
1" = ть Ы/(Й+ 1) + 0(т"). Она достаточна для того, чтобы определить член наивысшего порядка в выражении для числа операций. Реализация гауссова исключения отнюдь не сводится к записи вложенных циклов в алгоритме 2.2. В самом деле, в зависимости от конкретных компьютера, языка программирования и порядка матрицы, простой акт перестановки двух последних циклов по 1' и Ь может изменить время работы алгоритма на несколько порядков. Мы подробно обсудим эти вопросы в равд.
2.6. Глава 2. Решение линейных уравнений Вторая цель заключается в том, чтобы вывести практичные оценки ошибок, которые были бы одновременно легко вычисляемыми и «реалистичными», т. е. близкими к подлинным ошибкам. Выясняется, что наилучшие оценки для 'ЯВАМИ, которые мы сможем обосновать формально, все же, в общем случае, много больше ошибок, наблюдаемых в реальных задачах. Поэтому наши практичные оценки ошибок (выводимые в равд.
2.4.4) опираются на вычисленную невязку г = Ах — Ь и оценку (2.5), а не оценку (2.3). Нам нужен, кроме того, недорогой способ оценивания числа к(А); эта задача обсуждается в разд. 2.4.3. К сожалению, мы не располагаем оценками ошибок, которые бы всегда удовлетворяли обеим нашим установкам на экономичность и реалистичность, т. е.
оценками, которые бы одновременно 1. могли быть вычислены с затратой числа операций, пренебрежимо малого в сравнении со стоимостью решения системы Ал = Ь (например, с затратой 0(п ) флопов, тогда как стоимость гауссова исключения составляет 0(пз) флопов); 2. всегда были не меньше подлинной ошибки н превосходили ее не более чем в заранее установленное число (скажем, 100) рвз. Практичные оценки ошибок из разд. 2.4.4 стоят 0(пз) операций, но в некоторых, очень редких случаях они оказываются чересчур малы или чересчур велики по сравнению с действительной ошибкой. Вероятность таких случаен столь мала, что эти оценки широко используются на практике. Единственный способ получить подлинно гарантированные оценки состоит в использовании интервальной арифметики, или арифметики с очень высокой разрядностью, или и того, и другого; стоимость этого способа в несколько раз выше стоимости решения системы Ах = Ь (см.
равд. 1.5). В действительности, выдвинута гипотеза о том, что оценки, удовлетворяющей обеим нашим установкам на экономичность и реалистичность, не существует; вопрос этот, однако, пока остается открытым. 2.4.1. Необходимость выбора главного элемента Давайте применим 1Л1-разложение без выбора главных элементов к матрице .0001 1 А = ' 1 1, используя трехразрядную десятичную арифметику с плавающей точкой, и разберемся в причинах, почему получается неверный результат.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
