Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 23

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 23)

Все остальные операции были выполнены точно. Каким же тогда образом смогла эта "малюсенькая" ошибка привести к столь катастрофическому отклонению вычисленного значения х1 от точного решения? Ответ на этот вопрос можно получить с помощью принципа обратного анааиэа ошибок, являющегося одной из наиболее важных концепций научного программирования. Основная идея обратного анализа ошибок состоит в выяснении не того, какая допущена ошибка, а того, какая же задача на самом деле решена.

Мы используем здесь этот принцип следующим образом. Обратите внимание, что величина 0,000001 Х 10, которая была отброшена при вычислении аЩ в (3.4.2), есть просто элемент ар, исходной матрицы. Так как это единственное место, где а22 участвует в вычислениях, то сосчитанное решение было бы тем же самым, если бы значение а., равнялось нулю. Другими словами, при вычислении на нашем четырехразрядиом компьютере мы нашли точное решение системы (3.4.4) Интуитивно мы чувствуем,что системы (3.4.1) и (3.4.4) могут иметь совершенно различные решения, и именно так и обстоит дело. Но почему же все это произошло? Виновником является большой множитель 1„, который из-за недостаточной длины машинного слова не позволил элементу а2, включиться в сумму (3.4.2).

Этот большой множитель возник иэ-за малости а„по сравнению с и2,, и "лекарством"„как и прежце, оказывается перестановка уравнений, Действительно, если бы мы на нашем гипотетическом четь1рехраэрядном компьютере решали систему (3.4.5) то имели бы следующие результаты: 0,1. 10-4 = — 05 10~, 0,2. 10' а(2',) = 0,1 10' — (-0,5 10 ') (1) = 0,1 10', Ь <~) = 0,1 . 10' — ( — 0.5 . 1О а) (0) = 0,1 . 10, 0,1 . 1О' х2=, =1,О, 0.1 .

1О' — (О,1 10')(1) х, ' — 05. 02 10' Как видно, полученное решение теперь прекрасно согласуется с точным решением. Следуя сравнительно простой стратегии, всегда можно добиться, чтобы используемые в процессе исключения множители были по абсолютной величине меньше или равны единицы. Эта стратегия, называемая стратегией частичного упорядочивания, состоит в следующем: если это необходимо, то на )с-м шаге процесса исключения осуществляется перестановка строк, так чтобы на главной диагонали оказался наибольший по абсолютной величине элемент из элементов 1с-го столбца, лежащих ниже главной диа- гонали и на ней самой.

Включая эту схему перестановок в алгоритм прямо- го исключения (3.3.11), получаем Прямое исключение с частичным упорядочиванием Дпя. )с=1,..., и — 1: найти т>Й такое,что !а,„»! = снах(1ат» ~: 1>)с1; если а,„» = О, то А вырождена, прекратить алгоритм, иначе поменять местами а»1 и а„,) (1 = Й, Й + 1,..., и), поменять местами Ь» и Ь,„. Дпя 1=1+1,1+2,..., и: (3.4.б) !с» -аи~а»», дпя у = )с + 1, )с + 2,..., и: ау+-а,у — 1;» а»;, Ь| Ь| — 1;» Ь». Метод гауссова исключения с частичным упорядочиванием зарекомендовал себя на практике как исключительно надеждный алгоритм. Тем не менее существуют два основных опасных момента, которые следует иметь в виду.

Во-первых, перед использованием алгоритма матрица должна быть масштабирована. Чтобы проиллюстрировать это, рассмотрим систему (3.4.7) Стратегия частичного упорядочивания не требует здесь никаких перестановок, так как элемент а„уже является наибольшим в первом столбце. 93 Однако если мы проведем исключение на нашем гипотетическом четырех- разрядном компьютере (см. упражнение 3.4.5), то встретимся точно с той же самой проблемой, с которой мы столкнулись при решении системы (3.4.5) .

В самом деле, (3.4.7) — это просто система (3.4.1), первое уравнение которой умножено на — 10'. Использование стратегии частичного упорядочивания обосновано в том случае, если матрица коэффициентов масштабирована так, что максимальные элементы в каждой строке и каждом столбце имеют одинаковый порядок величины. Такое масштабирование иногда называют уравковеиивдпием матрицы. К сожалению, нам не известна какая-либо простая и универсальная процедура для такого масштабирования, однако в конкретных практических задачах обычно бывает ясно, какие строки и столбцы нуждаются в масштабировании, причем оно может быть проведено до начала процесса исключения.

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

