Мозес (Moses J.). Algebraic Simplification: A Guide for the Perplexed, CACM, 14, 8, pp. 527–537, 1971.
Мозес (Moses J.). Symbolic Integration: The Stormy Decade, CACM, 14, 8, pp. 548–560, 1971.
Этот выпуск САСМ целиком посвящен символьной алгебре и ее приложениям. Две статьи Мозеса — хорошие обзоры, но остальные статьи тоже интересны. Библиография в этих статьях должна помочь в исследовании любой темы во всей этой области.
21.
Превратное обратное,
или Ошибки при работе с плавающей точкой
Многие из методов, которые сейчас изучаются в средней школе, создавались величайшими математиками в течение столетий. Среди них — методы решения системы линейных уравнений, которые неявно включают методы обращения квадратных матриц. Начинающий алгебраист, изучая эти алгоритмы, может усомниться в том, что они всегда будут работать; но, испробовав метод на двух-трех примерах* наш скептик отбросит всякие сомнения. Он даже себе не представляет, какой его ждет удар: программа, написанная им в соответствии с простым и обоснованным алгоритмом, дает совершенно неверные результаты. Разве можно заподозрить, чтобы метод обращения матриц, придуманный королем математиков Гауссом, оказался несостоятельным?
Прежде всего освежим в памяти основные положения. Матрица — это квадратный массив вещественных чисел, в котором по горизонтали и вертикали располагается по n ≥ 1 элементов. Произведение С матрицы А справа на матрицу В записывается в виде С = АВ и задается формулой
Здесь подразумевается, что А, В и С — матрицы размера n × n. Умножение некоммутативно; можно найти такие матрицы А и В, что АВ ≠ ВА. Обратной матрицей к матрице А будет такая матрица А−1, что
AА−1 = А−1A = I
где I — единичная матрица, определяемая формулами Iii = 1 и Iij = 0 для i ≠ j. Большинство матриц имеет обратные, но не все. К сожалению, простейший способ обнаружить такие вырожденные матрицы состоит в том, чтобы попытаться вычислить обратную матрицу и потерпеть неудачу.
Как вычислить обратную матрицу? Следующий алгоритм принадлежит Гауссу.
Во-первых, положите матрицу X равной матрице I. В процессе вычислений матрица А будет в конце концов преобразована в I, матрица X, которая изначально была единичной матрицей, станет обратной к матрице А.
Для каждого столбца А, начиная со столбца 1 слева и кончая столбцом n справа, выполните следующее: Обозначим столбец, который будет обрабатываться на каждом этапе, символом j.
Пусть 
есть наибольший по абсолютной величине элемент в столбце j ниже строки j − 1. Если М равно нулю, то А — вырожденная матрица и продолжать обращение не имеет смысла. В противном случае поменяйте местами в обеих матрицах А и X строку j и строку, в которой находится М. И наконец, разделите каждый элемент в строке j матриц А и X на новое значение Ajj.
Теперь для всех строк i, i ≠ j, выполните все вычитания:
Aik = Aik − AijAjk, j ≤ k ≤ n,
Xik = Xik − AijXjk, 1 ≤ k ≤ n.
Результатом всех этих поэлементных вычитаний будет вычитание из строки i строки j с коэффициентом Ajj в обеих матрицах А и X. После выполнения этого шага для всех j все элементы сверху и снизу от Ajj станут нулевыми, а сам элемент Ajj будет равен единице. В матрице А не нужно выполнять вычитания слева от столбца j, поскольку все элементы строки j слева от Ajj равны нулю.
Для любого алгебраиста будет одно удовольствие доказать, что этот алгоритм всегда работает правильно и после его остановки X = А−1, если матрица А невырожденна. Вы едва ли найдете алгоритм, более приспособленный для структурной реализации. Почему бы нам, исключительно ради забавы, не провести небольшую проверку? Матрица Гильберта Hn порядка n определяется формулой
Вычислите обратную к Hn матрицу для n = 1, 2, …, 20, 25, 30, 35, 40, 45, 50. Вы, несомненно, понимаете, что результат получится не вполне точным из-за небольших погрешностей машинной арифметики, но он должен быть очень близок к точной обратной матрице. Мерой погрешности служит левая остаточная матрица L = (Hn)−1Hn − I и правая остаточная матрица R = Hn(Hn)−1 − I; обе эти матрицы должны быть нулевыми, но, вероятно, не будут.