Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравнений. Под ред. А.А.Абрамова (1986) (1095855), страница 22
Текст из файла (страница 22)
Объем работы, выполняемой на этапах 2 и 3, только тогда приближается к объему работы на шаге факторизации (по крайней мере для заполненных матриц), когда значение т становится близким к и. Отметим, что при вычислении А ', когда т = и, полное количество операций по-премснему имеет порядок 0(п ), т.е, тот же порядок, что и при решении линейной системы. Отметим также, что при реализации схемы (3.3.25) этапы 1 и 2 можно выполнять либо одновременно, модифицируя правые части в ходе процесса гауссова исключения, либо сначала провести факторизацию и запомнить множители 1;;, а затем перейти к этапу 2.
Все алгоритмы этого раздела основывались на предположении, что элементис с и все последующие используемые в процессе исключения диаго- в6 нальные элементы преобразованных матриц не обращаются в нуль. При решении практических задач этого оказывается недостаточно; эти элементы должны быть в некотором смысле достаточно велики, так как в противном случае могут возникнуть большие ошибки округления. Связанные.
с этой проблемой вопросы, а также модификации, необходимые для того, чтобы процесс исключения оказался жизнеспособной процедурой, будут рассмотрены в разделе 3.4. Дополнительные замечания и ссылки 3.3 Имеется большое количество книг, посвященных численным методам линейной алгебры, и в каждой из них рассматривается проблема решения систем линейных уравнений; смотрите, например, [78, ВВ). Более детальное и глубокое исследование дано в книгах (84,91 ~. Сведения о программном обеспечении решения линейных уравнений можно найти в дополнительных замечаниях к разделу 3.5. Хотя метод гауссова исключения и является чрезвычайно эффективным.
существуют алгоритмы с лучшими показателями по числу операций как функции от л. Так. например, в методе, предложенном в статье 1751, число умножений, необходимое для решения линейной системы порядкал, естьО(л ' '' ). Однако постоянный 2,В... множитель, стоящий' перед этим членом, значительно превосходит множитель при старшем члене в методе гауссова исключения. В силу этого, а также большой сложности упомянутый метод не смог составить серьезной конкуренции в практических расчетах алгоритму гауссова исключения. УПРАЖНЕНИЯ 3.3 3.3.!.
Решите вручную следующую систему: 3х, +'х, + 2х, = 1, 4.т, + 3х, + 2хз = 2, х, +4х, +2х, =4. 3.3.2. Выпишите явно матрицы ь н К факторизующие матрицу А иэ упражнения 3.3.1. Пересчитайтс решение снова, используя А и (/. 3.3.3. Покажите справедливость следующих утверждений относительно числа операций в методе гауссова исключения: а) число сложений (умножений) для вычисления новой правой части раино л1л — 1)/2; б) для вычисления множителей /;я число делений равно л(л — 1)/2; в) число делений при обратной подстановке равно л; г) число сложений (умножений) при обратной подстановке равно и (л — 1) /2.
3.3.4. Покажите, что шаг процесса исключения, обнуляюший элементы,первого столбца матрицы А, находящиеся в строках 2, 3,..., л, эквивалентен умножению А слева на матрицу 1 — !,, 1 О /л1 Локажите более общее утверждение, чтоЬл 1... 1.,1,,А = К где 1 0 1 /1+1 А 89 3.3.5. Убедитесь, что произведение матриц ь и У из (3.3.20) является трехдиагональной матрицей. 3.3.6.
Составьте программу, реализующую алгоритм гауссова исключения для трехдиагональных систем. Для хранения матрицы используйте одномерные массивы. 3.3.7. Составьте программу, реализующую алгоритм гауссова исключения для ленточной матрицы ей диагоналями выше главной диагонали и р диагоналями ниже главной диагонали. Для хранения матрицы используйте схему, представленную на рис. 3.5, 3.3.8. Докажите следующие (уже использованные нами) свойства определителей: а) определитель треугольной матрицы равен произведению ее диагональных элементов; б) определитель произведения двух матриц равен произведению определителей зтих матриц.
3.4. Перестановки а,, (( — () а;„ Π— () пи О О-1) ()-1) а„з (( — 1) и что аы = О. Если все остальные элементы 1-го столбца, лежащие ниже аг(, равны нулю, то эта матрица является вырожденной (см. упраж(( — 1) пение 3.4.2). Преобразованная матрица получена из исходной с помощью не изменяющей определителя операции добавления к одной строке матрицы другой строки, предварительно умноженной на некоторый множитель, и операции перестановки строк, которая изменяет только знак определителя.
Таким образом, если преобразованная матрица вырождена, то, вопреки предположению, и исходная матрица является вырожденной. Следова((- з) тельно, по крайней мере один из элементов иат (7с = 1 + 1,..., л) не равен нулю, и мы можем переставить строку, содержащую этот элемент, 90 Рассматривая в предыдущем разделе алгоритм гауссова исключения, мы предполагали, что а,, и все последующие делители отличны от нуля. Теперь опишем модификацию этого алгоритма, которая позволяет, проводя при необходимости перестановку уравнений, избавиться от указанного ограничения.
Предположим, как обычно, что матрица коэффициентов А является невырожденной. Допустим, что а,, = О. Тогда в первом столбце матрицы А обязательно найдется другой отличный от нуля элемент, так как в противном случае матрица А оказалась бы вырожденной (см. упражнение 3.4.1). Пусть, например, аа, чь О. Поменяв теперь местами первое и )с-е уравнения, мы, очевидно, получим эквивалентную систему, у которой элемент с индексами (1,1) отличен от нуля, так что можем выполнить первый шаг процесса исключения. Такую перестановку уравнений можно осуществить и на любом последующем шаге, если вычисленный диагональный элемент, который должен служить в качестве делителя, обратится в нуль. Действительно, предположим, что мы довели процесс исключения до стадии с 1-й строкой, обеспечив тем самым отличие от нуля нового элемента с индексами (1, ~').
Подчеркнем еще раз, что перестановка строк не изменяет решения системы. Мы уже отмечали, что перестановка строк приводит к изменению знака определителя матрицы коэффициентов. Поэтому если в нашу задачу входит вычисление определителя, то необходимо зафиксировать, какое было сделано число перестановок: четное или нечетное. Так как в точной арифметике описанный процесс гауссова исключения дает решение линейной системы за конечное число шагов и с этим процессом не связано никакой ошибки дискретизации, то единственное обстоятельство, которое может повлиять на точность вычисленного решения,— это наличие ошибок округления. Здесь имеются две возможности. Первая связана с накоплением ошибок округления при выполнении большого числа арифметических операций.
Так, например, если и = 1000, то, как показано в предыдущем разделе, для решения системы потребуется порядка л = 10 операций. Лаже если ошибка на каждой отдельной операции мала, полная погрешность могла бы быть весьма значительной. Позднее увидим, что это потенциальное накопление ошибок округления оказывается не столь серьезным, как можно было бы ожидать. Вторая возможность связана с возникновением катастрофических ошибок округления.
Если для алгоритма характерно это неприятное свойство, то такой алгоритм называют численно неустойчивым и его нельзя использовать в качестве универсального метода. Хотя описанный выше процесс перестановок с математической точки зрения. гарантирует применимость алгоритма гауссова исключения для решения любой невырожденной системы, этот алгоритм может привести к возникновению катастрофических ошибок округления и является численно неустойчивым. Чтобы понять, как это может произойти, мы проанализируем простой пример размера2Х2. Рассмотрим систему (3.4.1) с точным решением х, = — 0,4999975..., хэ = 0,999995. Предположим теперь, что в нашем распоряжении имеется компьютер.
выполняющий операции над десятичными числами, причем длина слова соответствует четырем цифрам. Это означает, что числа представляются в форме О. * е * е Х 10~. Лавайте выполним алгоритм гауссова исключения на этом гипотетическом компьютере. Во-первых, заметим, что а,, Ф 0 и никакой перестановки не требуется. Множитель — 0,2 Х 10' 1эь= ~ = — 02Х10 0,1Х10~ вычисляется точно, а значение нового а~э получается как а~~э =0,1 1О' — ( — 0,2 10 )(0,1. 10')=0,1 10' +0,2 10~ =0,2. 10 . (3.4.2) Точная сумма в (3.4.2) равна, разумеется, 0,200001 10', но, так как в слове компьютера помещаются только четыре цифры, это число должно быть представлено в виде 0,2000.10 . Это первая ошибка округления, возникшая в ходе вычислений, Новое Ь~ есть ь201= .(-о,г- ю')(о,1 ю')=о,г. ю'.
(3.4.3) Отметим, что при вычислении этого значения ошибки округления не возникают. Не возникают они и при обратной подстановке: ь|01 ог ю' — — — — =О1 ~О', а '1121 0,2 10' 0,1. 10' — 0,1. 10' х, — — О. 0,1.10 4 Найденное значение х2 прекрасно согласуется с точным значением, но значение х1 не имеет ни одного верного знака. Подчеркнем, что была допущена только одна ошибка при вычислении а 22, причем эта ошибка произошла (т) в шестом десятичном знаке.