Нижняя треугольная матрица в этом случае модифицируется описываемым ниже образом. Магрицей перестановок называется такая матрица Р размера пХп, в каждой строке и в каждом столбце которой имеется ровно один элемент, равный единице, а все остальные элементы равны нулю. Пример такой матрицы размера 4 Х 4 дает матрица 1000 000 Р= 001 (3.4.8) 1 0 Перестановка строк матрицы может быть осуществлена с помощью умножения слева на матрицу перестановок. Например, умножение матрицы размера 4 Х 4 на матрицу перестановок (3.4.8) переставляет вторую и четвертую строки, оставляя первую и третью строки на месте (см.

упражнение 3.4.6) . Таким образом, требуемая стратегией частичного упорядочивания перестановка строк матрицы коэффициентов А может быть представлена как результат умножения матрицы А слева на соответствующие матрицы перестановок. Если обозначить через Р; матрицу перестановок, соответствующую перестановке строк на ~'-м шаге процесса исключения, то фактически получаем треугольную факторизацию матрицы Р„1 Р„з ° ° Р2Р, А = 1. Ю) (3.4.9) 94 а не самой матрицы А. Отсюда А = (Р 'А) К Так как произведение матриц перестановок является матрицей перестановок и обратная к матрице перестановок матрица является матрицей перестановок (см. упражнение 3.4.7), первый множитель в разложении А является перестановкой нижней треугольной матрицы, а второй множитель по-прежнему является верхней треугольной матрицей. Отметим, что если на 'г-м шаге никакой перестановки строк не требуется, то матрица перестановок Р~ является просто единичной матрицей. Перестановки строк требуют дополнительного машинного времени, а в случае ленточных матриц, кроме того, несколько усложняют хранение.

Рассмотрим сначала трехдиагональную систему. Если на первом шаге будет проведена перестановка строк, то элементы первых двух строк будут иметь вид *э*О. ~ *О... Тогда в процессе исключения в позиции (2,3) в общем случае появится ненулевой элемент, и преобразованная матрица размера (л — 1) Х (л — 1) будет опять трехдиагональной. Таким образом, эффект перестановки строк сводится к возможному появлению ненулевых элементов во второй наддиагонали преобразованной треугольной матрицы, т.е.

множитель У в разложении А больше не будет бидиагональной матрицей, а будет в общем случае иметь три ненулевые диагонали. По-видимому, самый простой путь работы с памятью в этом случае заключается в выделении дополнительного одномерного массива для хранения этих элементов второй наддиагонали. Та же самая проблема возникает и для ленточных матриц, имеющих форму, показанную на рис. 3.4.

Если на первом шаге потребуется перестановка первой строки и строки р ~ 1, то это приведет к появлению рдополнительных элементов в первой строке, что в свою очередь в процессе исключения приведет к появлению дополнительных. элементов в строках 2, 3, ..., р+ 1. Таким образом, необходимо предусмотреть место для хранения р возможных дополнительных наддиагоналей. Простейший способ заключается в выделении с самого начала дополнительного массива памяти размера пХр. Естественно, для этого потребуется лр дополнительных ячеек памяти. Другой путь основан да том, что нужное количество дополнительной памяти не превосходит объема памяти, используемого для хранения ненулевых подциагоналей.

Так как в ходе алгоритма лежащие на поддиагоналях элементы исключаются, то занимаемая ими память освобождается и ее можно использовать для хранения новых наддиагональных элементов. Однако обычно это пространство под главной диагональю используется для запоминания множителей, так что если эти множители нужны, то у нас нет никакой другой альтернативы, кроме выделения дополнительной памяти. Хотя для произвольных невырожденных матриц использование стратегии частичного упорядочивания является необходимым, имеются некоторые типы матриц, относительно которых известно, что они не требуют никаких перестановок. Наиболее важными иэ них являются диагонально доминирующие матрицы (см. уравнение (3.2.21)) и положительно опре- 95 ! 0 а, /и (3.4.11) = /!! . /и !!...

а,. „! ... а„„ /„! ... /„ Приравнивая элементы первого столбца слева и справа от знака равенства в (3.4.11), видим, что а/! =/!,/!!, так что первый столбец матрицы Ь находится по формулам /„=(ап)!/2 /!, -ап//и !-=2 и (3.4.12) Аналогично получаем соотношения ! / аи= ~ //ь, а!у= Х /!ь//а, /( /, (3.4.13) ь=! а=! которые являются основой дпя последовательного определения столбцов матрицы А по следующему алгоритму.

Характеристики

Список файлов книги