Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 12
Текст из файла (страница 12)
Из этого основного результата можно тотчас вывести неравенства, относяшиеся к невязке и ошибке в вычисленном решении. Невязка определена как Ь вЂ” Ах,== Ех„ и, следовательно, !! Ь вЂ” Ах, ,'! = — !! Ех„~, '( ), Е (1! х,!! . В определении невязки участвует произведение Ах„, так что уместно рассматривать относительную невязку, которая сравнивает норму вектора Ь вЂ” Ах„с нормами А и ха. Из приведенных выше неравенств сразу следует, что !~ Ь вЂ” Ах. !! !! А !!!! х„)! Если А не вырождена, то с помощью обратной матрицы ошибку можно записать в виде х — х,= А '(Ь вЂ” Ах,), а тогда (!х — х„'!- !)А-'!!!Ег!!х,!!. Проще всего сравнивать норму ошибки с нормой вычисленного решении. Таким образом, относительная ошибка удовлетворяет неравенству "!! (р)!А(()!А "'Р ".
!! х,/! Оказывается, что !!А-'!)=-)ут и потому сопд(А)=)!А!! !!А '!!. з. подпгог хммы онсоме и зошз Итак, ~*) (рсопс((А)р '. //х ) Реальное вычисление числа сопб (А) предполагает знание А-'. Если а; — столбцы А, аз — столбцы А-', то для векторной норьсм, которую мы используем, сопб(А)=шах~,'а,.(Ьгпах()а~!,'.
Легко вычислить !(А~1, однако вычисление йА-Ц( примерно удваи- вает время, нужное для гауссова исключения. К счастью, точ- ное значение сопб (А) требуется редко. Обычно бывает достаточно любой разумной оценки для него. Подпрограмма РЕСОМР, описываемая в следующем парагра- фе, оценивает число обусловленности матрицы таким образом: „)г( сопб(А) ж шах(,а ( — '„ ! '' 1У) ' где у и г — векторы, определяемые подпрограммой„такие, что й4Я~у~)ж)~А-')~. Для этого решаются две системы уравнений Агу=с, Аз=у, где А' — транспоиированная для матрицы А, а е — вектор с компонентами ~Е выбираемый так, чтобы максимизировать рост н процессе обратной подстановки для у.
Зта оценка является лишь нижней границей для действительного числа обусловленности, однако есть определенные теоретические основания ожидать, что это очень точная оценка. 3.3. Подпрограммы ОЕСОМР и 5Оз.ЧЕ Почти в любой машинной библиотеке имеются подпрограммы, основанные на вариантах гауссова исключения с частичным выбором ведущего элемента для решения систем линейных уравнений. Детали реализации разных подпрограмм очень различны.
Зги детали могут серьезно отразиться на времени выполнения данной подпрограммы, но они не должны сильно влиять на ее точность, если подпрограмма правильно составлена. В этом параграфе мы опишем две такие подпрограммы— РЕСОМР и ЗОРЧЕ. ОЕСОМР выполняет ту часть гауссова исключения, которая зависит лишь от матрицы. Она сохраняет множители и информацию о ведущих элементах. БОЕВОЕ использует эти результаты, чтобы получить решение для произвольной правой части. 3. Линейные системы уРАВнений 62 РЕСОМР вычисляет также оценку обусловленности матриц Такая оценка является намного более надежной и полезной м рой близости к вырожденности, чем такие величины, как определитель или наименьший ведущяй элемент. Зта оценка дает нижнюю границу для действительной обусловленности, но она вычисляется таким образом, что почти всегда' отличается от действительной обусловленности не более чем в и раз, а обычно гораздо ближе к ней.
Иными словами, почти для всех матриц РЕСОМР определяет величину СОХО, удовлетворяющую неравенствам < СОКР ( сопд (А). В тех случаях, когда СОМР(сопд(А)/л, это число все же правильно оценивает чувствительность решений для большинства правых частей. Ошибки округлений обычно не дают возможности подпрограмме РЕСОМР или любой другой подпрограмме гауссова исключения определить, вырождена или нет входная матрица. Если в процессе исключения встретится ведущий элемент„равный точному нулю, то РЕСОМР присваивает СОКР значение !О", чтобы сигнализировать, что обнаружена вырожденность.
Число 10" находится между ~)' и ()и во всех современных плавающих системах, так что оно больше величины, обратной к машинному эпсилон, и меньше уровня переполнения. Однако появление нулевого ведущего элемента не обязательно означает, что матрица выро>кдепа, так же как вырожденная матрица не обязательно порождает нулевой ведущий элемент. В действительности самым обычным источником нулевых ведущих элементов яв.ляются ошибки в вызывающей программе.
Нужно осознать, что при наличии частичного выбора ведущего элемента любая матрица имеет треугольное разложение. На самом деле РЕСОМР работает даже быстрее, когда встречается нулевой ведущий элемент, поскольку это означает, что соответствующий столбец уже находится в треугольной форме. Единственная трудность с таким элементом состоит в том, что зОЕУЕ будет делить на него в процессе обратной подстановки.
Поэтому ЗОРЧЕ не должна использоваться, если РЕСОМР присвоила СОРТИР значение, много большее, чем ~'. Некоторые подпрограммы, имеющиеся в машинных библиотеках, включают в себя прием, известный как итерационное улучшение или итерационное уточнение.
Зто процесс, состоящий в вычислении певязки в арифметике повышенной точности и решении системы уравнений, правой частью которой является эта невязка, чтобы получить поправку к первоначально вычисленному решению. Уточненный результат зачастую имеет меньшую 3 3. гюдпеоггхммы оасомР и зосча 63 ошибку, но не обязательно меньшук> невязку. Кроме того, величина поправки дает еще одну меру чувствительности решения к ошибкам в исходных данных и прн вычислениях. .Мы решили не включать в эту книгу программу итерационного улучшения по нескольким причинам. Во-первых, для большинства приложений решение, полученное без улучшения, является удовлетворительным.
Во-вторых, ошибки во входнои информации обычно влияют на решение сильнее, чем округления в процессе вычислений. В-третьих, наша оценка обусловленности дает ту же информацию, что и величина поправки. Наконец, что, вероятно, самое важное, наличие и доступность требуемой высокоточной арифметики неодинаковы для разных машин. Универсальная программа решения линейных уравнений, которая эффективно включает в себя итерационное улучшение, не может быть написана на стандартном ФОРТРАНе. На машинах типа 1ВМ 350 в обычных вычислениях, как правило, пользуются удвоенной точностью, так что арифметика еще более высокой точности должна реализовываться специальными подпрограммами.
Чтобы прокомментировать некоторые детали в подпрограммах ОЕСОМР и ЗОРЧЕ, мы должны напомнить, как фортран- системы хранят матрицы. Если программа содержит оператор Р1МЕК51ОЫ А(3, 5), то в памяти будет зарезервировано Зхб=-15 мест для элементов матрицы А. Они будут храниться в следующем порядке: А (1,1) А (2,1) А (3,1) А (1,2) А (2,2) А (3,2) А (1,3)... Другими словами, элементы каждого столбца хранятся подряд, Элементы одной строки отделены друг от друга числом мест, равным первому индексу в операторе размерности.
Это соглашение включено в спецификации Американского национального института стандартов для ФОРТРАНа. Многие употребительные матричные операции наиболее естественно описываются в терминах строк. Например, в гауссовом исключении кратное одной строки вычитается из другой строки. Реализованные на ФОРТРАНе, такие операции обычно приводят к внутренним циклам,меняющим второй индекс массивов. Вследствие этого возможны два потенциально неблагоприятных воздействия на эффективность программы.
Вычисления индексов могут оказаться дорогостоящими, поскольку для них нужна информация, содержащаяся в операторе размерности. Операционным системам, управляющим обменом между быстродействующей и внешней памятью в ходе вычисления, возможно, придется проделать чрезмерную работу. По этим причинам мы реализонали гауссово исключение несколько необычным образом, так, з линейные системы уРАВнений что все внутренние циклы меняют первый индекс.
Такая реали зация для некоторых типов операционных систем может ока- заться значительно более эффективной, Относительно деталей см. статью Моулер (1972). Большинства диалектов ФОРТРАНа (хотя и не все) содержат, условие относительно переменньи размерностей массивов, явля-1 ющихся параметрами подпрограмм. В основную программу .* можно включить описание О1МЕЮ1ОЛ) А (50, 50), но в действительности работать с Л' х Л'-матрицей, где У меня- ется от задачи к задаче.
Подпрограммы типа ОЕСОЛ1Р и ЗО(.ЧЕ требуют задания как действительного рабочего порядка Л', так и величины 50, содержащейся в операторе размерности, посколь- ку эта величина дает приращение памяти между последовате- льными элементами строки. Такая информация о размерности в подпрограммах ОЕСОМР и 501ХЕ обозначается КР1М. Подпрограмму ОЕСОМР можно использовать для вычисле- ния определителей. Зта возможность вытекает нз трех основных свойств определителей.
Прибавление кратного одной строки к другой не меняет определителя. Перестановка двух строк изме- няет знак определителя. Определитель треугольной матрицы равен попросту произведению ее диагональных элементов. РЕСОМР использует последнюю компоненту вектора ведущих элементов, чтобы поместить туда значение +1, если было про- изведено четное число перестановок, и значение — 1, если не- четное. Чтобы получить определитель, это значение нужно ум- ножить на произведение диагональных элементов выходной мат- рицы. Неприятной особенностью вычисления определителей явля- ется то, что промежуточные произведения, а зачастую и сам определитель, имеют тенденцию быть очень большими или очень малыми числами и, следовательно, легко могут повести к пере- полнению или машинному нулю. Одна из простых предупреди- тельных мер — вычисление логарифма абсолютной величины определителя как суммы логарифмов его сомножителей.