Численные методы

Министерство образования и науки Российской Федерации Южно-Уральский государственный университет Кафедра информационных технологий в экономике Б.М. Суховилов ЧИСЛЕННЫЕ МЕТОДЫ Конспект лекций Одобрено учебно-методической комиссией факультета экономики и управления Рецензенты: д.э.н., доцент С.А. Аристов к.т.н., доцент А.Г. Комирев Численные методы: конспект лекций / Б.М. Суховилов – Челябинск: ЮУрГУ. – 51 с. Материалы конспекта лекций соответствуют Федеральному государственному образовательному стандарту высшего профессионального образования третьего поколения. В первом разделе пособия рассмотрены вопросы машинного представления целых и действительных чисел, методы сравнения действительных чисел, способы корректной реализации финансовых расчетов, неустойчивые алгоритмы, чувствительные задачи и анализ сложности алгоритмов. Во втором разделе рассматривается математический пакет программ с открытыми кодами GNU Octave, предназначенный для моделирования решения инженерных задач в специализированной вычислительной среде, включающей встроенный язык программирования и набор профессионально реализованных алгоритмов в виде готовых стандартных функций. Каждый раздел снабжен примерами, иллюстрациями и контрольными упражнениями. Для бакалавров, обучающихся по направлению «Прикладная информатика». ВВЕДЕНИЕ «ЭВМ многократно увеличивает некомпетентность вычислителя» По мотивам принципа некомпетентности Питера. Решение инженерных задач часто нельзя выполнить в аналитической форме. К аналитическому решению надо стремиться, но это не всегда возможно в теоретическом плане или может оказаться нецелесообразным, например, из-за временных ограничений. В этом случае, задачу пытаются решить приближенно численно, применяя компьютер. Цель данного учебного материала – научить читателя решать вычислительные задачи эффективно. Под эффективностью будем понимать, прежде всего, максимизацию точности вычислений при минимизации сложности вычислительного алгоритма и затрат на его реализацию. Отметим, что затраты на реализацию можно сократить моделированием решения инженерной задачи в специализированной вычислительной среде, включающей несложный встроенный язык программирования, и большой набор профессионально реализованных алгоритмов в виде готовых стандартных функций. Указанное моделирование позволяет сосредоточиться на содержательной стороне поиска и исследования решения задач. Последние параграфы учебного пособия посвящены математическому пакету с открытыми кодами GNU Octave. GNU Octave с успехом используется для моделирования решений вычислительных задач. После выполненного моделирования, если требуется увеличить производительность вычислений, можно транслировать решение на любой подходящий язык программирования, например С или С++. Конспект лекций разработан в рамках ФГОС третьего поколения по направлению бакалавриата «Прикладная информатика» для дисциплин «Численные методы» и «Математические пакеты программ». В конспекте рассматриваются вопросы применения вычислительной математики и инструментальных средств реализации численных методов решения инженерных задач, которые отражают (по мнению автора) специфику подготовки бакалавров прикладной информатики, и пока не достаточно освещены и систематизированы в учебной литературе. Примеры, на которые есть ссылки в пособии, можно найти по адресу: http://inf.susu.ac.ru:8000/sukhovilov/na/samples.zip. Как правило, решение вычислительных задач осуществляется на множестве действительных чисел. Знания об их представлении и преобразованиях в компьютере чрезвычайно важно. Поэтому изложение начнем с описания машинного (компьютерного) представления действительных чисел. 3 1. МАШИННОЕ ПРЕДСТАВЛЕНИЕ ДЕЙСТВИТЕЛЬНЫХ ЧИСЕЛ Действительные числа при их использовании в расчетах на ЭВМ аппроксимируются машинным представлением. Наиболее распространен метод представления – числа с плавающей точкой. Данное представление обеспечивает компромисс между точностью и диапазоном представления действительных чисел. Число с плавающей точкой хранится в виде набора разрядов (позиционной системы счисления с некоторым основанием b), разделенных на знак, экспоненту (порядок) и мантиссу в следующем виде: Доказано, что числа с плавающей точкой при b = 2 (двоичная система счисления) достаточно устойчивы к ошибкам округления. Двоичная система счисления удобна для аппаратной реализации. Поэтому в дальнейшем изложении полагаем b = 2, и формулу числа с плавающей точкой запишем в виде: s (-1) × m × bE , где s – знак; b = 2 – основание; E – экспонента (порядок); m – мантисса. Множество чисел с плавающей точкой описывается параметрами:  |m| – количество цифр мантиссы, определяющее точность представления чисел;  Emin, Emax – нижняя и верхняя границы экспонент, определяющие диапазон представления чисел. Множество чисел с плавающей точкой обладает определенными свойствами. Числа по модулю большие максимального числа из множества и меньшие наименьшего числа из множества не могут быть отображены числом с плавающей точкой. Множество чисел с плавающей точкой не является непрерывным, бесконечным множеством подобно множеству действительных чисел, образующих континуумом. Отсюда следует вывод, что множество чисел с плавающей точкой лишь приближенно представляет множество действительных чисел. Рассмотрим пример с b = 2, |m| = 3, Emin = +2, Emax = +5. Число 4 можно представить в следующем виде: 410 = 1002 = 1.00 × 2+2 = 0.100 × 2+3 = 0.010 × 2+4 = 0.001 × 2+5. Из примера видно, что одно и то же число можно представить по-разному. Чтобы избежать этого, числа с плавающей точкой представляют в нормализованном виде, когда 1-й бит мантиссы всегда равен 1: 4 Нормализованное число имеет следующий вид: s (-1) × 1.M × 2E, где M – это остаток мантиссы, |M| = n. Так как неявную единицу не нужно хранить в памяти, такое представление экономит один бит и обеспечивает уникальность представления числа. В нашем примере «4» имеет единственное представление 1.00 × 2+2, а мантисса хранится в памяти как «00», т.к. старшая единица подразумевается неявно. Нормализованное представление чисел с плавающей точкой имеет свои проблемы, например, невозможность в данной форме представить 0. Также имеет место неравномерность заполнения действительной оси множеством нормализованных чисел с плавающей точкой, например, наличие «околонулевой ямы». Все это требует дальнейших уточнений и дополнений, отраженных в соответствующих стандартах. Упражнение 1. Доказать, что в нормализованном множестве чисел с плавающей точкой ровно 2× (b – 1) × bn (Emax – Emin + 1) + 1 отрицательных и положительных чисел вместе с 0. В данном случае n – длина остатка мантиссы. Упражнение 2. Показать, что нормализованное множество чисел с плавающей точкой расположено неравномерно на действительной оси чисел. 1.1. Стандарт представления чисел с плавающей точкой В компьютерных системах для вычислений с плавающей точкой широко используется стандарт IEEE 754, разработанный ассоциацией Institute of Electrical and Electronics Engineers (IEEE). Стандарт IEEE 754 определяет:  представление нормализованных чисел с плавающей точкой;  представление положительного и отрицательного нулевого числа с плавающей точкой (+ 0, – 0);  представление положительной и отрицательной бесконечности (+ ∞, – ∞);  представление специальной величины "Не число (Not-a-Number)" (NaN);  представление денормализованных чисел с плавающей точкой («решение» проблемы «околонулевой ямы»);  режимы округления чисел с плавающей точкой. IEEE 754 также определяет четыре формата представления чисел с плавающей точкой:   с одинарной точностью (single-precision) 32 бита; с двойной точностью (double-precision) 64 бита; 5  с одинарной расширенной точностью (single-extended precision) >=43 бит (редко используемый);  с двойной расширенной точностью (double-extended precision) >= 79 бит (обычно используют 80 бит). Ниже представлен формат хранения чисел с плавающей точкой в стандарте IEEE 754: При S = 0 – положительное число, S = 1 – отрицательное число; E – смещенная экспонента двоичного числа; M – остаток мантиссы двоичного числа (с учетом того, что для нормализованных чисел есть еще не хранимая 1). Формула вычисления десятичных чисел с плавающей точкой, из чисел, представленных в стандарте IEEE 754: ( m 1) 1)] M F  ( 1) S 2[ E  ( 2 (1  n ) , 2 ( m 1)  1) – смещение экспоненты (например, в 32-битном числе при m = 8 где ( 2 оно равно +127). Смещение позволяет определять как отрицательные, так и положительные экспоненты, без выделения специального бита, хранящего знак. На рис. 1, 2 даны представления чисел стандарта IEEE 754 одинарной (float) и двойной (double) точности и формулы вычисления соответствующих десятичных чисел: Рис. 1. Формат хранения float – числа одинарной точности (single-precision) 32 бита Рис. 2. Формат хранения double – числа двойной точности (double-precision) 64 бита Рассмотрим использование стандарта IEEE 754 на примере среды разработки Microsoft Visual Studio1. В Microsoft Visual C++ существует три внутренних типов хранения действительных чисел в формате представления чисел с плавающей точкой из стандарта IEEE 754: действительные числа размером в 4 байта (тип данных float), 8 байт (тип данных double) и 10 байт 1 Описание взято из MSDN Microsoft Visual Studio 2012. 6 (тип данных double-extended поддерживается на уровне ассемблера). Значения хранятся в следующем виде (табл. 1). Таблица 1 Значение Хранится как float разряд знака, 8-разрядная экспонента, 23-разрядная мантисса double разряд знака, 11-разрядная экспонента, 52-разрядная мантисса double-extended разряд знака, 15-разрядная экспонента, 64-разрядная мантисса Степени являются степенями числа 2 и смещены на половину своего возможного значения. Это означает, что для того, чтобы получить актуальное значение степени, следует вычесть данное смещение из хранимой степени. Если хранимая степень меньше, чем смещение, то это отрицательная степень. Экспоненты типов float и double определяются следующим образом (табл. 2). Таблица 2 Format Экспонента Смещение Emin Emax float 8-разрядная 127 –126 127 double 11-разрядная 1023 –1022 1023 При E = 1, Emin = 1 – 127 = -126. При E = 254, Emax = 254 – 127 = 127. E = 0 и E = 255 зарезервированные значения для представления специальных случаев хранения нуля и бесконечности. В форматах хранения типов float и double вначале мантиссы предполагается (для нормализованных чисел) наличие единицы «1», которая не хранится в памяти, поэтому мантисса на самом деле 24- или 53-разрядная, хотя хранятся только 23 или 52 разряда. В формате double-extended этот разряд хранится. Таким образом, мантисса хранится в виде двоичной дроби формы 1.ХХХ... Х (для нормализованных чисел). Значение мантиссы больше или равно 1 и меньше или равно 2. В таком случае для различных размеров используется следующий формат (табл. 3). Таблица 3 Format BYTE 1 BYTE 2 BYTE 3 BYTE 4 float SXXX XXXX XMMM MMMM MMMM MMMM MMMM MMMM double SXXX XXXX XXXX MMMM MMMM MMMM MMMM MMMM 7 .. . BYTE n .. MMMM . MMMM doubleextended SXXX XXXX XXXX XXXX 1MMM MMMM MMMM MMMM .. MMMM . MMMM где S – знаковый разряд, X – разряды степени, а M – разряды остатка мантиссы. Далее приведены несколько примеров представления чисел типа float. Представление числа 2. Знаковый разряд равен нулю; смещенная степень имеет значение 12810 = (127 + 1)10 = 100000002; мантисса равна (1.) 0000000...0000 00002 = 110. 210 = (1 * 21) 10 = 0100 0000 0000 0000 ... 0000 00002 = 4000 000016. Представление числа -2. -210 = (-1 * 21) 10= 1100 0000 0000 0000 ... 0000 00002 = C000 000016. Представление числа 4. Та же мантисса, степень увеличена на единицу (смещенное значение 12910 = 100000012). 410 = (1 * 22) 10= 0100 0000 1000 0000 ... 0000 00002 = 4080 000016. Представление числа 6. Та же степень, мантисса больше на половину и равна (1.) 100 0000 ... 0000 00002 = 1.510. 610 = (1.5 * 22) 10= 0100 0000 1100 0000 ... 0000 00002 = 40C0 000016. Представление числа 1. Степень равна 12710 = 011111112. Мантисса, как у всех чисел являющихся степенями 2. 110 = (1 * 20) 10= 0011 1111 1000 0000 ... 0000 00002 = 3F80 000016. Представление числа 0.75. Смещенная степень имеет значение 12610 = 011 1111 02. Мантисса имеет значение (1.) 100 0000 ...0000 00002 = 1.510. 0.7510 = (1.5 * 2-1) 10= 0011 1111 0100 0000 ... 0000 00002 = 3F40 000016. Представление числа 2.5.То же представление, что и для 2, за исключением того, что в мантиссе задан разряд, представляющий 0.25. 2.510 = (1.25 * 21)10= 0100 0000 0010 0000 ... 0000 00002 = 4020 000016. Представление числа 0.1. Одна десятая (0.1) представляет собой периодическую дробь в двоичной форме. Причина того, почему 1/10 и 1/100 не представимы точно в двоичной форме, сходна с причиной того, почему 1/3 не представима точно в десятичной форме. Мантисса представима только как 1.6 [(1).+0.6]. Для получения 0.1 требуется деление 1,6 на 16 (степень равна 12310 = 011110112, т. к. 123 – 127 = –4). Следует учесть, что хранимая мантисса округляется в большую сторону в последнем разряде, чтобы таким образом представить непредставимое число с максимально возможной точностью: 0.110 = (1.6 * 2-4) 10 = 1.100110011001100110011012 × 2-4 = 0011 1101 1100 1100 ... 1100 11012 = 3DCC CCCD16. Примечание. Перевод дробного числа из десятичной системы счисления в двоичную осуществляется по следующему алгоритму:     вначале переводится целая часть десятичной дроби в двоичную систему счисления; затем дробная часть десятичной дроби умножается на основание двоичной системы счисления; в полученном произведении выделяется целая часть, которая принимается в качестве значения первого после запятой разряда числа в двоичной системе счисления; алгоритм завершается, если дробная часть полученного произведения равна нулю или если достигнута требуемая точность вычислений. В противном случае вычисления 8 продолжаются с предыдущего шага при отбрасывании целой части полученного произведения. Таким образом, перевод 0.610 в двоичную систему счисления представляет периодическую дробь: 0.6*2=1.2; 0.2*2=0.4; 0.4*2=0.8; 0.8*2=1.6; 0.6*2=1.2; … Упражнение 3. Показать, что формы представления чисел с плавающей n точкой: ( 1) s (1   bn i 2i )2( E offset ) и ( 1) s (1  M  2  n )2( E  offset ) эквивалентны. i 1 Здесь b = 0/1 – цифры мантиссы, М – остаток мантиссы, |M |= n. 1.2. Специальные случаи представления чисел с плавающей точкой Рассмотрим представление нуля и бесконечности [1]. Поскольку в нормализованном представлении чисел с плавающей точкой нельзя получить ноль, в стандарте IEEE 754 число «±0» определяется специальным значением с экспонентой, равной E = Emin – 1 (для float это – 127) и нулевой мантиссой. Также в IEEE 754 предусмотрено представление для специальных чисел, работа с которыми вызывает исключение. К таким числам относится бесконечность (± ∞) и неопределенность, «не число» (NaN). Эти числа позволяют вернуть адекватное значение при переполнении. Бесконечности представлены как числа с порядком E = Emax + 1 и нулевой мантиссой. Получить бесконечность можно при переполнении и при делении ненулевого числа на ноль. Бесконечность при делении разработчики определили исходя из существования пределов, когда делимое и делитель стремятся к какому-то числу. Соответственно, c/0 = ± ∞ (например, 3/0 = + ∞, а –3/0 = – ∞). При 0/0 предел не существует, поэтому результатом будет неопределенность. Неопределенность или NaN – это представление, придуманное для того, чтобы арифметическая операция могла всегда вернуть какое-то не бессмысленное значение. В IEEE 754 NaN представлен как число, в котором E = Emax + 1, а мантисса не нулевая. Любая операция с NaN возвращает NaN. При желании в мантиссу можно записывать информацию, которую программа сможет интерпретировать. Стандартом это не оговорено и мантисса чаще всего игнорируется. NaN можно получить одним из следующих способов: ∞+(–∞); 0×∞; 0/0, ∞/∞; sqrt(x), где x < 0. По определению NaN ≠ NaN, поэтому, для проверки значения переменной нужно просто сравнить ее с собой. Таким образом, в описанном представлении существуют «+0» и «–0»: «+∞» и «–∞»: 9 NaN, «не число» (цифра в битах с 0…22 не может быть 0, т.е. «+∞», «–∞»): В стандарте знак сохранили умышленно, чтобы выражения, которые в результате переполнения или потери значимости превращаются в бесконечность или в ноль, при умножении и делении все же могли представить максимально корректный результат. Например, если бы у нуля не было знака, выражение 1/(1/x) = x не выполнялось бы верно при x = ±∞, так как 1/∞ и 1/–∞ были бы равны 0. Еще один пример: (+∞/0) + ∞ = +∞, тогда как (+∞/– 0) +∞ = NaN. Чем бесконечность в данном случае лучше, чем NaN? Тем, что если в арифметическом выражении появился NaN, результатом всего выражения всегда будет NaN. Если же в выражении встретилась бесконечность, то результатом может быть ноль, бесконечность или обычное число с плавающей точкой. Например, 1/∞ = 0. Иллюстрирующий пример к разделу – проект FloatPoint. В листинге FloatPoint.cpp описаны рассмотренные выше специальные случаи, вычисляется порядок следования байтов, описывается и моделируется очередность чисел в IEEE 754, приводится вычисление машинного эпсилон, характеризующего относительную точность плавающей арифметики, как наименьшего числа с плавающей точкой ε, такого, что 1 + ε > 1. 1.3. Представление денормализованных чисел Что такое денормализованные (subnormal) числа рассмотрим на примере [1]. Пусть имеем нормализованное представление с длиной остатка мантиссы |M| = 2 бита (+ один бит нормализации) и диапазоном значений порядка –1 ≤ E ≤ 2. В этом случае получим 16 чисел: Крупными штрихами показаны числа с мантиссой, равной 1,00. Видно, что расстояние от нуля до ближайшего числа (0 – 0,5) больше, чем от этого числа к следующему (0,5 – 0,625). Это значит, что разница двух любых чисел от 0,5 до 1 даст 0, даже если эти числа не равны. Что еще хуже, в пропасть между 0 и 0,5 попадает разница чисел, больших 1. Например, «1,5 – 1,25 = 0». Разработчики стандарта предложили добавить числа с минимально возможным порядком E = Emin, у которых неявный старший бит мантиссы равен 0. Такие числа называются денормализованными. Они накрывают диапазон между 0 и ближайшими нормализованными числами справа и слева от 0. Как это работает? Если при попытке нормализации числа получен порядок E = Emin – 1, а мантисса не нулевая, то число считается денормализованным, его порядок полагается E = Emin, а неявный старший бит 10 мантиссы полагается равным нулю. Таким образом, числа с плавающей точкой теперь имеют вид: s (-1) × 1.M × 2E, если Emin ≤ E ≤ Emax (нормализованные числа); s (-1) × 0.M × 2Emin, если E = Emin – 1 (денормализованные числа) . (1) (2) Для типа данных float формулы (1), (2) запишутся в виде: s (-1) × (1+M/223) × 2E-127 (нормализованные числа) ; s s (-1) × M/223 × 2-126 = (-1) × M × 2-149 (денормализованные числа) . (1а) (2а) Вернемся к примеру. Интервал от 0 до 0,5 теперь заполняют денормализованные числа: 0.01×2-1, 0.10×2-1, 0.11×2-1. В результате получаем новое представление чисел: Теперь, например, (1,5 – 1,25) = 1.1×20 – 1.01×20 = 0.1×2-1  0, что дает возможность не проваливаться в 0, как в рассмотренных выше примерах (0,5 – 0,25 и 1,5 – 1,25). Это сделало представление более устойчивым к ошибкам округления для чисел, близких к нулю. Программа FloatPoint1.m (среда Matlab) демонстрирует расположение на действительной оси нормализованных и денормализованных чисел, накрывающих «околонулевую» яму. Проект Subnormal и скрипт subnormal.m для среды Matlab демонстрируют работу с денормализованными числами. В проекте na/subnormal.c, на примере операционной системы Freebsd, представлена возможность компилятора gcc управлять представлением чисел с плавающей точкой, в частности, исключить денормализованные числа из расчетов путем присвоения им нуля. 1.4. Границы множества чисел с плавающей точкой На примере чисел одинарной точности (тип float) вычислим границы представления действительных чисел. Для этого подставляем значения минимальных и максимальных мантисс и порядков в формулы (1а) и (2а):  минимальное нормализованное число (абсолютное значение) 00 80 00 00 = 2-126∙(1+0/223)= 2-126 ≈ 1,17549435∙e-38 80 80 00 00 = -2-126∙(1+0/223)=-2-126 ≈ -1,17549435∙e-38 11  максимальное денормализованое число (абсолютное значение) 00 7F FF FF = 2-126∙(1-2-23) ≈ 1,17549421∙e-38 80 7F FF FF = -2-126∙(1-2-23) ≈ -1,17549421∙e-38  минимальное денормализованное число (абсолютное значение) 00 00 00 01 = 2-126∙ (0+2-23) = 2-149 ≈ 1,40129846∙e-45 80 00 00 01 = -2-126∙ (0+2-23) = 2-149 ≈ -1,40129846∙e-45  максимальное нормализированное число (абсолютное значение) 7F 7F FF FF = 2127∙(1+1-2-23) ≈ 3,40282347∙e+38 FF 7F FF FF = -2127∙(1+1-2-23) ≈ -3,40282347∙e+38 Из расчетов следует, что минимальное денормализованное число граничит с нулем, а максимальное денормализированное число граничит с минимальным нормализированным числом. Таким образом, использование денормализованных чисел решает проблему «околонулевой ямы». Как и следовало ожидать, максимальное нормализированное число граничит с бесконечностью. На рис. 3, 4 представлен весь диапазон чисел с плавающей точкой типов данных float и double [3]. Рис. 3. Диапазон чисел формата float (32 бита) по стандарту IEEE 754 12 Рис. 4. Диапазон чисел формата double (64 бита) по стандарту IEEE 754 1.5. Округление чисел с плавающей точкой Стандарт IEEE 754 требует, чтобы результат любой арифметической операции совпадал с результатом операции, выполненной над точными значениями, и округленным до ближайшего представимого в этом формате числа. Математически доказано [5], что если 0.5 округлять до 1 (в большую сторону), то существует набор операций, при которых ошибка округления будет возрастать до бесконечности. Поэтому, если округляемое число находится ровно посередине между меньшим и большим числом с плавающей точкой (возможными претендентами), в стандарте IEEE 754 применяется правило округления до четного. Т.е. последний бит округленного числа в этом случае всегда устанавливается в 0. Так, 12,5 будет округлено до 12, а 13,5 – до 14. 1.6. Абсолютная и относительная погрешность округления Обозначим точное значение числа через x , а его приближенное значение через x̂ . Тогда абсолютная (предельная абсолютная) погрешность есть: x̂  x   , где в качестве  стремятся указать наименьшее число, удовлетворяющее неравенству. Абсолютная погрешность подходит для сравнения денормализованных чисел, поскольку порядок этих чисел – Emin одинаков. 13 Для нормализованных чисел (как будет показано ниже) подходящей мерой близости чисел является относительная (предельная относительная) погрешность  : xˆ  x      x(1   )  xˆ  x(1   ) . x x Упражнение 4. Докажите следующие утверждения. 1. Абсолютная погрешность суммы и разности приближённых чисел равна сумме абсолютных погрешностей операндов. 2. Если погрешности приближённых чисел малы, то относительная погрешность их произведения приближённо (с точностью до членов более высокого порядка малости) равна сумме относительных погрешностей сомножителей. Подсказка. Используйте разложение функции произведения сомножителей в ряд Тейлора с точностью до членов первого порядка. 3. Если погрешности приближённых чисел малы, то относительная погрешность частного двух чисел приближённо (с точностью до членов более высокого порядка малости) равна сумме относительных погрешностей этих чисел. 1.7. Точность представления действительных чисел Функция погрешности представления числа в IEEE 754 показана на рис. 5: Рис. 5. Функция погрешности представления числа по стандарту IEEE 754 Шаг чисел удваивается с увеличением экспоненты двоичного числа на единицу. То есть, чем дальше от нуля, тем шире шаг чисел в формате IEEE 754 по числовой оси. Из рисунка видно, что абсолютная максимальная погрешность для числа в формате IEEE 754 равна половине шага чисел и увеличивается с ростом экспоненты числа с плавающей точкой. В листинге FloatPoint.cpp описан расчет абсолютной ошибки числа типа float вблизи 1. Программа err_f.m (пакет Matlab) использует файл данных err_.dat, полученный из расчетов проекта FloatPoint, для построения графика абсолютной ошибки числа типа float вблизи 1. Шаг нормализованного числа типа float в соответствии с (1а) равен: (1+1/223) × 2E-127 – (1+0/223) × 2E-127 = 2E-150 , а шаг денормализованного числа типа float в соответствии с (2а) равен: (0+1/223)× 2-126 – (0+0/223)× 2-126 = 1 × 2-149 – 0 × 2-149 = 2-149. 14 Соответственно абсолютная погрешность составит (тип float):  для нормализованного числа (½)2E-150 = 2E-151;  для денормализованного числа (½)2-149 = 2-150. Относительные погрешности нормализованных и денормализованных чисел составят (тип float): -n  2E-151/(1+ M /223) ×2E-127 = 1/(224+2× M ) ≤ (0.5×2 = 5.960464e-8, при M = 0, n = 23),  2-150/(0+ M /223) ×2-126 = 1/(2×M ) ≤ (0.5, при M = 1), при длине остатка мантиссы n = 23. Шаг нормализованного и денормализованного числа типа double равен: (1+1/252) × 2E-1023 – (1+0/252) × 2E-1023 = 2E-1075 , 1 × 2-1074 – 0 × 2-1074 = 2-1074. Соответственно абсолютная погрешность составит (тип double):  для нормализованного числа (½)2E-1075 = 2E-1076;  для денормализованного числа (½)2-1074 = 2-1075. Относительные погрешности нормализованных и денормализованных чисел составят (тип double): -n  2E-1076/(1+ M /252) × 2E-1023 = 1/(253+2× M ) ≤ (0.5×2 = 1.110223e-16, при M = 0, n = 52),  2-1075/(0+ M /252) × 2-1022 = 1/(2× M ) ≤ (0.5, при M = 1), при длине остатка мантиссы n = 52. Т. о. максимальное значение относительной погрешности представления нормализованных чисел с плавающей точкой зависит только от длины мантиссы, а абсолютная погрешность представления денормализованных чисел с плавающей точкой постоянна. В табл. 4, 5 приведены расчеты абсолютных и относительных погрешностей представления (ошибок округления) конкретных значений, охватывающих весь возможный диапазон, чисел (1-й столбец таблицы) типа float и double. Из расчетов следует, что денормализованные числа имеют фиксированную абсолютную погрешность и изменяющуюся относительную, а нормализованные числа наоборот имеют почти фиксированную относительную погрешность и изменяющуюся абсолютную. 15 Таблица 4 Погрешности представления чисел типа float 0000 0001 Абсолютная погрешность Число IEEE 754, hex 2-149 ≈1,401298e-45 Относительная погрешность Примечание 7.006492e-046 5.000000e-001 Subnormal, min 7.006492e-046 2.500000e-001 Subnormal 0000 0002 2 0000 0032 ≈7,00649e-44 7.006492e-046 1.000000e-002 Subnormal 007F FFFF ≈1,175494e-38 7.006492e-046 5.960465e-008 Subnormal, max 0080 0000 ≈1,175494e-38 7.006492e-046 5.960464e-008 Normal, FLT_MIN 0DA2 4260 ≈1,0e-30 4.701977e-038 4.701977e-008 Normal 1E3C E508 ≈1,0e-20 4.038968e-028 4.038968e-008 Normal 2EDB E6FF ≈1,0e-10 3.469447e-018 3.469447e-008 Normal 3F80 0000 =1,0 5.960464e-008 5.960464e-008 Normal, 4120 0000 ≈10,0 4.768372e-007 4.768372e-008 Normal 42C8 0000 ≈1,0e+2 3.814697e-006 3.814697e-008 Normal 5015 02F9 ≈1,0e+10 5.120000e+002 5.120000e-008 Normal 60AD 78EC ≈1,0e+20 4.398047e+012 4.398046e-008 Normal 7149 F2CA ≈1,0e+30 3.777893e+022 3.777893e-008 Normal 7F7F FFFF ≈+3,40282e+38 1.014120e+031 2.980232e-008 Normal, FLT_MAX -148 ≈2,802597e-45 FLT_EPSILON/2 Таблица 5 Погрешности представления чисел типа double Число IEEE 754, hex Абсолютная погрешность Относительная погрешность Примечание 00000000 00000001 2-1074 ≈4,940656e- 2.470328e-324 324 5.000000e-001 Subnormal, min 00000000 00000002 2-1073 ≈9,881313e- 2.470328e-324 324 2.500000e-001 Subnormal 00000000 00000032 ≈2,470328e-322 2.470328e-324 1.000000e-002 Subnormal 000FFFFF FFFFFFFF ≈2,225073e-308 2.470328e-324 1.110223e-016 Subnormal, max 00100000 00000000 ≈2,225074e-308 2.470328e-324 1.110223e-016 Normal, DBL_MIN 2B2BFF2E E48E0530 ≈1,0e-100 6.344855e-117 6.344855e-017 Normal 3FF00000 00000000 =1,0 1.110223e-016 1.110223e-016 Normal, 54B249AD 2594C37D ≈1,0e+100 9.713344e+083 9.713344e-017 Normal 6974E718 D7D7625A ≈1,0e+200 8.498208e+183 8.498208e-017 Normal 7FEFFFFF FFFFFFFF ≈1,79769e+308 9.979202e+291 5.551115e-017 Normal, DBL_MAX DBL_EPSILON/2 16 Из расчетов относительных погрешностей следует, что тип float позволяет (в среднем) правильно представить нормализованные действительные числа с 7-8 значащими цифрами, а тип double с 15-16. Упражнение 5. Написать программу для расчета табл. 4,5. Использовать multiprecision арифметику библиотеки Boost для вычисления ошибок представления чисел double в табл. 5. 1.8. Машинный эпсилон Точность плавающей арифметики удобно характеризовать весом младшего знака мантиссы вблизи числа 1, т. е. наименьшего числа с плавающей точкой ε, -n такого, что 1 + ε > 1. Значение ε составляет ε = 2 , где n – длина остатка мантиссы. Машинный эпсилон используется при теоретическом анализе ошибок округления алгоритмов, т. к. он имеет порядок значения относительной погрешности округления нормализованных чисел. На практике, машинный эпсилон применяется при назначении относительных порогов при сравнении нормализованных чисел с плавающей точкой (см. значения, присвоенные переменной maxRelDiff, в проекте Compare). В файле float.h определены константы машинных эпсилон для базовых типов данных float и double: соответственно FLT_EPSILON и DBL_EPSILON. Ниже приведен листинг подпрограммы из проекта FloatPoint, вычисляющей ε непосредственно как для базовых типов данных, так и для multiprecision типов cpp_dec_float_100 и cpp_dec_float_50 из библиотеки Boost: // Вычисление машинного эпсилон для различных типов данных с плавающей точкой template< typename T > T Epsilon (void) { T one = (T)1.; T two = (T)2.; T epsilon = one; while ((T)(one + epsilon) > one) epsilon = epsilon/two; return epsilon = two * epsilon;} Упражнение 6. Используя значения машинных эпсилон для типов данных cpp_dec_float_50 и cpp_dec_float_100, полученные «прямыми» вычислениями (в проекте FloatPoint), определить длины мантисс соответствующих типов. 1.9. Сравнение чисел с плавающей точкой Необходимо быть внимательным при сравнении чисел с плавающей точкой. Рассмотрим фрагмент программы из проекта Compare: float f = 0.1; double d = 0.1; 17 cout << cout << (f == 0.1) << endl; (d == 0.1) << endl; // print 0 // print 1 При первом и втором сравнении будут выведены разные значения: 0 и 1. Дело в том, что десятичное число 0.1, во-первых, не имеет точного представления в двоичном коде а, во-вторых, по умолчанию, 0.1 – это константа двойной точности (0.1d). При присваивании переменной одинарной точности f = 0.1 происходит отсечение части мантиссы (о чем может предупредить компилятор). Для выбора правильной процедуры сравнения чисел с плавающей точкой обратимся к табл. 4, 5. Как отмечалось ранее, денормализованные числа имеют фиксированную абсолютную погрешность и изменяющуюся относительную, а нормализованные числа наоборот имеют почти фиксированную относительную погрешность и изменяющуюся абсолютную. Поэтому процедура сравнения денормализованных чисел должна опираться на значение абсолютной погрешности, а нормализованных чисел на значение относительной погрешности. На основании сделанных заключений разработаны функции сравнения чисел с плавающей точкой (проект Compare): // Функция сравнения с абсолютной погрешностью maxAbsDiff // Применяется для денормализованных чисел, у которых абсолютная погрешность ~ постоянна bool IsFloatEqualAbs(float a, float b, float maxAbsDiff) { return fabs(a - b) <= maxAbsDiff; } /* Функция сравнения с относительной погрешностью maxRelDiff Применяется для нормализованных чисел, у которых относительная погрешность ~ постоянна Переменная maxRelDiff=iNum*FLT_EPSILON, где iNum определяет количество двоичных разрядов (начиная с младшего) в сравниваемых числах, которые разрешается упустить. Например, iNum=16, означает, что младшие 4 бита = log2(16) (это примерно 1,5 последней значащей десятичной цифры) могут не совпадать, а числа все равно будут считаться равными */ bool IsFloatEqualRelative(float a, float b, float maxRelDiff) { return fabs(a - b) <= fmax(fabs(a),fabs(b))*maxRelDiff; } /* Функция сравнения по maxUlps (от Units-In-Last-Place) – это максимальное количество чисел с плавающей точкой, которое может лежать между проверяемым и ожидаемым значением. Другой смысл этой переменной – это количество двоичных разрядов (начиная с младшего) в сравниваемых числах, которые разрешается упустить. Например, maxUlps=16, означает, что младшие 4 бита = log2(16) могут не совпадать, а числа все равно будут считаться равными. http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ */ union Float_t { Float_t(float num = 0.0f) : f(num) {} unsigned char sign() const { return (i >> 31); } unsigned __int32 i; float f; }; bool IsFloatEqualUlps(float a, float b, int maxUlpsDiff) { Float_t uA(a), uB(b); // Различные знаки -> неравные числа if (uA.sign() != uB.sign()) return false; int delta = uA.i - uB.i; 18 return abs(delta) <= maxUlpsDiff; } // Универсальная функция сравнения нормализованных и денормализованных чисел типа float bool IsFloatEqual(float a, float b, float maxAbsDiff, float maxRelDiff, int maxUlpsDiff = -1) { // Сравниваем по Ulps if (maxUlpsDiff > -1) return IsFloatEqualUlps(a, b, maxUlpsDiff); // Сравниваем по относительной погрешности // Сравнение денормализованных чисел, а также нормализованных из диапазона // c минимальной экспонентой if (fabs(b - a) < FLT_MIN) return IsFloatEqualAbs(a, b, maxAbsDiff); // Сравнение нормализованных чисел return IsFloatEqualRelative(a, b, maxRelDiff); } Отметим, что функция сравнения IsFloatEqualUlps [2] использует особенность представления чисел в формате IEEE 754, состоящую в том, что порядок и мантисса расположены друг за другом таким образом, что вместе образуют последовательность целых чисел {n}, для которых выполняется: n < n+1 ⇒ F(n) < F(n+1), где F(n) – число с плавающей точкой, образованное от целого n, разбиением его битов на порядок и мантиссу. Поэтому если взять положительное число с плавающей точкой, преобразовать его к целому, прибавить «1», мы получим следующее число, которое представимо в этой арифметике [1]. Универсальная функция IsFloatEqual, основываясь на значениях сравниваемых чисел, позволяет правильно выбрать метод сравнения. Указанные функции позволили правильно выполнить тестовые расчеты методом половинного деления корней уравнений, значения которых (для проверки) назначались во всем диапазоне абсолютных значений чисел типа float, как денормализованных, так и нормализованных (проект Compare). 1.10. Сравнение чисел с плавающей точкой с нулем Рассмотрим фрагмент кода из проекта CompareWithZero: float summa = 0; float step = 0.1f; for (int i = 1; i <= 10; i++) summa = summa + step; float delta = abs(1.f - summa); // Абсолютная ошибка представления денормализованных чисел float maxAbsDiff = pow(2.,-149); float maxRelDiff = 16*FLT_EPSILON; int maxUlpsDiff = 16; // 1) Вывод - 1 cout << IsFloatEqual(summa,1.f, maxAbsDiff, maxRelDiff, maxUlpsDiff) << endl; // 2) Вывод - 0 cout << IsFloatEqual(delta,0.f, maxAbsDiff, maxRelDiff, -1) << endl; // 3) Вывод - 0 cout << IsFloatEqual(delta,0.f, maxAbsDiff, maxRelDiff, maxUlpsDiff) << endl; 19 Сложив переменную step = 0.1f десять раз и получив результат в переменной summa, мы ожидаем, что значение переменной summa близко к 1. Так и есть – значение summa отличается от 1 на величину ошибки округления вблизи 1, т. е. на FLT_EPSILON или Ulps = 1. Это подтверждает результат 1-го вызова функции IsFloatEqual. Далее мы можем предположить, что очень близким к 0 должно быть значение переменной delta = abs(1.f - summa). Но это не так. Это подтверждает результат 2-го вызова функции IsFloatEqual, когда выполняется функция сравнения с относительной погрешностью IsFloatEqualRelative, и результат 3-го вызова функции IsFloatEqual, когда выполняется функция сравнения IsFloatEqualUlps. Оказывается, между значением delta = FLT_EPSILON и нулем расположены 872 415 232 действительных чисел, как показывает соответствующий вывод на печать из проекта CompareWithZero. Попробуем обосновать полученные результаты. Рассмотрим значения переменных: summa  1   , где  — погрешность округления результата вычислений, имеющая порядок машинного эпсилон чисел типа float и delta  summa  1   . Сравним delta с 0. Сравнивая числа "a =  " и "b = 0" с относительной погрешностью k   , где k — количество чисел с плавающей точкой между "a" и "b" имеем | a  b |  k    max(| a |, | b |) |   0 |  k   2  k  1   8  10 6 — не удовлетворительно. Сравним переменную summa с 1: | summa  1 |  k    max( summa , 1) |  |  k    k  1 — удовлетворительно. Приведенные рассуждения показывают сложности с принятием решения о порядках величин при их сравнении с 0. Поскольку, часто при сравнении некоторой величины с 0 сама величина получается в результате операции вычитания двух величин (или может быть приведена, в конечном счете, к результату такого вычитания), то можно посоветовать (вместо сравнения величины разности с 0) выполнять сравнение самих величин, участвующих в вычислении разности. 1.11. Финансовые расчеты и числа с плавающей точкой При финансовых (денежных) расчетах нельзя использовать числа с плавающей точкой. Проблема заключается в невозможности зафиксировать экспоненту в проводимых операциях с разно-порядковыми числами при ограниченности длины мантиссы. При добавлении малых величин к большим, младшие разряды могут потеряться, т. к. не попадут в мантиссу. Тип double «отодвигает» эту проблему, но не решает ее принципиально. Если в языке программирования нет типов данных с фиксированной точкой (например, в VBA для финансовых расчетов применяется тип Currency, а в VB.NET Decimal), можно выйти из положения и хранить деньги в виде целой, и 20 дробной частей как отдельных целых чисел (сбрасывая определённый старший разряд из копеек в рубли при переполнении). Иллюстрирующий пример – проект Money, представляет класс CMoney, хранящий деньги в виде целой, и дробной частей как отдельных целых чисел и класс CBadMoney хранящий деньги в виде переменной с плавающей точкой типа double. Пример демонстрирует «правильное» сложение разно-порядковых чисел в денежных операциях при использовании объектов класса СMoney и сложение с потерями при использовании объектов класса CBadMoney. Проект Round демонстрирует катастрофическую потерю точности при вычитании двух «длинных» чисел ((float)123456789, (float)123456788), имеющих разницу в один младший бит у мантисс и равные экспоненты. Указанную операцию можно интерпретировать как финансовые вычисления с неудачно выбранным типом данных в виде чисел с плавающей точкой. Необходимо отметить, что операции с плавающей точкой коммутативны, но не ассоциативны, и дистрибутивный закон для них тоже не выполняется. Например, (1 + 1020) – 1020 = 0 ≠ 1 + (1020 – 1020) = 1. 1.12. Машинное представление целых чисел В отличие от чисел с плавающей точкой, не существует формального стандарта хранения и обработки целых чисел в вычислительных системах. Работа с целыми числами значительно проще, чем с числами с плавающей точкой, и целочисленные типы данных и целочисленные операции устроены практически одинаково во всех вычислительных системах. Вычислительные системы предоставляют два класса целочисленных типов – знаковые и беззнаковые. Беззнаковые типы занимают n бит памяти и допускают в качестве значений целые числа от 0 до 2n – 1. Вычислительные системы предоставляют типы для значений n = 8, 16, 32, 64, поскольку такие типы занимают целое число байт. Максимальные представимые числа для этих размеров приведены в табл. 6. Таблица 6 2n – 1 n 8 255 16 65535 32 4 294 967 295 64 18 446 744 073 709 551 615 (≈ 2  1019 ) Значение беззнакового типа представляется в памяти двоичным числом, дополненным до n бит нулевыми битами слева. Числу 0 соответствует представление из n нулевых битов, числу (2n – 1) – из n единичных. 21 Минимальной единицей хранения чисел является байт, и значит, чтобы полностью определить представление числа в памяти, надо указать отображение n-битного числа на 8-битные байты. Определение порядка байтов важно для правильной записи структур bitfloat и bitdouble (см. проект FloatPoint), а также для понимания принципов отображения данных при использовании ассемблера, отладчика и т. п. Порядок байтов – метод записи байтов многобайтовых чисел [http://en.wikipedia.org/wiki/Endianness]. В общем случае, для представления числа M, большего 255 – максимального целого числа, записываемого одним байтом, приходится использовать несколько байтов. При этом число M записывается в позиционной системе счисления по основанию 256: M n  Ai  256i  A0  2560  A1  2561  A2  2562  ...  An  256n i 0 Набор целых чисел A0, …, An, каждое из которых лежит в интервале от 0 до 255, является последовательностью байтов, составляющих M. При этом A0 называется младшим байтом, An – старшим байтом числа M. В вычислительных системах используют два способа представления порядков байтов: от старшего к младшему и от младшего к старшему. Порядок от старшего к младшему или (англ. big-endian, дословно: «тупоконечный»): An, …, A0 – запись начинается со старшего и заканчивается младшим. Этот порядок является стандартным для протоколов TCP/IP, он используется в заголовках пакетов данных и во многих протоколах более высокого уровня, разработанных для использования поверх TCP/IP. Поэтому, порядок байтов от старшего к младшему часто называют сетевым порядком байтов (англ. network byte order). Порядок от младшего к старшему или (англ. little-endian, дословно: «остроконечный»): A0, …, An, – запись начинается с младшего и заканчивается старшим. Этот порядок записи принят в памяти персональных компьютеров с x86-процессорами, и иногда его называют интеловский порядок байтов (по названию фирмы-создателя архитектуры x86). Существенным достоинством little-endian по сравнению с big-endian порядком записи считается возможность «неявной типизации» целых чисел при чтении меньшего объёма байт (при условии, что читаемое число помещается в диапазон). Так, если в ячейке памяти содержится число 0x00000022, то прочитав его как int16 (два байта) мы получим число 0x0022, прочитав один байт – число 0x22. Однако это же может считаться и недостатком, потому что провоцирует ошибки потери данных. Обратно, считается, что у little-endian, по сравнению с big-endian есть «неочевидность» значения байтов памяти при отладке [последовательность байтов (A1, B2, C3, D4) на самом деле значит 0xD4C3B2A1, для big-endian эта последовательность (A1, B2, C3, D4) читалась бы «естественным» образом: 0xA1B2C3D4. Термины big-endian и little-endian первоначально не имели отношения к информатике. В сатирическом произведении Джонатана Свифта «Путешествия 22 Гулливера» описываются вымышленные государства Лилипутия и Блефуску, в течение многих лет ведущие между собой войны из-за разногласия по поводу того, с какого конца следует разбивать варёные яйца. Тех, кто считает, что их нужно разбивать с тупого конца, в произведении называют Big-endians («тупоконечники»). Споры между сторонниками big-endian и little-endian в информатике также часто носят характер т. н. «религиозных войн». Термины big-endian и little-endian ввёл Кохэн (англ. Danny Cohen) в 1980 году в своей статье «On Holy Wars and a Plea for Peace» – «Священные войны и призыв к миру». Что касается знаковых типов данных, то они обычно определены для тех же размеров n, что и беззнаковые. Диапазоны знаковых типов от – 2n - 1 до 2n - 1 – 1 включительно представлены в табл. 7. Таблица 7 – 2n - 1 n 2n - 1 – 1 8 – 128 127 16 – 32768 32767 32 – 2 147 483 648 2 147 483 647 64 – 9 223 372 036 854 775 808 9 223 372 036 854 775 807 Знаковые числа в памяти хранятся в дополнительном коде, что облегчает аппаратную реализацию операций сложения и вычитания для знаковых типов. При выполнении операций с целыми числами, результат которых превышает диапазон хранения чисел, происходит переполнение. Это приводит к ошибкам в программах, в том числе возможному нарушению безопасности кода. Например, при сложении двух 16-битных знаковых положительных целых можно в результате получить отрицательное число: 30000 + 30000 = −5536. При переполнении беззнаковых типов происходит округление результата по модулю числа 2n. Иллюстрирующий пример – проект IntOverflow. 1.13. Неустойчивые алгоритмы Рассмотрим задачу вычисления интегралов [4]: 1 x n x 1 e dx, n = 1, 2, … Интегрируя по частям, получим: 23 1 x n x 1 e dx  x n e x 1 |10 1   nx n 1e x 1dx, или En = 1 – n En-1 , n = 2, 3,…, где E1 = 1/e. (3) Используем это рекуррентное выражение для вычисления En (проект Algorithm) с одинарной точностью (float), двойной точностью (double) и с точностью не менее 50 десятичных знаков (boost::multiprecision::cpp_dec_float_50) (табл. 8). Таблица 8 n En (float) En (double) En (float_50) 1 0.3678795 0.3678794411714423 0.3678794411714423 2 0.2642411 0.2642411176571153 0.2642411176571154 3 0.2072767 0.207276647028654 0.2072766470286539 4 0.1708932 0.170893411885384 0.1708934118853843 5 0.145534 0.1455329405730801 0.1455329405730786 6 0.1267958 0.1268023565615195 0.1268023565615285 7 0.1124296 0.1123835040693635 0.1123835040693008 8 0.100563 0.1009319674450921 0.1009319674455933 9 0.09493256 0.09161229299417073 0.09161229298966058 10 0.05067444 0.0838770700582927 0.08387707010339416 11 0.4425812 0.07735222935878028 0.0773522288626642 12 –4.310974 0.07177324769463667 0.07177325364802956 13 57.04266 0.06694777996972334 0.0669477025756157 14 –797.5973 0.06273108042387321 0.06273216394138015 15 11964.96 0.05903379364190187 0.05901754087929777 16 –191438.3 0.05545930172957014 0.0557193459312356 17 3254453 0.05719187059730757 0.05277111916899476 18 –5.858015e+007 -0.02945367075153627 0.05011985495809426 … … … … 54 -1.#INF -2.869099451145142e+054 -0.1564473789098101 Хотя, на интервале [0, 1] подынтегральное выражение x n e x1 > 0, E12 (float) < 0, E18 (double) < 0, E54 (float_50) < 0. В чем причина ошибок? Отметим, что при вычислении начального значения E1 = 1/e невозможно избежать ошибки 24 округления, вызванной ограниченным размером мантиссы назначенного типа данных, соответственно составляющей: E1 (float) ≈ 9.14976e-09, E1 (double) ≈ 1.24288e-17, E1 (f loat_50) ≈ 1.25083e-72. Ошибку при вычислении En , вызванную ошибкой округления при вычислении E1 , можно оценить на основе (3) по формуле: En  (1) n1 n!E1 , что для конкретных значений n составит: E12 (float) ≈ –12! * 9.14976e-09 ≈ – 4.3827, E18 (double) ≈ –18! * 1.24288e-17≈= – 0.0796, E54 (f loat_50) ≈ –54! * 1.25083e-72 ≈ – 0.2887. Полученные ошибки согласуются с расчетами, представленными в строках 12, 18, 54 табл.8. Дополнительно, для проверки приведенных утверждений о влиянии ошибки округления E1 на En , был проведен расчет с повышенной точностью (50 значащих десятичных цифр – float_50) при начальной ошибке округления E1 (double) ≈ 1.24288e-17. При этом данные расчетов совпали со значениями столбца 2 табл.8, а не 3, что подтверждает превалирующее влияние на точность расчетов значения погрешности E1 , а не точности применяемой арифметики с плавающей точкой в промежуточных расчетах. Приведенный алгоритм является примером неустойчивого алгоритма. Он катастрофически увеличивает ошибки округления в начальных данных. Применение более «совершенных» типов данных (с увеличенным размером мантиссы) только отодвигает, но решает проблему. Решением проблемы является нахождение устойчивого алгоритма, который не увеличивает ошибки округления, а, в идеале, их подавляет. Сконструируем устойчивый алгоритм, решающий поставленную задачу. Перепишем формулу (3) в следующем виде: En-1 = (1 – En)/n , n = … 3, 2. (4) Теперь на каждом шаге вычисления ошибка в En делится на n. Начав с n >> 1 и вычисляя En-1 в обратном направлении, мы добиваемся уменьшения начальных 25 и промежуточных ошибок округления. Это есть пример устойчивого алгоритма. Чтобы определить начальное значение En, получим его оценку: 1 1 En   x n e x 1dx   x n dx  1  0 , n  1 n откуда следует, что взяв, например, в качестве начального значения для E100 ноль, мы совершим начальную ошибку не превышающую 0.01. При 100! вычислении En начальная ошибка уменьшится в раз. Например, при n = 54 n! 100! в ≈ 4.0428e+086 раз. Проведем вычисления En на основе устойчивого 54! алгоритма (табл. 9). Таблица 9 n En (float) En (double) En (float_50) 100 99 0.0 0.01 0.0 0.01 0.0 0.01 … … … … 54 0.01786274 … 18 … 0.05011985 … 12 … 0.07177325 … 1 … 0.3678795 0.01786274234481424 … 0.05011985495809426 … 0.07177325364802957 … 0.3678794411714423 0.01786274234481424 … 0.05011985495809426 … 0.07177325364802956 … 0.3678794411714423 После 10 вычислительных шагов по формуле (4) начальная ошибка уменьшается в 100! ≈ 6.9028e+017 раз и становится сопоставимой с ошибкой 91! округления числа типа double в рассматриваемом диапазоне значений. Поэтому дальнейшие расчеты будут выполняться с точностью возможной ошибки округления в последней цифре. Упражнение 7. Изучите алгоритм компенсационного суммирования: (http://en.wikipedia.org/wiki/Kahan_summation_algorithm). Этот алгоритм вычисления суммы последовательности чисел с плавающей точкой значительно уменьшает вычислительную погрешность по сравнению с наивным подходом. Уменьшение погрешности достигается введением дополнительной переменной для хранения нарастающей суммы погрешностей. 26 Напишите программу, реализующую данный алгоритм. Подберите данные наглядно демонстрирующие эффективность алгоритма. 1.14. Чувствительные задачи Существуют задачи, в которых «хорошие» ответы нельзя получить никаким алгоритмом, потому что сама задача чувствительна к малым ошибкам в исходных данных и ошибкам округления. Вместо чувствительности иногда используют термин обусловленность задачи, означающий меру чувствительности решения задачи к изменениям (возмущениям) её входных данных. Важно отличать чувствительные (плохо обусловленные) задачи от неудачных неустойчивых алгоритмов, примененных для решения. Анализ чувствительности задачи можно выполнить на основе вычисления коэффициентов чувствительности (влияния) исходных данных задачи, подверженных малым ошибкам, на результат решения. Пусть результат решения есть функция y от вектора данных x = [x1, x2 ,… xn]. В первом приближении, с точностью до малых второго порядка, представим функцию y отрезком ряда Тейлора в точке заданного значения вектора x0 в виде: n f |x0 xi . y = y – f(x0) ≈  i 1 dxi Здесь приращение функции y , зависит от частных производных функции, вычисленных в заданной точке x0, умноженных на малые приращения вектора данных xi , i = 1, 2, …, n. Отсюда следует, что указанные частные производные являются коэффициентами чувствительности малых изменений в данных на результат. В качестве примера чувствительной задачи рассмотрим следующий полином [4]: P(x) = (x – 1) (x – 2) … (x – 11) (x – 12) = x12 – 78 x11+ …= x12 – a x11 + … . Проведем анализ чувствительности корней r = 1, 2 , …, 12 полинома P(x) к изменению коэффициента a при x11. Полином представим функцией от a: P(x, a) = x12– a x11 + … = 0. Чувствительность будем оценивать частной производной по a для каждого корня полинома. Используя формулу производной неявной функции, получим: P x x11   a  12 12 . P a ( x  j) x i 1 j 1 j i Вычисляя последнее выражение для каждого корня, получаем: 27 x |x  i  a i11 12 , i  1, 2, ..., 12 .  (i  j ) j 1 j i Эти числа являются коэффициентами чувствительности каждого из корней к изменению в коэффициенте a. Рассчитаем их (проект Polinom) (табл. 10). Посмотрим, согласуются ли данные по коэффициентам чувствительности с погрешностями вычисления корней при малом «возмущении» коэффициента a, равном 1.0e-5. Для этого вычислим корни исходного и «возмущенного» полинома (корни вычислены в среде Matlab, файл polinom.m) (табл. 11). Расчеты показывают, что малое изменение (1.e-5) в коэффициенте a = –78 привело к тому, что 6 корней из 12, имеющих наибольшие коэффициенты чувствительности, стали комплексными, а остальные, в соответствии со своими коэффициентами чувствительности, отодвинулись от номинальных значений. Таблица 10 Корень, ri Коэффициент чувствительности, 1 -2.50521e-008 2 0.000564374 3 -0.244085 4 17.3376 5 -403.672 6 4199.04 7 -22885.7 8 71014.7 9 -129717 10 137787 11 -78624.2 12 18613.9 x | x i a Таблица 11 № Корни полинома P(x) + 1.e-5 x11 Корни полинома P(x) 12 11.999999990248927 11.671877829513981 + 0.204250367297240i 11 11.000000058400660 11.671877829513981 - 0.204250367297240i 10 9.999999847496268 9.480533999962713 + 0.746046193254474i 9 9.000000227655217 9.480533999962713 - 0.746046193254474i 8 7.999999786097445 7.365089539520506 + 0.217312596608508i 7 7.000000131091507 7.365089539520506 - 0.217312596608508i 28 6 5.999999947608275 5.961076403539962 5 5.000000013198870 5.004081689151476 4 3.999999998067797 3.999826734133822 3 3.000000000137872 3.000002440821697 2 1.999999999997095 1.999999994358464 1 0.999999999999972 1.000000000000223 Средствами «борьбы» с чувствительными задачами является повышение точности машинной арифметики при выполнении расчетов участков программы, где возможна потеря точности. Например, можно использовать типы данных cpp_dec_float_50 и cpp_dec_float_100 библиотеки Boost. В качестве упражнения проанализируем известный итерационный алгоритм Герона, вычисляющий квадратный корень из положительного числа «a»: 1 a . (5)  xn   xn 1  2 xn 1  Формулу (5) можно получить, применяя метод Ньютона к решению уравнения f ( x )  a  x 2  0 в виде: f ( xn 1 )  f ( xn 1 )( xn  xn 1 )  a  xn21  2 xn 1( xn  xn 1 )  0 . Докажем, что алгоритм является устойчивым и нечувствительным к погрешностям в исходных данных и ошибкам округлений. Определим значения для приближений: xn  a (1   n ), xn 1  a (1   n 1 ) , (6) где  n ,  n 1 – относительные погрешности вычисления квадратного корня на шаге n и n – 1. Подставив (6) в (5), после преобразований получим: 1   n21  .  n   2   n 1  1  (7) Принимая в качестве начального приближения x0  0 , из (5), (6) получаем x  n 1 1 и x1, x2 ,..., xn 1  0 и (1   n 1 )  n 1  0 ; из (7)  n 1  0 ,  n 1  1 a 1  n   n 1 . Таким образом,  n убывает с ростом n на каждом шаге не менее 2 чем в 2 раза и xn  a при n   . В итоге, взяв произвольное положительное число в качестве исходных данных мы, тем не менее, получаем правильный ответ, что означает устойчивость рассматриваемого алгоритма. 29 Оценим влияние ошибок округления, используя вместо (5) формулу: 1 a  1   M  , xn   xn 1  2 xn 1  (8) где  M – относительная ошибка округления на шаге n, имеющая порядок значения машинного эпсилон. Подставляя в (8) формулы (6) получаем:  1 a (1   M )  (1   n ) a   (1   n 1 ) a  2 a (1   n 1 )  1   n21  1  2 n 1  2  1   n21  (1   M )   (1   M )   M . (1   M )  1    n   2   n 1  1  2   n 1  1  2   n 1  1  Из последней формулы следует, что  n   M при n   , т.е. погрешность вычисления a не превышает погрешности округления и, в итоге, алгоритм является устойчивым и нечувствительным как к погрешностям в исходных данных, так и к ошибкам округлений. Упражнение 8. Напишите программу, реализующую алгоритм Герона, вычисляющий квадратный корень из положительного числа. Упражнение 9. Постройте алгоритм и напишите программу, вычисляющую корень целой степени n > 0 из положительного числа. Упражнение 10. Вычитание близких чисел может значительно увеличить относительную погрешность вычисления из-за ошибок округления. Для многих распространенных математических формул есть алгоритмы вычисления, позволяющие значительно уменьшить погрешности округлений. Докажите, что относительная ошибка расчета формулы (x-y)(x+y) меньше, чем у (x2-y2). Полезную информацию по этой теме можно найти в [5]. Упражнение 11. Укажите потенциальные проблемы, связанные с потерей точности вычисления корней квадратного уравнения из-за ошибок округления. Проведите анализ чувствительности задачи вычисления корней квадратного уравнения. Разработайте алгоритм вычисления корней квадратного уравнения с максимальной точностью. Полезную информацию по этой теме можно найти в [5,6]. 1.15. Сложность алгоритмов Предположим, что выполнен анализ чувствительности задачи, и надо выбрать алгоритм для ее решения. Безусловно, нам необходим устойчивый алгоритм. Но кроме этого, алгоритм должен быть эффективным, т.е. обладать минимально возможной или приемлемой степенью сложности. 30 Например, как было показано ранее, множество чисел с плавающей точкой является конечным, и потому мы можем вычислить приближённые значения корней некоторой функции (или убедиться в их отсутствии), перебрав все числа и вычисляя в них значения функции. Но, будучи устойчивым, такой алгоритм не эффективен и для практики бесполезен. Сложность алгоритма будем оценивать временем выполнения и требуемой памятью для его реализации при различных входных данных. Параметры время и память часто бывают взаимозависимы. Например, в задаче поиска в массиве элементов, время поиска можно сократить, создав индексы сортировки элементов, что увеличивает объем необходимой алгоритму памяти. В свою очередь, чрезмерное увеличение требуемой алгоритму памяти может негативно сказаться на быстродействии программы, например, из-за необходимости использования виртуальной памяти, хранящейся на внешних носителях. В дальнейшем, будем рассматривать временную сложность алгоритмов, за исключением тех случаев, когда размер используемой памяти непосредственно влияет на возможность решения задачи. Нахождение точной зависимости времени t(n) от размера исходных данных n для конкретной программы сложно и вряд ли необходимо. Поэтому используют асимптотические оценки t(n) при «большом» n, в пределе стремящимся к бесконечности. Примером n в задачах поиска, сортировки и других, связанных с обработкой наборов объектов, может быть размер набора. Для задачи вычисления корней функции это может быть количество корней. Оценивать t(n) целесообразно не по астрономическому времени выполнения, зависящему от используемого «железа», а по количеству базовых операций, вносящих наибольший вклад в общее время выполнения алгоритма. Базовыми операциями, как правило, являются наиболее длительные операции, выполняемые во внутреннем цикле алгоритма. Необходимо отметить, что существуют алгоритмы, у которых время выполнения зависит не только от размера входных данных, но и от их значений. Например, у алгоритма последовательного поиска завершение может произойти как вначале, так и в конце просмотра списка. Поэтому эффективность алгоритма можно оценивать в наихудшем, наилучшем или среднем случаях. Анализ наихудшего случая важен, поскольку позволяет получить оценку эффективности алгоритма сверху. Анализ эффективности в среднем, позволяет получить оценку для типовых или случайных данных и очень важен при рассмотрении алгоритмов, эффективность которых гораздо выше для среднего, чем для наихудшего случаев. Рассмотрим алгоритм простой сортировки (сортировка пузырьком) с эффективностью t(n) = кС(n), где к – постоянный множитель, равный времени выполнения базовой операции, С(n) =n(n-1)/2 – количество базовых операций. Как изменится эффективность алгоритма при увеличении размера данных в m раз? В асимптотике, при n, стремящимся к бесконечности, справедлива формула t(mn)/t(n) ≈ m2 , которая не зависит ни от к, ни от других постоянных множителей. По этим причинам в процессе анализа эффективности оценивают 31 порядок роста количества основных операций с точностью до постоянного множителя. Порядок роста вычисляют для больших объемов данных, потому что при малых объемах, как правило, не удается найти разницу между сравниваемыми алгоритмами. Для того чтобы можно было сравнивать порядки роста в математике введено обозначение O – прописная «O». Обозначение O(g(n)) – это множество всех функций, порядок роста которых при достаточно больших n не превышает некоторую константу, умноженную на значение функции g(n). Более строгое определение O(g(n)) заключается в следующем. Говорят, что функция t(n) принадлежит множеству O(g(n)), если для всех больших значений n функция t(n) ограничена сверху некоторой константой, умноженной на функцию g(n), т. е существует положительная константа c и неотрицательное целое число n0 такое, что t(n) < c g(n) для всех n ≥ n0 [7] . Справедливо следующее утверждение. Если и t1 (n)  O( g1 (n)) t2 (n)  O( g2 (n)) , то t1 (n)  t2 (n)  O(max{ g1 (n), g 2 (n)}) . Это означает, что общая эффективность алгоритма зависит от его части, имеющей наибольший порядок роста. Например, проверить содержит ли массив повторяющиеся элементы можно следующим двухэтапным алгоритмом. На первом этапе сортировка, на втором проверка рядом стоящих дубликатов. Если алгоритм сортировки выполняется за n(n-1)/2 базовых операций, т. е. принадлежит O(n2), а алгоритм 2-го этапа выполняется за n – 1 операций сравнения, т. е. принадлежит O(n), то суммарная эффективность общего алгоритма будет принадлежать множеству O(n2). Приведем еще одно полезное для анализа утверждение. Если во внутреннем цикле одного алгоритма вызывается другой алгоритм, то общая эффективность t1 (n)  t2 (n)  O( g1 (n))  O( g 2 (n)) . При сравнении порядков двух функций t(n) и g(n) существует удобный метод, основанный на вычислении предела их отношения: t (n) lim g (n) . n  Если рассматриваемый предел равен 0, то t(n) имеет меньший порядок роста, чем g(n). Если рассматриваемый предел равен константе, то t(n) имеет тот же порядок роста, что и g(n). При вычислении пределов можно воспользоваться, например правилом Лопиталя: t (n) t ' (n) ,  lim ' lim n    g (n) n    g (n) формулой Стирлинга: n n!  2n ( )n для достаточно больших n, e и др. Для примера, сравним порядки роста функций log2n и n , функций n! и n 2 : 32 ' log 2 n (log 2 n )  lim  lim ' n n n ( n ) n lim (log 2 e) 1 1 n  0; 2 n n 2n ( )n n! n n lim 2n  lim 2n e  lim 2n ( 2e )   . n  n  n  Перечислим основные функции, которые используются при вычислении эффективности алгоритмов, ранжируя их по возрастанию сложности: 1. log2(n) – логарифмическая (поиск в отсортированном массиве, метод половинного деления при нахождении корня функции); 2. n – линейная (поиск в несортированном массиве, поиск корня функции прямым перебором); 3. n log2(n) – «n-log-n функция» (сортировка слиянием, быстрая сортировка); 4. n2 – квадратичная (простые алгоритмы сортировки, содержащие два вложенных цикла, ряд операций над матрицами размером n x n); n 5. 2 – экспоненциальная (например, алгоритм, в котором могут участвовать все натуральные числа размером n бит); 6. n! – факториальная (алгоритм, выполняющий обработку всех перестановок некоторого множества, состоящего из n элементов). Переборные задачи из 5, 6 пунктов являются самыми трудоемкими и получили название задач неполиномиальной трудности (NP-трудных задач). Эти задачи изучает теория NP-трудности. Теория часто не отвечает на вопрос о трудоемкости тех или иных задач, но позволяет утверждать, что некоторые задачи имеют одинаковую трудность, как другие известные задачи. Если известно, что некоторая задача не проще, чем известные переборные задачи, то имеет смысл и для рассматриваемой задачи применить алгоритм переборного типа, имеющий экспоненциальную трудоёмкость. В заключение, отметим, что существуют методы экспериментальной оценки сложности, основанные на генерации типовых, случайных наборов входных данных разных размеров и оценки функции, аппроксимирующей зависимость времени вычисления от размера данных [7]. Упражнение 12. Для шести представленных функций, которые используются при оценке эффективности алгоритмов, вычислить их значения 2 3 4 5 6 при размерах входных данных равных: 10, 10 , 10 , 10 , 10 , 10 . Упражнение 13. Оцените сложность алгоритма Герона, вычисляющего квадратный корень из положительного числа. 33 Дополнительные упражнения Упражнение 14. Вычислите методами половинного деления и Ньютона: arcsin(a), arccos(a), arctan(a) на ЭВМ, умеющей вычислять только функции sin(), cos(), tan() соответственно. Упражнение 15. Вычислите (с наивысшей точностью для выбранного типа данных: float, double) методами прямоугольников, трапеций, парабол 1 (Симпсона) интеграл  exp(  x 2 )dx . Примечание. Число разбиений интервала интегрирования не должно задаваться, а должно вычисляться автоматически для достижения наивысшей точности. Для этого: 1. Берем вначале 1 часть (диапазон от 0 до 1) и вычисляем интеграл - s1. 2. Увеличиваем количество частей в 2 раза и вычисляем интеграл - s2. 3. Сравниваем s1 и s2 (как числа с плавающей точкой):  если они "равны", то прекращаем расчет, принимая в качестве значения интеграла - s2;  если нет, то s1=s2 и переходим к пункту 2. 34 2. Математические пакеты программ Цель следующих параграфов – научить читателя моделировать и решать вычислительные задачи в математическом пакете программ GNU Octave. GNU Octave – это свободно распространяемая система для математических вычислений, использующая совместимый с Matlab (Matrix Laboratory) встроенный язык программирования. Matlab широко используемый коммерческий математический пакет. При соблюдении некоторых простых правил, описанных в [http://ru.wikipedia.org/wiki/GNU_Octave, раздел «Совместимость с Matlab»], программы Octave будут совместимы c Matlab. 2.1. Установка GNU Octave Рассмотрим установку GNU Octave на примере версии Octave3.6.4_gcc4.6.2 для платформы Windows. Для установки необходимо выполнить следующую последовательность действий: 1. Выбрать диск для установки, например, C:. 2. Создать установочный каталог, название которого не должно содержать пробелы, например, C:\Octave. 3. Разархивировать в установочный каталог основной архив Octave3.6.4_gcc4.6.2_20130408.7z, например, с помощью бесплатной программы 7-zip, так, чтобы образовалась следующая иерархия для каталога Bin: C:\Octave\Octave3.6.4_gcc4.6.2\bin. Аналогично bin, располагаются все остальные каталоги: doc, etc и т. д. 4. Разархивировать в установочный каталог дополнительные Octave-forge пакеты из архива Octave3.6.4_gcc4.6.2_pkgs_20130402.7z, сохраняя иерархию каталогов: C:\Octave\Octave3.6.4_gcc4.6.2\bin и др. 5. Запустить пакет Octave (ярлык для запуска – Octave3.6.4_gcc4.6.2.lnk) и в соответствии с инструкцией из http://wiki.octave.org/Octave_for_Windows выполнить 5 команд: 1. 2. 3. 4. 5. pkg rebuild -auto pkg rebuild -noauto ad % may crash octave when loaded and 'clear all' is executed pkg rebuild -noauto nan % shadows many statistics functions pkg rebuild -noauto gsl % shadows some core functions pkg rebuild -auto java 6. Для удобства в работе можно отредактировать содержимое файла C:\Octave\Octave3.6.4_gcc4.6.2\share\octave\site\m\startup\octaverc добавив строки, определяющие кодовую страницу Windows (win1251)  system("chcp 1251"); и текущий рабочий каталог  cd("c:/Octave/work/numerical_methods"). 35 Когда потребуется переустановить актуальную на текущую дату версию Octave для Windows, воспользуйтесь указаниями по ссылке http://wiki.octave.org/Octave_for_Windows. 2.2. Введение в GNU Octave В следующих параграфах основное внимание уделим работе с командным интерпретатором и встроенному языку программирования Octave. Подробное изложение возможностей Octave и примеры его использования можно найти в книге [8]. Для запуска пакета GNU Octave используйте ярлык Octave3.6.4_gcc4.6.2.lnk, который (после установки пакета) находится в папке C:\Octave. При запуске пакета создается консольное окно Octave, предназначенное для ввода и интерпретации команд встроенного языка Octave. В Octave есть также возможность обработки групп команд, подготовленных в виде файлов с расширением .m (m-файлов). Последнее позволяет создавать на встроенном языке Octave программы, выполняющие пользовательские алгоритмы, использующие арсенал реализованных в составе пакета и его расширений функций. Справку по встроенным функциям можно получить по команде: help <имя_функции>. Например, команда «help fzero» выводит справку по встроенной функции вычисления корней уравнений. При вводе команд учитывается регистр. Выход из пакета осуществляется командами exit или quit. 2.3. Представление данных в Octave Octave, как и Matlab, работает с единственным видом объектов – матрицами, элементами которых могут быть действительные или комплексные числа, символьные переменные, участвующие в символьных вычислениях, а также текстовые символы. Скаляр – это матрица размером 1  1. Матрица, состоящая из одной строки или столбца, является вектором. Комментарий в Octave предваряется символами «%» или «#». Для совместимости с Matlab используйте символ «%». Ввод данных может быть выполнен через присвоение значений элементам соответствующих матриц. Элементами матриц могут быть любые выражения Octave. Квадратные скобки используются для формирования матриц и векторов. В качестве разделителей элементов в строках матриц используются пробелы или запятая. Точка с запятой отделяет строки матриц. Например, [6.9 9.64 sqrt(-1)] – вектор с тремя элементами, разделенных пробелами, а [6.9, 9.64, sqrt(-1)] – тот же самый вектор, где в качестве разделителей использованы запятые. [1+i 2-i 3] и [1 +i 2 -i 3] разные векторы с комплексными числами 36 (символы i, I, j, J используются для обозначения мнимой единицы). Первый состоит из трех элементов, второй из пяти. В матрице [1 2 3; 4 5 6], размером 2x3, точка с запятой отделяет строки. Векторы и матрицы могут использоваться внутри квадратных скобок ([ ]). Запись [A B; C] допустима, если количество строк А равно количеству строк B и в сумме количество столбцов A и B составляет количество столбцов матрицы C. Это основное правило для правильного конструирования матриц. A = [] – запоминание пустой матрицы в переменной A. Полезной операцией при генерации векторов и при работе с индексами матриц является двоеточие. Запись m:k генерирует вектор-строку [m, m+1, …, k]. Вектор будет пустым, если m > k. Запись m:h:k генерирует вектор-строку [m, m+h, m+2h, …, k]. Вектор будет пустым, если h > 0 и m > k или h < 0 и m < k. Запись A(:) есть вектор-столбец, составленный из всех элементов матрицы A. Если A(:) стоит слева от оператора присваивания, то это приводит к заполнению матрицы A данными, генерируемыми выражением, стоящим справа от оператора присваивания, в соответствии с её размерами. Запись A(:, j) есть j-й столбец, а A(j, :) j-я строка матрицы A. Круглые скобки используются для назначения индексов матриц. Так, если X и V векторы, то X(V) есть [X(V(1)), X(V(2)), …, X(V(n))]. Компоненты вектора V должны быть целыми положительными числами и используются как индексы. Например, при V=[1, 2, 3] – X(V) первые три элемента вектора X. Выражение X(length(X) : -1 :1) осуществляет перестановку элементов вектора X в обратном порядке. Если величины индексов меньше 1 или больше размерности вектора X, генерируется сообщение об ошибке. Подобное индексирование используется и в матрицах. Пусть вектор V имеет m компонент, а вектор W имеет n компонент. Тогда матрица A(V,W) размером m  n будет сформирована из элементов матрицы A, индексы которых определяются значениями элементов векторов V и W. Например, выражение A([1,5], :) = A([5,1], :) меняет строки 1 и 5 матрицы A местами. 2.4. Функции формирования матриц Приведем основные функции, формирующие матрицы или вычисляющие их параметры. Функция zeros(N) формирует матрицу размером N  N, состоящую из 0. Функция zeros(M,N) формирует матрицу нулей размером M  N. Функция ones(N) формирует матрицу размером N  N, состоящую из 1. Функция ones(M,N) формирует матрицу из единиц размером M  N. Функция eye(N) формирует единичную матрицу размером N  N. Функция eye(M,N) формирует матрицу размером M  N, состоящую из 1 по диагонали и 0 вне её. 37 Функция diag(X) формирует квадратную диагональную матрицу с вектором X на главной диагонали. Функция diag(X[, k]) формирует квадратную диагональную матрицу с вектором X на k-ой диагонали. Функция diag(A[, k]), где A матрица, формирует вектор-столбец, составленный из элементов главной (k = 0), k-ой (k > 0) по порядку вверх, k-ой (k < 0) по порядку вниз диагонали матрицы A. Функция rand([M,N, ... ]) формирует матрицу с элементами, распределенными по равномерному закону в диапазоне (0, 1). Функция rand без аргументов возвращает одно случайное число. Функция randn([M,N, ... ]) формирует матрицу с элементами, распределенными по закону Гаусса с нулевым средним и единичной дисперсией. Функция randn без аргументов возвращает одно случайное число. Функция repmat(A, N[, M]) формирует матрицу, составленную из N  N или N  M копий матрицы A. Функции tril(A[, k]), triu(A[, k]) формируют из матрицы A, соответственно, нижнюю или верхнюю треугольную матрицу начиная с главной или k-ой диагонали. Функция [m, n] = size(A) определяет число строк (m) и столбцов (n) матрицы A. Функция sum(A[, k]) – формирует вектор-строку (k = 1) или вектор-столбец (k = 2), каждый элемент которого является, соответственно, суммой соответствующего столбца или строки матрицы A. Функция sum(A) вычисляет сумму столбцов. Выражение sum(sum(A)) вычисляет сумму всех элементов матрицы A. Функция diff(A) – формирует матрицу, количество строк которой на 1 меньше количества строк матрицы A. Строки сформированной матрицы представляют разность соседних строк матрицы A. При этом соседняя строка с младшим номером вычитается из строки со старшим номером. Функция max(X) возвращает наибольший элемент вектора X. Для матриц, [Y, I] = max(A) – вектор Y, состоит из максимальных элементов каждого столбца матрицы A, вектор I, состоит из номеров строк, в которых располагаются максимальные элементы столбцов. Функция max(A,B) возвращает матрицу такого же размера, как A и B, с наибольшими элементами, взятыми из A или B. Вызов max(A,[],k) позволяет управлять поиском. При k = 1 ищутся максимальные элементы столбцов, а при k = 2 строк матрицы A. Выражение max(max(A)) возвращает максимальный элемент матрицы A. Описание функции «max» применимо для функции «min» при поиске минимальных элементов матрицы вместо максимальных. Функция mean(A, [k]) – формирует вектор-строку (k = 1) или векторстолбец (k = 2), каждый элемент которого является, соответственно, средним значением столбца или строки матрицы A. Функция mean(A) вычисляет 38 средние значения столбцов матрицы A. Выражение mean(mean(A)) вычисляет среднее значение всех элементов матрицы A. Функция [Y, I] = sort(A) возвращает матрицу Y, состоящую из столбцов матрицы A, отсортированных по возрастанию. Матрица I состоит из столбцов, содержащих номера строк, расставляемых по порядку элементов столбцов матрицы A, так что A(I(:,k),k)=Y(:,k), где k – номер столбца матриц A и Y. Все введенные или вычисленные переменные-матрицы хранятся в рабочей области Octave, существующей до завершения работы Octave. 2.5. Сервисные функции Функции who, whos формируют список текущих переменных из рабочей области. Функция whos печатает также размеры переменных. Функция clear [x, …] удаляет все или указанные переменные из рабочей области. Функция clear functions удаляет все скомпилированные m-функции. Функции save и load используются для сохранения в файле и считывания из файла переменных рабочей области Octave. Функции pwd, dir, cd работают с файловой системой. Функция pwd отображает имя текущего каталога, функция dir содержимое каталога. Вызов функции cd path делает текущим каталог с именем path. 2.6. Функции матричных вычислений Функция det(A) вычисляет определитель квадратной матрицы A. Функция trace(A) вычисляет след матрицы A, равный сумме элементов её главной диагонали. Функция norm(A[, k]) вычисляет нормы матрицы A. При k = 1, вычисляется 1-я норма. Эта норма равна максимальному значению вектора, составленного из сумм абсолютных величин каждого столбца матрицы A. С помощью рассмотренных ранее функций эта норма может быть вычислена следующим выражением: max(sum(abs(A))). При k = 2, вычисляется 2-я норма. Эта норма равна максимальному сингулярному числу матрицы A и может быть вычислена следующим выражением: max(svd(A)), где функция svd выполняет сингулярное разложение (описание svd см. далее). При k = inf, вычисляется бесконечная норма. Эта норма равна максимальному значению вектора, составленного из сумм абсолютных величин каждой строки матрицы A. С помощью рассмотренных ранее функций эта норма может быть вычислена следующим выражением: max(sum(abs(A),2)). При k = 'fro', вычисляется Евклидова норма. Эта норма равна корню квадратному из суммы квадратов элементов матрицы A и может быть вычислена следующим выражением: sqrt(sum(sum(A.^2))), где «.^» операция поэлементного возведения в степень матрицы A. 39 Функция cond(A[, k]) вычисляет число обусловленности (большее или равное 1) матрицы A, используя соответствующую норму матрицы A – norm(A[, k]). Рассматривая системы линейных уравнений Ax  b и A( x  x)  b  b , где b – вектор ошибки вектора правой части линейного уравнения, а x – соответствующий вектор ошибки вектора решения линейного уравнения, получаем знаком x x  cond ( A) b b . В последнем неравенстве обозначены нормы соответствующих векторов. Величина относительное изменение правой части, а величина x x b b есть есть относительная ошибка, вызванная этим изменением. Указанное неравенство показывает, что число обусловленности исполняет роль коэффициента чувствительности (обусловленности) получаемого решения, поскольку изменения правой части могут повлечь за собой изменения в решении большие в cond(A) раз [4, с. 54]. Функция inv(A) вычисляет матрицу обратную A. Функция eig(A) вычисляет собственные значения квадратной матрицы A. Вызов функции [V, L] = eig(A) вычисляет матрицу V, столбцы которой есть собственные векторы матрицы A, и диагональную матрицу L, содержащую её собственные значения. Собственный вектор квадратной матрицы – это вектор, который будучи умноженным на матрицу даёт коллинеарный вектор, т. е. тот же вектор, умноженный на некоторое скалярное значение, называемое собственным значением матрицы: A ·V(:, i) = L(i, i) ·V(:, i), где i – номер собственного значения L(i, i), соответствующего собственному вектору V(i, :). Функция rref(A) выполняет приведение квадратной матрицы A к треугольному виду методом исключения Гаусса. Функция [L, U, P] = lu(A) вычисляет матрицы: L – нижняя треугольная с единицами на главной диагонали, U – верхняя треугольная матрица общего вида и P – матрица перестановок, получаемая из единичной матрицы путём перестановки строк или столбцов, так что P ·A = L ·U. Такое разложение можно осуществить для любой невырожденной матрицы для вычисления обратной матрицы по компактной схеме или решения системы линейных уравнений. Функция [Q, R, P] = qr(A) выполняет QR – разложение, вычисляя ортогональную матрицу Q, верхнюю треугольную матрицу R и матрицу перестановок P, так что A ·P = Q ·R. Матрица перестановок P выбирается так, чтобы элементы вектора abs(diag(R)) располагались по убыванию. QR – разложение применяется для решения переопределенных систем линейных уравнений методом наименьших квадратов. Функция chol(A) выполняет разложение Холецкого положительноопределенной матрицы A в виде , где – нижняя треугольная матрица со строго положительными элементами на диагонали. Разложение 40 Холецкого всегда существует и единственно для любой симметричной положительно-определённой матрицы. Это разложение может применяться для решения системы линейных уравнений , если матрица симметрична и положительно-определена. Такой способ решения иногда называется методом квадратных корней. По сравнению с более общими методами, такими как метод исключения Гаусса или LU-разложение, метод квадратных корней устойчивее численно и требует примерно вдвое меньше арифметических операций. Функция [U, S, V] = svd(A) выполняет сингулярное разложение матрицы A размером M  N, вычисляя U – ортогональную M  M матрицу, V – ортогональную N  N матрицу, и S – диагональную матрицу размером M  N, у которой Sij = 0, если i не равно j, и Sij  0, если i равно j. Положительные числа, стоящие на диагонали матрицы S называют сингулярными числами матрицы A, а столбцы матриц U и V – левыми и правыми сингулярными векторами. При этом A=U ·S ·VT. Сингулярное разложение – метод матричной факторизации, выполняющий наиболее надежное решение систем линейных уравнений методом наименьших квадратов. Надежность понимается в плане учета ошибок в исходных данных и округления, а также учета возможной линейной зависимости уравнений [4]. 2.7. Встроенный язык программирования Рассмотрим язык программирования Octave. На этом языке составляются программы Octave, размещаемые в m-файлах. Для выполнения программы достаточно в среде выполнения Octave ввести имя m-файла без расширения. Простейшим оператором языка является оператор присваивания: имя_переменной = выражение; или просто: выражение;. В последнем случае создается специальная встроенная переменная с именем ans. Её значение будет переопределено, если выполнится следующее выражение без присваивания значения конкретной переменной. Если в конце оператора присваивания опустить символ «;», то значение вычисленной переменной будет выведено на консоль. Имена переменных должны начинаться с буквы. Необходимо учитывать регистр. Например, A и a – это разные переменные. По умолчанию, числовое «выражение» вычисляется с удвоенной точностью, как тип double. Функция single(x) конвертирует переменную «x» к типу данных с одинарной точностью (float). 41 Octave поддерживает работу со значениями чисел с плавающей точкой, определенными в стандарте IEEE 754: inf (бесконечность, например, получается в результате 1/0) и nan («не число», например, получается в результате 0/0). Функции isinf и isnan позволяют определить элементы матриц, принимающих значения inf или nan. Встроенная переменная eps определяет машинный эпсилон числа типа double. Функции format short (формат установлен по умолчанию), format long и др. позволяют управлять форматом выводимых чисел с плавающей точкой. Так short устанавливает вывод с 5-ю значащими цифрами, long c 15. Функция disp(X) печатает значения матрицы с новой строки. Например, disp("The value of pi is:"), disp (pi) напечатает значение числа пи. В выражениях Octave определены арифметические матричные операции: + сложение; – вычитание;  умножение; \ левое деление; / правое деление; ^ степень. Сложение, вычитание, умножение выполняются по правилам, принятым для матриц. Матричное деление слева A\B, приближенно можно рассматривать как inv(A)·B, где функция inv обращает матрицу A. Если A матрица размером N  N и B матрица с одним или несколькими столбцами размером N, то x = A\B есть решение системы линейных уравнений A·x = B, вычисленное методом исключения Гаусса. Если матрица A плохо обусловлена, то печатается предупреждение. Если A матрица размером M  N с M меньшим или большим N и B матрица с одним или несколькими столбцами размером M, то то x = A\B есть решение методом наименьших квадратов, соответственно, недоопределенной или переопределенной системы линейных уравнений A·x = B. Правое матричное деление B/A, приближенно можно рассматривать как B·inv(A) и все сказанное про левое деление справедливо для правого деления, поскольку B/A = ( A'\B')', где « ' » операция транспонирования матриц. Возведение в степень A^p, где p скаляр. Если p целое число большее 1, то степень вычисляется последовательным умножением. Для других величин «p», вычисление производится с использованием техники разложения матрицы A на матрицы собственных значений и векторов, возведением в степень p диагональных элементов матрицы собственных значений и умножению слева и справа матрицы собственных значений на матрицу собственных векторов. Кроме матричных операций существуют поэлементные арифметические операции : .+ сложение; .– вычитание; .  умножение; .^ степень. Например, результатом операции А.  B будет умножение соответствующих элементов матриц между собой, а результатом операции А.^p поэлементное возведение в степень. Существует 6 операций отношения, которые могут использоваться при сравнении матриц одного размера: < – меньше, <= – меньше или равно, > – больше, >= – больше или равно, == – равно, ~= – не равно. Сравнение осуществляется между парами соответствующих элементов и формирует матрицу из 0 и 1, где 1 соответствует истине, а 0 ложности результата сравнения. Функция find позволяет находить ненулевые элементы матриц. Так 42 для вектора X, функция find(X < 3) выдаст вектор, элементами которого являются индексы всех элементов вектора X, меньшие 3. В Octave определены три логические поэлементные операции, использующиеся в основном на 0-1 матрицах: & – логическое И, | – логическое ИЛИ, ~ – логическое НЕ. Эти операции рассматривают любые ненулевые значения как истинные, возвращают 0-1 матрицы и работают в соответствии с правилами логических операций. Функции any и all применяются в сочетании с логическими операциями. Если X – 0-1 вектор, any(X) возвращает 1, если в X есть ненулевой элемент и 0, если X состоит из одних 0. Функция all(X) возвращает 1, если все элементы X ненулевые, и 0, если есть хотя бы один нулевой элемент. Эти функции бывают особенно полезны в операторах управления (например if), поскольку операторы управления проверяют скалярные, а не векторные условия. Для аргументов-матриц функции any и all действуют по столбцам и возвращают вектор-строку с результатом для каждого столбца. Управляющие конструкции языка Octave состоят из операторов условий if и switch, операторов цикла while и for. Структура оператора if: if условие 1 операторы 1 elseif условие 2 операторы 2 … elseif условие n операторы n else операторы end Работает оператор так. Если условие 1 истинно, то выполняются операторы 1, иначе проверяются все условия по порядку и для истинного условия выполняется соответствующий оператор. Если ни одно из if, elseif условий не выполнено, то выполняются операторы ветки else. Структура оператора switch: switch параметр case значение 1 операторы 1 … case значение n операторы n otherwise операторы end 43 Выполняется оператор той ветки case, значение которой совпадает со значением параметра в условии switch. При несовпадении значения параметра со значениями веток case, выполняются операторы ветки otherwise. Структура оператора цикла с предусловием while: while выражение операторы end Работает цикл следующим образом. Вычисляется «выражение». Пока оно истинно выполняются операторы, иначе управление передается оператору, следующему за end. Выражение вычисляется перед каждой итерацией цикла. Структура оператора цикла с известным количеством повторений for: for переменная_цикла = начальное_значение : шаг : конечное_значение операторы end Цикл выполняется пока значение переменной_цикла меньше или равно конечному_значению. На каждой итерации цикла переменная_цикла складывается со значением шага. Для досрочного завершения циклов while и for используется оператор break, а для досрочного перехода к следующей итерации оператор continue. Рассмотрим правила записи программ Octave, хранящихся в m-файлах. Как отмечалось ранее, есть два типа m-файлов: m-файлы сценарии и m-файлы функции. Сценарий – это последовательность операторов или команд Octave. В процессе интерпретации сценария могут вызываться встроенные в Octave или составленные пользователем m-файлы функции. Областью видимости переменных, созданных в процессе выполнения сценария, является рабочая область Octave – «глобальный» пул памяти, существующий с момента запуска и до завершения работы Octave. Сценарии позволяют автоматизировать выполнение последовательностей команд Octave или служат внешним интерфейсом, сохраняющим переменные в рабочей области Octave, при вызове функций, хранящихся в m-файлах функциях. Файлы-функции хранят разработанные пользователем функции, расширяющие функционал Octave. По правилам Octave, имя функции должно совпадать с именем m-файла, в котором она записана. В отличие от сценария, функция Octave имеет определенную структуру. Тело функции помещается между заголовком функции и необязательным оператором завершения end. Заголовок функции записывается в виде: function [Z1,Z2, …, Zn] = имя_функции(X1,X2, …, Xm), 44 где Z1,Z2, …, Zn – выходные параметры функции, X1,X2, …, Xm – входные параметры функции. В m-файле функции кроме основной функции, имя которой совпадает с именем файла и имеет глобальную область видимости, могут присутствовать внутренние функции с локальной (в рамках этого mфайла) областью видимости их имен. Переменные, определенные в теле функции, имеют локальную для этой функции область видимости и не сохраняются в рабочей области Octave после завершения выполнения функции. В заключении отметим, что Octave поддерживает рекурсивный вызов функций. 2.8. Графика в Octave Octave позволяет строить двумерные и трехмерные графики в декартовой и полярной системе координат. Мы рассмотрим работу с двумерными графиками в декартовой системе координат. Более подробную информацию можно получить в [8]. Простейший вид функции plot(X1, Y1, S1, …, Xn, Yn, Sn) строит в декартовой системе координат множество графиков зависимостей Yi как функции от Xi, i = 1, …, n, используя строки Si для обозначения типа, цвета линии и при необходимости «легенды». Символы, отвечающие за тип и цвет, следуют друг за другом без разделителя, а легенда отделяется от них точкой с запятой и завершается точкой с запятой. Весь набор символов берется в кавычки, формируя строковую константу. За сплошную линию графика отвечает символ «-», за точечную символы маркера точки: «.», «*», «+», «o» и т. д. (см. help на функцию plot). Цвет символа определяется буквами латинского алфавита: y – желтый, m – розовый, c – голубой, r – красный, g – зеленый, b – синий, w – белый, k – черный. Порядок следования символа типа и цвета не важен. Если одновременно присутствует, как символ сплошной линии («-»), так и символ маркера, то линия будет сплошная, а точки будут помечены символом маркера. Если символ сплошной линии отсутствует, то при наличии маркера линия будет точечная, а без символа маркера сплошная. За размер маркера отвечает параметр «markersize» функции plot. Например, вызов функции plot(x=0:0.5:10, x.^2, “-*;x^2;”, “markersize”, 4) нарисует график параболы, у которой точки размером 4 помечены символом «*» и соединены сплошной линией, а в правом верхнем углу присутствует описание графика в виде легенды. Приведем ряд полезных команд, управляющих графикой Octave. Команда grid on рисует на графическом окне сетку, а команда grid off её стирает. Вызов команды grid без параметра работает как переключатель присутствия или отсутствия сетки на графике. Команда hold on включает режим удержания текущего графика в графическом окне. В этом режиме следующие команды plot будут добавлять сформированные графики к текущему с использованием 45 уже существующего масштаба осей. Команда hold off выключает режим. Вызов команды hold без параметра работает как переключатель режима удержания графика. Функции title(“graph title”), xlabel(“axis X”), ylabel(“axis Y”) – выводят соответственно заголовок графиков, подписи под осями X и Y. Функция text(X,Y,“text”) выводит определенный “text” справа от точки с координатами X,Y. Функция subplot(m,n,p) разбивает графическое окно на m  n подокон и выбирает p-е окно для отображения текущего графика. Например, команды: subplot(2,1,1); plot(x=1:10,x.^2); subplot(2,1,2); plot(x=pi:0.5:pi,sin(x)) рисуют график параболы в верхнем подокне, а график синуса в нижнем. Вызов функции subplot(1,1,1) без параметров возвращает конфигурацию, установленную по умолчанию, когда для графики отводится все графическое окно. 2.9. Пример решения вычислительной задачи в Octave Пусть имеется набор точек на изображении. Требуется найти параметры прямой, наилучшим образом приближающей этот набор точек. Рассмотрим два критерия подбора параметров прямой. В качестве первого критерия выберем сумму квадратов отклонений y-координат точек от yкоординат прямой при совпадающих x-координатах. Для выбранного критерия оптимальные параметры прямой y  p1 x  p2 вычисляются минимизацией функционала: n F1   ( yi  p1 xi  p2 )2    min , p ,p i 1 1 2 где i, n – номера и количество рассматриваемых точек. Для демонстрации возможностей Octave решим задачу минимизации по 1му критерию тремя способами:  применением стандартной функции P = polyfit(x,y,N), вычисляющей методом наименьших квадратов вектор P – коэффициентов полинома степени N. При подборе прямой N = 1. Здесь x – вектор x-координат, а y – вектор y-координат набора точек;  приравниванием нулю двух частных производных от функционала F1 по параметрам прямой p1 и p2 и решения полученной системы 2-х уравнений с 2-мя неизвестными;  применением операции левого матричного деления: P=(X\y)’, где матрица X = [x' ones(n,1)]. Листинг программы, результаты вычисления коэффициентов прямой и её график представлены на рис. 6,7. Функция mnk_test1 является главной, и её имя совпадает с именем m-файла. Локальная функция generate_data вычисляет исходные данные по заданному вектору «p» коэффициентов модельной прямой 46 и среднеквадратическому отклонению «sigma» y-координат точек. Локальные функции estimate_line1, estimate_line2, estimate_line3, соответственно, 1-м, 2-м и 3-м способом оценивают параметры прямой. Вычисленные 3-мя способами параметры совпадают. Рис. 6. Набор точек и приближающая прямая (1-й критерий) Рис. 7. Листинг программы, результаты вычислений (1-й критерий) 47 Вернемся к постановке задачи. Предположим, что набор точек и, соответственно, оцениваемая прямая линия расположены на изображении почти вертикально. С увеличением угла наклона рассмотренный метод теряет точность (зеленый пологий график на рис. 8). Рис. 8. Набор точек и приближающие прямые (1,2-й критерии) Получить параметры прямой, наилучшим образом приближающей этот набор точек, можно, рассматривая вместо разницы ординат, расстояния от точек до прямой. Расстоянием является длина перпендикуляра, опущенного из точки на прямую. Пусть L – длина перпендикуляра, опущенного из начала координат на прямую, а θ – угол (измеренный в положительном направлении) между положительным направлением оси Ox и направлением перпендикуляра, опущенного из начала координат на прямую. Тогда уравнение прямой (нормальное уравнение прямой) имеет вид: x cos   y sin   L  0 . Введя обозначения: p1  cos  , p2  sin  , p0   L , запишем минимизируемый функционал в виде: n F2   ( p1 xi  p2 yi  p0 ) 2   min , при условии p12  p22  1 , p ,p ,p i 1 1 2 где p1 xi  p2 yi  p0 – длина перпендикуляра, опущенного из точки i на прямую, описываемую нормальным уравнением. Условная минимизация функционала может быть выполнена методом неопределенных множителей Лагранжа. В результате минимизации, вычисление коэффициентов p1 и p2 свелось к поиску компонент собственного вектора V=[p1, p2]T, соответствующего наименьшему собственному значению матрицы BTB, где B=[xT – x0 yT – y0], x – вектор-строка x-координат, а y – вектор-строка y-координат точек, x0, y0 – средние значения векторов x,y. 48 Отметим, что ненулевые собственные значения матрицы BTB и матрицы ковариации векторов x и y – BBT совпадают. Коэффициент p0 = – [x0 y0] ·V. График приближающей прямой, подобранной по второму критерию, представляет вертикальную линию (представлен на рис. 8 красным цветом) и согласуется с набором точек. Листинг программы и результаты вычисления коэффициентов прямой по 1-му и 2-му критериям представлены на рис. 9. Для вычисления собственных векторов и собственных значений использовалась стандартная функция eig, для вычисления средних значений векторов функция mean, а для поиска индекса наименьшего собственного числа функция min из пакета Octave. Рис. 9. Листинг программы, результаты вычислений (1,2-й критерии) Упражнение 16. Численно решить дифференциальное уравнение: dy  2 xy , y (0)  1 модифицированным методом Эйлера, методами Рунгеdx 49 Кутта, Кутта-Мерсона, Адамса, Милна. Решение сравнить с аналитическим 2 решением y ( x )  e  x и с решением, выполненным с помощью стандартной функции ode45 из пакета расширений Octave – odepkg. Загрузку пакета выполнить с помощью команды: pkg load odepkg. Построить графики решений. Полезной будет информация из главы 9 книги [8]. Упражнение 17. Численно решить дифференциальное уравнение (сведя его к системе дифференциальных уравнений 1-го порядка), моделирующее поведение осциллятора Ван Дер Поля: d 2x 2 dx   ( 1  x )  x  0, x (0)  2, x (0)  0 при   0, 1, 100, 1000 . dt dt 2 При   1000 система дифференциальных уравнений становится «жесткой». Убедиться, что стандартная функция Octave – ode45 не дает решения. Вместо функции ode45 применить функцию решения жестких систем – ode5r. Полезной будет информация из параграфа 5 главы 9 книги [8]. БИБЛИОГРАФИЧЕСКИЙ СПИСОК 1. Ющенко, Р. Что нужно знать про арифметику с плавающей запятой [Электронный ресурс] / Р. Ющенко. – Режим доступа: http://habrahabr.ru/post/112953/ (17.11.2013). 2. Брюс, Д. Comparing Floating Point Numbers, 2012 Edition [Электронный ресурс] / Д. Брюс. – Режим доступа: http://randomascii.wordpress.com/2012/02/25/comparing-floating-pointnumbers-2012-edition/ (17.11.2013). 3. Яшкардин, В. IEEE 754 – стандарт двоичной арифметики с плавающей точкой [Электронный ресурс] / В. Яшкардин. – Режим доступа: http://www.softelectro.ru/teoriy.html (17.11.2013). 4. Форсайт, Д. Машинные методы математических вычислений / Д. Форсайт, М. Малькольм, К. Моулер. – М. : Мир, 1980. – 280 с. 5. Goldberg, D. What Every Computer Scientist Should Know About Floating-Point Arithmetic / D. Goldberg // ACM Computing Surveys. – 1991. – V. 23, № 1. 6. Forsythe, G. How do you solve a quadratic equation? Technical report № CS40. June 16, 11966 . [Электронный ресурс] / Forsythe G. – Режим доступа: ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/66/40/CS-TR-66-40.pdf (17.11.2013) 7. Левитин, А.В. Алгоритмы: введение в разработку и анализ / А.В. Левитин. – М.: Вильямс, 2006. – 576 с. 8. Алексеев Е.Р. Введение в Octave для инженеров и математиков / Е.Р. Алексеев, О.В. Чеснокова. – М.: ALT Linux, 2012. – 368 c. 50 ОГЛАВЛЕНИЕ ВВЕДЕНИЕ .................................................................................................................. 3 1. МАШИННОЕ ПРЕДСТАВЛЕНИЕ ДЕЙСТВИТЕЛЬНЫХ ЧИСЕЛ .................. 4 1.1. Стандарт представления чисел с плавающей точкой .................................... 5 1.2. Специальные случаи представления чисел с плавающей точкой ............... 9 1.3. Представление денормализованных чисел .................................................. 10 1.4. Границы множества чисел с плавающей точкой ......................................... 11 1.5. Округление чисел с плавающей точкой........................................................ 13 1.6. Абсолютная и относительная погрешность округления............................. 13 1.7. Точность представления действительных чисел ......................................... 14 1.8. Машинный эпсилон ........................................................................................ 17 1.9. Сравнение чисел с плавающей точкой ......................................................... 17 1.10. Сравнение чисел с плавающей точкой с нулем ......................................... 19 1.11. Финансовые расчеты и числа с плавающей точкой .................................. 20 1.12. Машинное представление целых чисел ...................................................... 21 1.13. Неустойчивые алгоритмы ............................................................................ 23 1.14. Чувствительные задачи ................................................................................. 27 1.15. Сложность алгоритмов ................................................................................. 30 2. МАТЕМАТИЧЕСКИЕ ПАКЕТЫ ПРОГРАММ................................................. 35 2.1. Установка GNU Octave ................................................................................... 35 2.2. Введение в GNU Octave .................................................................................. 36 2.3. Представление данных в Octave .................................................................... 36 2.4. Функции формирования матриц .................................................................... 37 2.5. Сервисные функции ........................................................................................ 39 2.6. Функции матричных вычислений ................................................................. 39 2.7. Встроенный язык программирования ........................................................... 41 2.8. Графика в Octave ............................................................................................. 45 2.9. Пример решения вычислительной задачи в Octave..................................... 46 БИБЛИОГРАФИЧЕСКИЙ СПИСОК...................................................................... 50