Амосов А.А., Дубинский В.А., Копченова Н.В. Вычислительные методы для инженеров (1994) (1095853), страница 34
Текст из файла (страница 34)
Заметим, что где А~ ~а 1> = Т 1,„А<а-г) 6~~1> = Т„1 $(а-г) У ! Введем обозначение зз для полученной верхней треугольной матрицы А~ а ы . Она связана с исходной матрицей А равенством Я= ТА, где Т = Т„1~ ... Тзм...ТззТпа ... Т|зТ1з — матрица результирующего вращения. Заметим, что матрица Т ортогональна как произведение ортогональных матриц. Обозначая Я = Т 1 = Тт, получаем ЯЯ- разложение матрицы й. Обратный ход метода вращений проводится точно так же, как и для метода Гаусса.
Метод вращений обладает замечательной численной устойчивостью. Для него оценка (5.43) справедлива с коэффициентом 1(т) = бтп. Однако этот метод существенно более трудоемок по сравнению с методом Гаусса. Получение ЯЯ-разложения для квадратной матрицы А общего вида требует примерно 2гпз арифметических операций. Пример 5.18. Используя метод вращений, решим на 6-разрядной десятичной ЭВМ систему уравнений 2з1 — 9з~ + 5зз = -4, 1.2х1 — 5.3999л~ + Вгз = 0.6001, — з~ — 75гз = 8 5 ° П р я м о й х о д. 1-й ш а г. Исключим з1 из второго уравнения.
Для этого вычислим с1г и ~з по формулам (5.80): 2 1.2 с1г = и 0.857493, з~г = и 0.514495. ,~21 г~ 12и,ф ~ Д2 Преобразуя коэффициенты первого и второго уравнений по формулам (5.83), приходим к системе 2.33238з1 — 10. 4957зд + 7.37444зв = — 3. 12122, 7.85493 10 за + 2.57248зз = 2.57256, х~— зз — 7.5зз = — 8 5 Далее вычислим коэффициенты с1з и з1з по формулам (5.84): 2. 33238 1 йз= и 0.919087 з~з = и 0.394055.
лн2Я3 ~. ~г л332щ5 + 1г 169 Заменяя первое н третье уравнения нх линейными комбинациями с коэффи- циентами с(з, з)з н -8(з, с)з соответственно, получим систему 2.53772х) — 10.0405хг + 3.82234хз = -6.21814, 7.85493'10 зхг + 2.57248хз = 2.57256, 3.21680хг — 9.79909хз = -6 58231. 2-й ш а г. В полученной системе имеем агг = 7.85493 10 з, азг (1! - (1) = 3.21680. Поэтому (1) ()гг сгз )з 2.44185.10, з (з 1.00000, Заменяя второе и третье уравнения системы нх линейными комбинациями с коэффициентами сгз, згз н -згз, сгз соответственно, приходим к системе 2.53772х) — 10.0405з)г + 3.82234хз = -6.21814, 3.21680хг — 9.79903хз = -6.58225, — 2 57272хз = -2 57272. О б р а т н ы й х о д дает последовательно значения хз = 1, хг = 0.999994, х) = -1.58579 10 з.
2. Метод отраисений. Матрицами Хаусхолдера 1или отражений) называются квадратные матрицы вида ~= Я~ — 2и)и)т. Здесь и) — вектор-столбец в й", имеющий единичную длину. Матрица 1 является ортогональной и симметричной. Умножение на нее называют преобразованием Хаусхолдера (или отпражениел(). Действие матрицы 1 на вектор х можно интерпретировать как ортогональное отражение вектора в Ж' относительно гиперплоскости, проходящей через начало координат и имеющей нормальный вектор, равный и). Как и вращения, отражения используются для обращения в нуль элементов преобразуемой матрицы. Однако здесь с помощью одного отражения можно обратить в нуль не один элемент матрицы, а целую группу элементов некоторого столбца или строки.
Поэтому, являясь почти столь же устойчивым, как и метод вращений, метод отражений позволяет получить ЧЛ-разложение квадратной матрицы общего вида примерно за (4/3)тз арифметических операций, т.е. в полтора раза быстрее. Изложение самого метода можно найти, например, в Щ. 170 Ь 5.11. Итерационное уточнение В большинстве случаев метод Гаусса с выбором главного элемента позволяет найти приближенное решение с довольно высокой точностью.
Однако иногда возникает необходимость найти решение с большей точностью. Полезно знать, что существует метод, позволяющий найти приближенное решение с относительной точностью, сравнимой с ем, если только число обусловленности не слишком велико. Этот метод, называемый итерационным утпочнениел, требует небольшого (примерно на 25%) увеличения машинного времени по сравнению с затратами на получение решения методом Гаусса. Пусть х<о> — найденное на ЭВМ приближенное решение системы Ах = Ь. Напомним (см. 5 5.1), что невязка г<о< = Ь вЂ” Ах<о> и погрешность е< о' = х — х< о> связаны равенством (5.87) А.<о1 = г<о< Если бы удалось найти с<о' как точное решение системы (5.87), то вектор х'" = х<® + е<о> дал бы точное решение системы Ах = Ь.
Однако в действительности вычисленное на ЭВМ значение х'" неизбежно будет содержать ошибку. Тем не менее можно ожидать, что х<'~ окажется лучшим приближением, чем х<о>. Используя приближение х' ", аналогичным образом можно найти приближение х< 2) . Опишем более подробно очередной (Ь + 1)-й шаг метода. 1о. Вычисляют г< "««Ь — Ах<~'. Исключительно важно, чтобы вычисление г'"' производилось с повышенной точностью.
Дело в том, что Ь м Ах<"' и поэтому при вычислении невязки неизбежно вычитание близких чисел, а следовательно, потеря большого числа значащих цифр. Одна из возможностей состоит в использовании для вычисления Ах< "> и Ь вЂ” Ах< <<> арифметики удвоенной точности. 2о. Вычисляют решение системы Ае< "1 = г< "1. Так как матрица А не меняется, то получение очередного приближения с использованием однажды вычисленного ЬГ-разложения матрицы А (с учетом необходимых перестановок строк) требует сравнительно небольшого числа (примерно т~) арифметических действий.
Зо. Вычисляют х<Ь" м х< "1 + е< ">. Если число обусловленности не очень велико (например, сов<1 (А) < я '), то метод довольно быстро сходится. Сходимость характеризуется постепенным установлением значащих цифр в приближениях х< "<. Если же процесс расходится, то в приближениях не устанавливаются даже старшие значащие цифры, Необходимо отметить одну интересную особенность. В процессе 171 ной формулой и~(х) в е„> 1е'"' ~/~х< "> ~. (5.88) Пример 5.19. Используя алгоритм итерационного уточнения, найдем решение системы 1. 03х> + 0.991х2 — — 2.51, О. 991х> + 0.943х2 = 2,41 (5.89) на 3-разрядной десятичной ЭВМ, считая, что режим вычислений с удвоенной точностью на ней эквивалентен использованию 6-разрядной десятичной ЭВМ. Вычислим множитель первого (и в данном случае последнего) шага прямого хода метода Гаусса: « = 0.991/1.03 и 0.962.
Так как при вычислении на 3- разрядной ЭВМ имеем 0.943 — 0.991 0.962 и -0.01 и 2.41 — 2.51 0.962 н О, то второе уравнение системы (5.89) приводится к виду -0.01х2 = О. Обратная подстановка дает, приближенное решение х< о> = 2,44, х2<~> = О. > Итерационное уточнение. 1-й шаг. Вычислим(судвоенной точностью) компоненты невязки: г<о> — 251 1 03х<о> 0991х<о> — О 00320 г< о> — 2А1 0 991х< о> 0 943х2<о> — 0 00840 Вычисляя методом Гаусса 'решение системы 1,03е< о> + 0.991е2<о> — 0 00320 0.991е< о> .+ 0.943е2<о> — О 00840 получим е<> о> = -ОА81, е< о> = 0.496. Завершается 1-й шаг вычислением х> — х<о> + с<о> х~ = х<о> + е<о>, приводящим на 3-разрядной ЭВМ к значениям х'<> = 1.96, х'<> = 0.496, 2-й ш а г дает значения г<>>> = -0,000336, г<2<> = -0.000088, е<>> = 0.0223, е< > = — 0.0235 и х'2> = 1.98, х2<~> = 0.473.
3-й ш.а г дает значения г<>т> = 0.00186, г'2> = 0.00178, е<>2> = 0.00084, е2<2> = 0.001 и х<з> = 1.98, х<з> = 0.474. Сравнивая х<2' и х< э>, замечаем, что последние значащие цифры практи- 172 итерационного уточнения невязки х< "> обычно не уменьшаются, а даже несколько возрастают по величине. Кроме того, как оказывается, приближения х<"' могут быть использованы для грубого по порядку оценивания естественного числа обусловленности и~(х). Для этого достаточно воспользоваться приближен- чески установились и, следовательно, процесс следует завершить, приняв х~ и ю 1.98, х2 й 0.474.
Использование формулы (5.88) с ам = 5 10 4 дает следующую оценку естественного числа обусловленности: и (х) ~ вм ~ е1о) ~[ У ~ х1 о1 ~ м 2. 10з,0.496/1 96 р 407 Приведем для сравнения действительные значения решения и естественного числа обусловленности: х1 м 1.9812, т2 м 0.4735, иЬ(з) в 273 (см. пример 5.5). 3 5.12. Дополнительные замечания 1. Более подробное и в то же время доступное изложение рассмотренных в Ътой главе вопросов можно найти, например, в [38], [39], [67], [87].
2. Предположим, что требуется решить систему и линейных алгебраических уравнений с т неизвестными, где т < и. Так как число уравнений превышает число неизвестных, то вполне вероятно, что рассматриваемая система не имеет решения. Хотя уравнения системы нельзя удовлетворить точно, можно попытаться удовлетворить их как можно точнее, минимизируя величину вектора невязки г = Ь вЂ” Ах . Выбор в качестве минимизируемой величины п евклидовой нормы невязки ~~г~2 — Е ] г;~2 приводит к методу каи иеньших 1=1 квадратов решения переопределенных систем линейных уравнений [50].
Отметим, что в гл. 11 метод наименьших квадратов рассматривается в связи с задачей приближения таблично заданной функции. 3. Интересной модификацией метода Гаусса является а41оритм Краута, основанный на изменении порядка выполнения арифметических операций. Существуют два случая, когда эта модификация может иметь преимущества. Во-первых, алгоритм Краута выгоден тогда, когда для вычислений используется микрокалькулятор, так как алгоритм позволяет избегать выписывания промежуточных результатов. Во-вторых, он приводит к меньшей вычислительной погрешности при условии, что вычисления ведутся на ЭВМ, обладающей сравнительно быстрой арифметикой с расширенной точностью или особенно быстрой операцией вычисления скалярных 'произведений (х, у) = Е т,у;. 1=1 Обсуждение модификации Краута можно найти, например, в [67], [87].