Добавлен: 28.11.2018
Просмотров: 8365
Скачиваний: 76
Тема 19: Ошибки при работе с плавающей точкой
Многие из методов, которые сейчас изучаются в средней школе, создавались величайшими математиками в течение столетий. Среди них — методы решения системы линейных уравнений, которые неявно включают методы обращения квадратных матриц. Начинающий алгебраист, изучая эти алгоритмы, может усомниться в том, что они всегда будут работать; но, испробовав метод на двух-трех примерах* наш скептик отбросит всякие сомнения. Он даже себе не представляет, какой его ждет удар: программа, написанная им в соответствии с простым и обоснованным алгоритмом, дает совершенно неверные результаты. Разве можно заподозрить, чтобы метод обращения матриц, придуманный королем математиков Гауссом, оказался несостоятельным?
Прежде всего освежим в памяти основные положения. Матрица — это квадратный массив вещественных чисел, в котором по горизонтали и вертикали располагается по 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; обе эти матрицы должны быть нулевыми, но, вероятно, не будут.
Конечно, если бы все элементы матриц L и R были порядка, скажем, 10−20, то мы бы не имели забот. Для всех практических целей 10−20 есть нуль, если элементы исходной матрицы равны в среднем 1/50 или больше. Существует, однако, точный способ оценки величины остаточных матриц L и R. Определим норму по строкам матрицы А как
Добавьте к своей программе, которая вычисляет обратную к матрице Гильберта, подпрограмму, печатающую таблицу |L|r и |R|r для каждой обратной матрицы. Проверьте вашу программу на отсутствие ошибок. Не сможете ли вы теперь объяснить, почему остаточные матрицы столь велики. Уверены ли вы в правильности программы?
Ваша программа правильна; причина неполадок — погрешность машинной арифметики. Матрицы Гильберта внешне выглядят вполне безобидно, однако они специально предназначены для демонстрации накопления ошибок в длинном ряду взаимосвязанных вычислений. Вы, быть может, считаете источником бед то, что ваш компьютер хранит недостаточное число цифр вещественных чисел. На многих ЭВМ имеется арифметика двойной точности. Предусмотрев в своем алгоритме двойную точность, вы сможете улучшить ситуацию, но заведомо не сможете полностью решить проблему. Весь этот этюд посвящен изучению влияния арифметики ограниченной точности на алгоритмы, которые являются абсолютно точными для «действительных» чисел (как их понимают математики). Прикладные математики и специалисты по численным методам в программистских лабораториях тратят большую часть времени на изменение теоретических алгоритмов, чтобы они могли работать на реальных ЭВМ.
Тема. Запрограммируйте алгоритм обращения матриц и проверьте его на матрицах Гильберта указанных выше порядков. Напечатайте таблицу или начертите график зависимости |L|r и |R|r от порядка n матрицы Hn. Если используемый вами язык допускает выбор точности чисел, то повторите вычисление обратных матриц с большей точностью, чтобы увидеть, улучшится ли в результате таблица или график ошибок. (Мудрый программист так составит программу, чтобы изменение точности достигалось путем замены небольшого числа деклараций.) Следите также за числом фактических перестановок строк при выборе ведущего элемента; оно будет показывать, насколько плохо алгоритм согласуется с теорией.
Указания исполнителю. В этом этюде требуется прямая реализация алгоритма. Единственная трудность — это проверка соответствия программы теоретическим определениям и алгоритму. Не делайте над алгоритмом никаких оптимизирующих преобразований; вы изучаете, до чего можно дойти, если слепо следовать советам математиков, забыв о важнейшем положении — о точности вычислений.
Развитие темы. Если в достаточной степени расширить задачу, то она послужит основой семестрового курса методов вычислений. Тем не менее вы можете получить дополнительную информацию о поведении ошибок, если вычислите |L| и |R| с использованием других норм, отличных от нормы по строкам, например нормы по столбцам:
Ниже определены нормы L1, L2 и L∞:
Дополните значениями этих норм вашу таблицу анализа ошибок. Наблюдаются ли какие-либо существенные отличия одной нормы от остальных по форме или приблизительному положению кривой ошибок?
Тема 20: Арифметические вычисления с высокой точностью
Математику часто считают сухой наукой, однако и математику творили люди. Одной из самых печальных была судьба Уильямса Шенкса, жившего в девятнадцатом веке и посвятившего себя вычислению числа π с высокой точностью. Закончив многолетний труд, Шенкс в 1837 г. опубликовал значение π до 707-го десятичного знака, впоследствии исправив некоторые знаки. Может быть, надо счесть за благо, что Шенкс умер в 1882 г., поскольку в 1946 г. было показано, что его вычисления ошибочны начиная с 528-го десятичного знака. Фактически Шенкс не продвинулся дальше своих предшественников.
Полученное Шенксом значение было проверено, вероятно, с помощью механических устройств, а компьютер был впервые использован для вычисления π только в 1949 г.; это была машина ENIAC. Даже тогда проект был монументальным. Джорж У. Рейтуиснер писал: «Поскольку получить машину в рабочее время было практически невозможно, мы воспользовались разрешением выполнить эту работу за 4 выходных дня в период летних отпусков, когда ENIAC стоял без дела». Собственно вычисления (не программирование!) заняли 70 часов: было получено несколько больше 2 000 цифр. Все это время приходилось постоянно обслуживать компьютер, поскольку из-за ограниченности его возможностей требовались постоянная перфорация и ввод промежуточных результатов. Те первые программисты так же далеки от нынешних, как Шенкс далек от них.
Как бы мы стали вычислять π? Во-первых, необходимо выражение, которое можно вычислять. Ряд
π/4 = 1 − 1/3 + 1/5 − 1/7 + 1/9 − …
довольно прост для понимания, но он ужасно медленно сходится. Гораздо лучше ряд для арктангенса
arctg x = x − х3/3 + х5/5 − х7/7 + …, |x| ≤ 1
Объединим его с формулой сложения для тангенса
tg (a + b) = (tg a + tg b)/(1 − tg a · tg b)
и выберем а и b так, чтобы tg (a + b) = 1 = tg π/4. (Учитывая, что tg (arctg x) = х для −π/2 < х < π/2, можно взять, например, а = arctg (1/2), b = arctg (1/3).)
Тогда
arctg (tg (a + b)) = a + b = arctg 1 = π/4,
и теперь можно использовать приведенный выше ряд для нахождения a и b. На практике чаще всего используются следующие суммы:
π/4 = 4 arctg (1/5) − arctg (1/239),
π/4 = 8 arctg (1/10) − 4 arctg (1/515) − arctg (1/239),
π/4 = 3 arctg (1/4) + arctg (1/20) + arctg (1/1985).
Теперь мы собираемся просуммировать эти ряды на ЭВМ. Как известно, все, что нужно для суммирования, — это простой итерационный цикл, но тут возникает одна проблема. Точность вычислений на ЭВМ ограничена, а весь смысл этого упражнения в том, чтобы найти много-много цифр числа π, значительно превзойдя обычную точность. Первое, что приходит в голову, — промоделировать ручные методы выполнения арифметических действий. Будем представлять числа очень большими целочисленными массивами (по одной десятичной цифре в каждом элементе), тогда ясно, как составить программы сложения, вычитания и умножения. Запрограммировать ручной метод деления несколько сложнее, но все же возможно. Неприемлемым, однако, оказывается время выполнения алгоритмов. Хотя на это редко обращают внимание, но при ручных методах для умножения или деления n-значных чисел требуется время, пропорциональное n². Если речь идет об операциях над числами из тысяч цифр, то такие расходы будут нам не по карману. К счастью, имеются лучшие алгоритмы.