Этюды для программистов [неполностью, главы 1–24]

Как вычислить обратную матрицу? Следующий алгоритм принадлежит Гауссу.

Во-первых, положите матрицу X равной матрице I. В процессе вычислений матрица А будет в конце концов преобразована в I, матрица X, которая изначально была единичной матрицей, станет обратной к матрице А.

Для каждого столбца А, начиная со столбца 1 слева и кончая столбцом n справа, выполните следующее: Обозначим столбец, который будет обрабатываться на каждом этапе, символом j.

Пусть есть наибольший по абсолютной величине элемент в столбце j ниже строки j − 1. Если М равно нулю, то А — вырожденная матрица и продолжать обращение не имеет смысла. В противном случае поменяйте местами в обеих матрицах А и X строку j и строку, в которой находится М. И наконец, разделите каждый элемент в строке j матриц А и X на новое значение A_jj.

Теперь для всех строк i, i ≠ j, выполните все вычитания:

A_ik = A_ik − A_ijA_jk, j ≤ k ≤ n,

X_ik = X_ik − A_ijX_jk, 1 ≤ k ≤ n.

Результатом всех этих поэлементных вычитаний будет вычитание из строки i строки j с коэффициентом A_jj в обеих матрицах А и X. После выполнения этого шага для всех j все элементы сверху и снизу от A_jj станут нулевыми, а сам элемент A_jj будет равен единице. В матрице А не нужно выполнять вычитания слева от столбца j, поскольку все элементы строки j слева от A_jj равны нулю.

Для любого алгебраиста будет одно удовольствие доказать, что этот алгоритм всегда работает правильно и после его остановки X = А−1, если матрица А невырожденна. Вы едва ли найдете алгоритм, более приспособленный для структурной реализации. Почему бы нам, исключительно ради забавы, не провести небольшую проверку? Матрица Гильберта Hⁿ порядка n определяется формулой

Вычислите обратную к Hⁿ матрицу для n = 1, 2, …, 20, 25, 30, 35, 40, 45, 50. Вы, несомненно, понимаете, что результат получится не вполне точным из-за небольших погрешностей машинной арифметики, но он должен быть очень близок к точной обратной матрице. Мерой погрешности служит левая остаточная матрица L = (Hⁿ)⁻¹Hⁿ − I и правая остаточная матрица R = Hⁿ(Hⁿ)⁻¹ − I; обе эти матрицы должны быть нулевыми, но, вероятно, не будут.

Конечно, если бы все элементы матриц L и R были порядка, скажем, 10⁻²⁰, то мы бы не имели забот. Для всех практических целей 10⁻²⁰ есть нуль, если элементы исходной матрицы равны в среднем 1/50 или больше. Существует, однако, точный способ оценки величины остаточных матриц L и R. Определим норму по строкам матрицы А как

Добавьте к своей программе, которая вычисляет обратную к матрице Гильберта, подпрограмму, печатающую таблицу |L|_r и |R|_r для каждой обратной матрицы. Проверьте вашу программу на отсутствие ошибок. Не сможете ли вы теперь объяснить, почему остаточные матрицы столь велики. Уверены ли вы в правильности программы?

Ваша программа правильна; причина неполадок — погрешность машинной арифметики. Матрицы Гильберта внешне выглядят вполне безобидно, однако они специально предназначены для демонстрации накопления ошибок в длинном ряду взаимосвязанных вычислений. Вы, быть может, считаете источником бед то, что ваш компьютер хранит недостаточное число цифр вещественных чисел. На многих ЭВМ имеется арифметика двойной точности. Предусмотрев в своем алгоритме двойную точность, вы сможете улучшить ситуацию, но заведомо не сможете полностью решить проблему. Весь этот этюд посвящен изучению влияния арифметики ограниченной точности на алгоритмы, которые являются абсолютно точными для «действительных» чисел (как их понимают математики). Прикладные математики и специалисты по численным методам в программистских лабораториях тратят большую часть времени на изменение теоретических алгоритмов, чтобы они могли работать на реальных ЭВМ[30].

Тема. Запрограммируйте алгоритм обращения матриц и проверьте его на матрицах Гильберта указанных выше порядков. Напечатайте таблицу или начертите график зависимости |L|_r и |R|_r от порядка n матрицы Hⁿ. Если используемый вами язык допускает выбор точности чисел, то повторите вычисление обратных матриц с большей точностью, чтобы увидеть, улучшится ли в результате таблица или график ошибок. (Мудрый программист так составит программу, чтобы изменение точности достигалось путем замены небольшого числа деклараций.) Следите также за числом фактических перестановок строк при выборе ведущего элемента; оно будет показывать, насколько плохо алгоритм согласуется с теорией.

30

Фактически поиск максимального в столбце элемента М на шаге 3 алгоритма обращения есть одно из таких изменений. М называется ведущим элементом, а сама операция — выбором ведущего элемента; на самом деле необходимо лишь, чтобы М был ненулевым. Максимальный элемент используется, чтобы уменьшить арифметическую погрешность ЭВМ. При обращении матрицы Гильберта ведущим элементом всегда должен оказываться ; если же алгоритм выбирает в качестве ведущего элемент, лежащий ниже, то это означает, что погрешность уже очень велика.