Содержание к диссертации
Введение
1 Алгоритмика вычислений 15
1.1 Основные определения 15
1.11 Варианты организации интервальной стандартной про цедуры 28
1.1.2 Базовый промежуток 29
1.1.3 Достаточные признаки для выбора аппроксимации 29
1.2 Применение рациональных аппроксимаций 31
1.2.1 Меюдика учета инструментальной погрешности. Анализ погрешности вычисления положиіельного полинома. 31
1.2.2 Формирование локализующего интервала для значения функции на базовом промежутке 35
1.2.3 Выбор параметров, от которых зависит точность 37
1.2.4 Различные случаи задания коэффициентов вычисляемого полинома 38
1.3 Погрешность вычисления ряда 42
1.3.1 Введение 42
1.3.2 Оценка относительной величины арифметической погрешности 44
2 Исследование некоторых элементарных функций 47
2.1 Вариант процедуры для '-экспоненты 47
2.1.1 Общие соображения, выбор базового промежутка 47
2,1.2 Аппроксимация и локализующий промежуток для вы числяемой функции 48
2.2 Вариант процедуры для натурального логарифма 49
2.2.1 Общие соображения, выбор базовою промежутка 49
2.2.2 Границы для логарифма в базовом промежутке 50
2.2.3 Алгоритм нахождения и: сохранение монотонности 51
2 2.4 Поіреіпносіь нахождения и 52
2.3 Вариант процедуры для арктангенса 53
2.3.1 Общие соображения, выбор базового промежутка 53
2.3.2 Граница для арктангенса в основной части базового промежутка 54
2 3.3 Арифмегическам погрешность аргуменіа пененною ряда. 55
2.4 Сравнение предлагаемых методов с разрабоїанньши ранее GO
Реализация программных модулей. 65
3.1 Общие сведения 65
3.2 Общие программные модули 66
3,2.1 Библиоіека вспомогательных методов 66
3.3 Базовый класс вычислителя 68
3 3.1 Методы, реализующие общую для вгех функций логику вычисления 68
3 3.2 Меюды, нуждающиеся в переопределении 70
3.3.3 Методы, которые могут быть переопределены в случае
необходимости 70
3.4 Базовый класс вычислителя для рядов 71
3.5 Программная реализация вычисления экспоненты 72
3.5 1 Численные результаты работы программы для экспоненты 73
3.6 Программная реализация вычисления натурального логарифма 73
3.6.1 Численные результаты работы программы для натурального логарифма 76
3.7 Программная реализация вычисления арктангенса 78
3.7.1 Численные результаты работы программы для арктангенса 80
Заключение 82
А Текст программных модулей: основные классы 84
АЛ Базовый класс вычислителя 84
А.2 Базовый класс для вычисления рядов 99
А.З Класс для вычисления экспоненты 101
А.4 Класс для вычисления натурального лоїарифма 108
А.5 Класс для вычисления арктангенса 116
В Текст программных модулей: дополнительные классы 125
ВЛ Запуск вычислительного процесса 125
8.2 Реализация дополниіельньїх функций 133
8.3 Обработка исключительных ситуаций 139
Литература
- Достаточные признаки для выбора аппроксимации
- Выбор параметров, от которых зависит точность
- Вариант процедуры для натурального логарифма
- Методы, реализующие общую для вгех функций логику вычисления
Введение к работе
В насюящий момент акіивное развиїие вьічислиіельной техники и использование ее в научных исследованиях и различных областях производства требу-юг наличия алгоритмов, гарантирующих получение результата с точностью, не меньшей наперед заданной.
Однако, так как не все вещественные числа могут быть точно представлены машинными, чаще всего различные форматы представления сопоставляют множество битовых последовательностей определенной длины некоторому конечному подмножеству вещественных чисел.
О числах, входящих в такое подмножество, говорят, что они представимы в машинном формате. Числа, не представимые в машинном формате точно, могут быть выражены несколькими способами. Чаще всего используется выражение вещественного числа ближайшим числом, представимым в машинном формате, и предеіавление содержащим число интервалом, границы которого входят в указанное выше подмножество.
На практике ситуации, когда исходные данные известны точно, а расчеты носят сугубо целочисленный характер и проводятся по точным явным формулам, достаточно редки. Гораздо чаще 'і очный результат не представим машинными числами, поэтому получить его на компьютере не представляется возможным.
В таком случае традиционные методы вычислений предусматривают ответ на вопрос «чему приблизительно равно?». Однако точность полученного результата в большинстве случаев определить методами традиционных вычислений весьма проблематично в связи тем, что количество возникающих в процессе расчета погрешностей увеличивается вмесіе с количесівом вычислительных операций.
B 1985 году Институтом IEEE был принят сіандарт на вычисления с плавающей точкой [58]. В этом докуменіе указаны требования, предъявляемые к округлению чисел, выполняемому в операциях сложения, умножения, вычитания, деления и извлечения корня, что позволило сделать компьютерные вычисления более корректными и доказательными [64].
Однако указанный стандарт никак не регламентирует вычисление элементарных функций (таких как синус, косинус, логарифм, и т.д.). Более того, по той же причине в некоторых случаях такие функции могут вычисляться некорректно (иначе говоря, погрешность вычисления функции не гарантирована и может быть слишком велика в контексте какой-либо решаемой задачи).
Основной причиной для такого результата является Table Maker's Dilemma [77]: если есть число х и некоторая злеменіарная функция /, вычисление правильно округленного к формату числа с плавающей точкой результата у = /(.г) может потребовать слишком большой точности промежуточных вычислений (причем требуемая точность промежуточных вычислений заранее неизвестна).
Для решения этой проблемы было предложено множество алгоритмов и методов. Условно их можно разделить на четыре типа:
Традиционные вычисления.
Интервальные вычисления.
«Гибридные» вычисления.
Другие методы.
Рассмотрим перечисленные методы.
Традиционные вычисления Этот метод предусматривает использование классических способов вычисления, т.е. использование компьютера в качестве калькулятора, идентично копируя вычислительный процесс так, как он происходил без использования автоматизированных систем. Обработка ошибок либо не предусмотрена вообще, либо крайне ограничена. Однако точность резулыата никак не гарантируется.
Для решения проблемы точности используются методы, которые предусматривают манипуляции внешними по отношению к вычислительному процессу вещами - объемом оперативной памяти, увеличением размера переменных для хранения промежуточных и окончательных результатов (т.е. расширением множества вещественных чисел, представимых в машинном формате), уточнением начальных данных, и т.д. Однако, такой подход плох іем, чю он в большинстве случаев достаточно дорог, а кроме того, все равно не дает абсолютной гарантии точности полученных результатов, предоставляя пользователю полагаться лишь на то. что изменений, сделанных им, хватит, для того, чтобы обеспечить необходимую ему точность.
Сравнительно недавно в литературе стали появляться описания алгоритмов, которые умеют «подстраиваться» под требуемую точность Один из них был предложен в работе [57]. Он предусматривает увеличение ЮЧН0С1И вычислений в случае необходимости. Для этого в процессе вычислений проводятся тесты на точность, и если они не выполняются, точноегь вычисления увеличивается путем использования переменных большей длины. Однако критерии, используемые для сигнализации о необходимости увеличения точности, преде і являются не совсем надежными.
Интервальные вычисления. Другой метод поиска решения, так называемые ишервально - локализующие вычисления [48], ищут ответ на друї ой вопрос - «в чем содержится?», т.е. вмесю того, чтобы искагь результат вычисления в виде одного числа, оп ищется в виде вилки - ишерка-ла вида [а, 6], в котором гарантировано будет лежать точный результат. Этот способ как правило довольно дешев - в большинстве случаев он может использоваться па тех аппаратных и программных средствах, которые имеются в наличии (в отличие от традиционных методов, которые требуют минимум создания собственных типов данных), а кроме того, позволяет манипулировать точностью полученных результатов, и требует минимум программирования - в некоторых случаях іребуетея только
переопрсделение стандартных функций.
Термины «интервал», «интервальный» впервые появились в работе Те-руо Сунаги [83] в 1956-58-м, двумя годами раньше была опубликонана работа Мечислава Вармуса [85]. В этих работах была предложена классическая интервальная арифметика и намечены ее приложения. Кроме того, Т.Супага заложил основы интервального алгебраического формализма и дал весьма нетривиальные примеры применений новой техники, к примеру, в численном решении задачи Коши для обыкновенных дифференциальных уравнений [70].
Первая систематическая монография [76], посвященная локализующим вычислениям, была написана Раймоном Э. Муром 1966-м году.
В СССР у исюков развития этого направления еще в 20-х годах прошлого века стоял В.М. Брадис, предложивший использование «метода границ» ([6], [7]).
В 1962-м году в одном из первых выпусков «Сибирского математического журнала» появилась статья Лауреата Ленинской, Сталинской и Нобелевской премии Леонида Витальевича Канторовича [22], обозначившего эту тематику как приоритетную для нашей вычислительной науки.2
Дальнейшее развитие интервальные методы получили в СССР в рабоїах академика Николая Николаевича Япенко и ею ученика Юрия Ивановича Шокина, а также ряда других авторов ([11], [21],[50]). В иностранной литературе также было опубликовано множество работ, посвященных интервальному анализу, например [2]
В настоящий момент развиїие интервально- локализующих вычислений продолжается, и эти методы находят применение во все большем количестве отраслей науки и производства (медицинская техника, экономика [10], распознавание изображений, задачи классификации и регрессии (в том числе и с нечеткими множествами) [28], металлургия [9], задачи управления [8], и т.д.).
В рамках интервального подхода разработано большое количество мето-
дов, позволяющих получать различные интервальные расширения стандартных функций {например, через декомпозицию графика [43]).
Сущесівуют библиотеки интервальных функций для различных языков программирования ([72], [74]}, а также для специализированных вычислительных сред типа MATLAB ([25], [66]).
К недостаткам этого метода можно отнести слишком большое в некоторых случаях расширение локализующего интервала.
«Гибридные» вычисления. К этой группе огносигся вычислительные алгоритмы, использующие как традиционные, так и интервальные методы.
Например, А. Зив в 1991 году опубликовал работу «Past evaluation of elementary mathematical functions with correctly munded last bit» [86]. Методы, примененные им, были основаны на вычислении значения функции не на всей области определения, а лишь на некотором сравнительно узком базовом промежутке. При этом для вычисления требовалось использование так называемой «двойной» арифметики (т.е. величин, составленных из двух чисел типа double)
В настоящий момент резулыаты, полученные А. Зивом, активно используются для построения алгоритмов доказательных вычислений ([55], И).
Другие. Существует еще некоторое количество подходов к вычислению па компьютере. Например, метод, основанный на использовании таблиц заранее посчитанных констант. Одной из первых работ в этом направлении была [65]. Такие алгоритмы получили распространение в связи со значительным удешевлением оперативной памяти, требуемой ими (например, в алгоритме АТА-М в упомянутой работе используются таблицы размером 2Мб). Однако этот алгоритм предназначен для приближенного вычисления (хотя и достаточно быстрого), а не для получения гарантированных оценок.
Конечно, существуют и другие, менее распространенные методы.
Кроме общей идеологии процесса вычисления, на точность результата оказывают влияние методы получения мащинно-представимого числа из чочного (иначе говоря, методы округления).
Известно, что в случае получения в качестве результата интервала в явной (указаны границы) или неявной (указан результат, и изнргпю. чю точное решение отличается от него не более чем на с, где с может вычисляться, например, как A/i-s, где А - максимальное значение отношения абсолютной величины округления к единице младшего разряда, }Г$ - величина единицы младшего разряда переменной результата) форме наиболее часто встречаюі-ся три указанных далее способа.
Пусть значение какой-то переменной х содержится в некотором интервале [а, 6], а < Ь. Предположим, что надо оценить значение еЛ Очевидно, чїо оно содержится в интервале [е'1,^], но поскольку числа еа, еь могут оказаться не машинно-предславимыми, для оценки значения ех используют следующие интервалы:
Интервал, границы которого - приближенные значения чисел ta и еь. Для вычисления приближенных значений можно использовать, например, функции стандартных математических библиотек.
Интервал, полученный так же, как и первый, расширенный на некоторое заранее сосчитанное є (например, таким значением может быть А/Г'1').
Гарантированный интервал, т.е. интервал, нижняя граница которого - максимальное машинно-представимое число, не превосходящее с", а верхняя граница - минимальное машинпо-представимое число, не меньше еь
Необходимо отметить, что в первом случае истишюе значение (х может и не содержаться в том интервале, коюрым оно представляется. Однако в случае, когда вычисления носят итерационный характер, положительные погрешности одних операций с большой вероятностью компенсируются отрицательными погрешностями других.
Второй способ применяют при оценке значений выражений, являющихся суперпозицией нескольких функций. В этом случае часто вычисляют каждую функцию первым способом, а затем расширяют полученный интервал. При этом положительные и отрицательные ііогрешносіи также компенсируют друг друга, и вероятное гь того, что истинное значение выражения лежиі в полученном интервале, очень высока, хотя и не равна единице.
Треіий способ назьіваеіся вычислением с гараніированной точностью Несмотря на то, что при его использовании в итерационном процессе будет происходить постепенное расширение интервала, этот способ позволяет гарантировать, что истинное значение некоторого выражения лежит в интервале, которым это значение представлено.
Анализ доступных интервальных магматических библиотек ([26]) показал, что они производят вычисления либо первым (например, [63]), либо вторым способом (например, [71]).
Нередко хороших результатов удается достичь в построении алгоритмов для вычисления некоторого класса функций (например, в данной работе рассматриваются неубывающие элементарные трансцендентные функции). Так как у этого класса функций можно выделить специфичные свойства, есть возможность сузить результирующий интервал с помощью введения некоторых ограничений, про которые известно, что для данного класса функций они будут вьшолняться. В связи с іем, что на насюящий момент накоплен обширный материал как по аппроксимации функций, в том числе и трансцендентных ([3],[4],[5], [45],[46], [51]) так и по их свойствам ( [18], [23]), такой подход представляется досіаточно продуктивным.
Целью данной работы является создание вычислительных алгоритмов для получения максимально узкой оценки погрешности результаїа вычисления для заданного формата переменных, и интервала такой ширины, іаранти-рованно содержащею точный результат, без использования расширенных (удвоенных) типов данных для некоторого класса функций, представимых полиномом (или рядом) с неотрицательными коэффициентами. Иначе говоря, строятся методы, позволяющие получить как можно меньший по ширине
результирующий интервал в имеющемся формате переменных, для элементарных трансцендентных функций (например, экспоненциальная функции, логарифм, арктангенс). В связи с тем, что при построении алгоритмов учитываются различные свойсгва вычислительной системы (величина единицы младшего разряда в определенном формаїе, «точность» арифметики, и т.д.), они применимы на большинстве используемых сейчас вычислительных машин, и позволяют провести быструю адапіацию к новым.
Вычисления проводятся в некотором базовом промежутке [0, -г$], что, с одной стороны, требует дополнительных исследований вычисляемой функции, но с другой стороны позволяет повысить качество вычислений. Для рассматриваемых функций приводяїся критерии выбора базового промежуїка, способы приведения аргумента к базовому промежутку (входа в базовый промежуток) и получения значения функции для исходного аргумента по результату для аргуменіа в базовом интервале (выхода из базового промежутка).
Для проверки полученных алгоритмов написан программный комплекс, осуществляющий вычисления согласно приведенным методикам.
Диссертационная работа состоит из введения, трех глав, заключения и двух приложений, включает в себя семь иллюсіраций и четыре таблицы. Общий объем работы - 149 страниц.
Первая глава содержит общее описание шпоритмов проведения вычислений. Приводятся общие соображения о процедуре вычисления для положительного полинома, а также ряда, оцениваемого геометрической прогрессией.
Во второй главе рассматриваю і ся формулы для расчета погрешности вычисления экспоненциальной функции, натурального логарифма и арктангенса, показаны методы входа в базовый иромежуюк и выхода из него по завершении вычисления, приводятся криіерии для выбора базового промежуїка.
Третья глава посвящена описанию программной реализации приведенных алюритмов.
Приложение А содержит тексты основных программных модулей (классов, осуществляющих вычисления).
Приложение В содержит тексты дополнительных программных модулей
(классов, реализующих некоторые функции, необходимые в процессе вычисления).
Основными научными результатами, полученными и обоснованными is данной рабоїе, являются:
Сформулирована и доказана Теорема 1.2.1 об оценке относительной погрешности вычисления полинома (раздел 1.2 1). Получена формула для расчета погрешности (1.27), разработаны способы применения теоремы к различным случаям задания исходного полинома (раздел 1.2.4), в юм числе сформулирована и доказана Теорема 1.2.2 о случае приближенного задания коэффициентов вычисляемого полинома.
Сформулирована и доказана Теорема 1 3 1 об оценке методической погрешности вычисления ряда.
Разработан алгоритм расчета инструментальной и методической погрешности вычисления функций, нредставимых положительными рядами с неотрицательными коэффициентами (разделы 1 3 1-1.3.2, формула для расчета (1.52)).
Получена модификация алгоритма расчета погрешности для экспоненциальной функции (раздел 2.1).
Получена модификация алгоритма расчета погрешности для логарифмической функции (раздел 2.2).
Получена модификация алгоритма расчета погрешности для функции арктангенса (раздел 2.3).
К основным практическим результатам работы можно отнесги программную реализацию вычислительного комплекса, выполняющего расчет с гарантированной точностью экспоненциальной функции, логарифмической функции и функции арктангенса по полученным алгоритмам.
По результатам диссертации опубликовано 6 работ. Основное содержание составляют результаты работ [13], [14], [16].
Результаты диссертации докладывались автором на XXXV Конференции «Процессы управления и усюйчивосгь» (Санкт-Петербург, Петергоф, 14-16 апреля 2004 г.), XXXVI Конференции «Процессы управления и устойчивость» (Саіікт-Пеїербург, Пеіерюф, 11-14 апреля 2005 г.), XXXVII Конференции «Процессы управления и устойчивость» (Санкт-Пеїербург, Пеіер-гоф, 10-13 апреля 2006 г.), Всероссийском (с международным участием) совещании по интервальному анализу и его приложениям Интервал-06 (Санкт-Петербург, Петергоф, 1-4 июля 2006 г.).
Достаточные признаки для выбора аппроксимации
Для вычисления фансцендентных функций в ЭВМ применяются те или иные рациональные аппроксимации р{.г). В данной работе рассматривают! монотонные трансцендентные функции/(ж). Для определенности будем полагать их неубывающими на рассматриваемом множестве X.
Можно указать два крайних варианта интервальной процедуры вычисления функции.
Вариант 1 Задать функцию /(.г) приближенно или точно в виде композиции арифметических дейсівий и находить ее интервальное расширение как композицию интервальных расширений этих действий. Арифметические же процедуры составляют ту или иную применяемую интервальную программную сисіему. Этот вариант рассмотрен, например, в [41]. Необходимо отметить, что в нем не используеіся отдельная процедура мажоризации результата вычисления всей функции, так как мажориза-ция входит почти в каждую арифметическую операцию. Эшг вариант хорош в случае точного задания /, в этой ситуации ошадаюі забоїы о монотонности по включению в соответствии с п. 24 [41]. Однако, для успешного конструирования интервальной процедуры, на функцию накладываются некоторые ограничения ([75]}.
Вариант 2 предусматривает неинтервальное вычисление аппроксимации ip функции /(.г) на концах отрезка-операнда X = [г,т], получение таким образом отрезка Ш,ф(х)\ Ф(Х), затем учет полной погрешности такого вычисления, т.е. величины {х) -ip(x)},x = x,x.
Затем по результату этого учета проводится мажорирующее преобразование М{) отрезка Y = Ф[Х). Здесь через Ф обозначается реально вычисляемая функция Ф. Таким образом, процедура в этом варианіе представляет реализацию отображения F(X) М(Ф(Х)). коюрое, как и во Варианіе 1, являеіся расширением функции /. Здесь монотонность по включению не обеспечивается самим алгоритмом, и нужны специальные меры по ее обеспечению. В работе [41] было показано, что отображения Ф(Х) = [ р(т)М )] и Ф(х) моноюнны по включению, если соответственно функция \р или ф повторяет нестрогий харакіер моноюнносіи исходной стандартной функции /.
Базовый промежуток. Рассматриваемые функции удобнее вы числять не во всей обласіи их существования, а на некотором сравнительно узком базовом промежутке (БП) (подобный метод использовался, например, в [86], [55]). Он будет назначаться гіак, чюбы исследуемая функция/ на нем была моноюнной и неотрицательной. Таким образом, вычислению/ предше ствует пересчет данного г (т.е. конца отрезка X) в значение ,т, приведенное к базовому промежутку. После нахождения приближенного f(x) осуществ ляется его перевод в приближенное /(г). Базовый отрезок будет иметь вид Xs = [0, Хц] с тем ИЛИ ИНЫМ Хц. Пусть аппроксимирующая композиция р(х) функции f(x) развертывается в после довательность рабочих формул: X] — X Г 2 = р2[Х\) = ( 1.) (1-11) У = тн+] = yj(l+2{ri,...r„). Очевидна
Лемма 1.1.2 Функция tp(x) сохраняет на отрезке х гаракт.ер нестрогой монотонности функции /, если правые части рабочих формул: (.охраняют тот же характер монотонности по всем вхождениям аргумента х(= Xi) не убывают по промежуточным переменным Г\,...хп.
Из эюй леммы и гипотезы о сохранении моноюнносіи машинными арифметическими операциями следует Лемма 1.1.3 В условиях: Леммы 1 функция ф сохраняет характер нестрогой монотонности функции /.
Отсюда следует Теорема 1.1.2 В уеловиях Леммы 1 отображение ф(х) монотонно по включению. Тогда в соответствии со Второй теоремой о композициях [36] ин і ернальное расширение F(r) функции / будет также монотонно по включению. Рассмотрим полином с неотрицательными коэффициентами (р(х) = й0 + сцх + ... + апхп ф 0. (1.12) Предположим, что коэффициенты а, точно задаются машинными числами. Пусть ір{т) вычисляется по схеме Хорнера: р(х) = (..({апх + an-i)x + ап-2х + .. + сц)х + а0. Справедлива следующая
Теорема 1.2.1 При вычислении по схеме Хорпера относительная погрешность вычисления полинома (1.12) оценивается соотношением ф(х) - р{х) р{х) -(1-/?)2Ит)Г где обозначено 2п +1 п 2га - 1 _ , М = (1-/У)2»Д»Т + (1- )2»-2а»- + - + (Парамегр /? = А/Г5 был введен в (1.10)).
Выбор параметров, от которых зависит точность
Несмотря на то, что в идеальном случае {бесконечно точная арифметика) моноюн-ность выражения (2.6) очевидна, машинная монотонность данного выражения требует дополнительного исследования. Лемму 1 из [15] можно расширить на этот случай Лемма 2.2.1 При использовании во всех случаях расчета одинаковых режимов округления машинная монотонность выражения с одной операцией, потенциально проводимой с ошибкой, при условии, что результат не выходит за пределы наибольшего машинного числа, выполняется.
Доказательство Пусть a = f(x l,x"), Ь = /(.т , г ]). В результате работы машинной функции / возможны три случая:
1. Каждое из чисел а,Ь лежит в разных не соседних промежутках между соседними машинными числами. В этой ситуации диапазоны возможных значений округленных результатов точно не пересекутся, а значит, лемма для этого случая выполняется.
2. Числа a, b лежат в соседних промежутках между соседними машинными числами. Очевидно, результаты округления будут лежать относительно друг друга также, как и точные результаты, и разве лишь при том условии, что выполняется округление в сторону ближайшего машинного числа и числа я, Ь лежат близко к г„, округленные резулыаты будут совпадать.
3. Числа а,Ь лежат в одном промежутке между соседними машинными числами. В этом случае при всех режимах округления кроме режима округления в сторону ближайшего машинного числа результатом округления будет одно и то же число и для а, и для Ь. В случае же округления к ближайшему машинному числу, округленные значения а и b могут либо совпасть, если оба числа лежат близко к одной границе, либо находится па разных границах интервала соседних машинных чисел, но в любом случае, порядок их поменяться не может. Лемма доказана. Очевидно, что для случая (2.6) x lj2 = #1,2 - 1, х\ = X\ i + 1, и, следовательно, машинная монотонность выполняется.
Погрешность нахождения и. В случае произвольного основания машинной арифметики числитель и знаменатель могу г вычисляться с погрешностью. Так как вместе с делением присутствуют максимум три опера,-ции, приводящие к потере точности, в соответствии с (1.10)
Далее можно воспользоваться результатами, полученными в 1.3.2 для вычисления полной погрешности (очевидно, что при всех возможных значениях «, определяемых базовым интервалом, нет необходимости предусматривать дополнительные ограничения на значения аргумента для обеспечения непрерывности (1.44). При этом для о необходимо искать только верхнюю границу, т.к. нижняя в силу вхождения аргумента во все слагаемые будет равна 0. Вариант процедуры для арктангенса.
Общие соображения, выбор базового промежутка. За основу алюритмики примем предсіавление рядом, сходящимся на всей числовой оси, из [19] (формула А 641.3): х +2? (2fc)! / J2 \к ,п , arctg = 7ТТ5 №ТТ) \ТЇ ) (216)
Существует еще несколько формул для аппроксимации арктангенса (например (55.4) из [20], также есть аппроксимации в [45]), однако в приведенной выше формуле влияние инструментальной погрешности существенно больше, чем Б остальных, в силу порядка чисел, используемых в ней (22\ факториалы), посему для демонстрации работоспособности метода избрана именно она.
Члены правой части одинаково монотонны для всех вхождений аргумента т. Можно переписать выражение (2.16) в виде: arctg(x) = ju(i)v(u), (2.17) где обозначено v{u) = аки\ (2.18) к о и = и(т) = г --2. (2.19) Эта функция возрастает, а с ней и члены правой части (2.1G) при х 0. Что же касается случая .г 0, то по нечетности арктангенса он сводится к противоположному.
В оценке скорости сходимости ряда по Даламберу участвует отношение соседних членов. Найдем ею асимптотику при к — +оо. Согласно (2.16) ашхш _ 2fc + 112fc+l х2 х2 акхк fc + 1 22& + 31+я2 1 + z2 Этот предел при всех х меньше единицы, но при больших близок к ней, и потому ряд буде г сходиться очень медленно.
Вариант процедуры для натурального логарифма
Это недостаток был в преодолен, например, в методах, основанных на методе разложения Бернштейна (Bernstein Expansion), подробно описанных в работах [61], [62], [73] [80]-[82]. Однако приведенные там алгоритмы по сравнению с полученными в настоящей работе достаточно сложны и требуют большою количества преобразований (например, для получения полиномов Бернштейна), активно используют операции над матрицами, что увеличивает время рабоїьі и затраты ресурсов вычислительной системы. При этом, безусловно, при разовых вычислениях, проводимых в специализированных системах, они остаются актуальными.
При сравнении предложенных в работе алгоритмов с интервальными алгоритмами надо учитывать тот факт, что использование интервальной ариф-меіики требует наличия специальных библиотек или создания собственных библиотек для работы с типом данных «интервал (фактически сейчас суще-сівуеі только один компилятор, имеющий встроенную интервальную арифметику - Fortran-95 от Sun Microsystems).
Кроме того, существует довольно большое количество алгоритмов расчета (в том числе и трансцендентных функций), реализованных в различных библиотеках. В [52] приводя і ся ссылки на 227 работ в этой области. Однако, как отмечалось, например, в [26], основным недосіатком опубликованных в Интернете библиотек, находящихся в свободном доступе, является отсутствие гарантий со стороны авторов коррекіносіи их работы. Например, библиотека fdlibm [71] не строится по причине отсутствия каких-то функций (то же самое можно сказать и про библиотеку filibt \ [72]), в C-XSC [87] не работают поставляемые в комплекте примеры. Другой проблемой, с которой приходится сталкиваться при проведении сравнения, является тот факт, что некоторые разработки (например, INTLAB) основаны на типах данных с нестандартным размером. В таких случаях выводы о свойствах вычислительной системы можно сделать юлько на основе анализа кода.
Рассмотрим некоторые я І имеющихся свободно распространяемых библиотек. Например, в известной интервальной библиотеке INTLAB [66] реализован метод расчета полиномов и трансцендентных функций.
Анализ докуменіации к INTLAB показывает, что основным отличием предлагаемых в наеюящеЙ работе мєіодов от используемых в INTLAB является юг факт,чю в INTLAB для поиска округленных значений используется таблица заранее посчитанных значений, комруго при реализации используемых алгоритмов на другой вычислительной системе придется, соответственно, і OJ овить заново перед реализацией вычислительною процесса В предложенных в настоящей работе методах такой необходимости нету, т.к. максимально возможная погрешность вычисления определяема автоматически исходя из разложения результата на составляющие в соответствии с формой представления числа в компьютере, при этом качество арифметики является одним из параметров, что позволяет быстро адаптировать их к новым вычислительным системам.
Полиномы в INTLAB вычисляются также по схеме Хорнера (хотя присутствуют и другие, менее точные в соответствии с [54]), причем учет погрешности проводится в интервальной манере посредством расширения результата до ближайших машинно-представимых чисел. Вычисления экспоненциальной функции, натурального логарифма и арктаніенса проводятся посредством стандартных функций MATLAB. Описания его внутренних ал-горишов работы в открытых источниках, к сожалению, найти не удалось, но эксперименты показьіваюі , что, скорее всею, внутренние вычисления в MATLAB проводятся с точностью большей, чем возможно вывести на экран, и лишь потом округляются к нужному формату результата. Исходя из вышесказанного, можно сделать вывод, что точиосль вычисления в INTLAB достигается в том числе за счет использования нестандартных размеров переменных и операций с чакими переменными. Явного вычисления погрешности не используеіся.
Кроме того, нужно 01 метить, что на качество вычислений в INTLAB ока.-зывают влияние ошибки, допущенные производителями MATLAB. Часть из них (но, скорее всего, не все) известна, судя по многочисленным комментариям типа «Matlab bug fix» в коде INTLAB. Рассмотрим другую известную интервальную библиотеку filib+ч- [72].
Алгоритмы эюй библиотеки также отличаются возможностью использования только в среде разработки, для которой они писались - адаптация их довольно сложна. Связано это не юлько и не столько с механизмами смены режимов округления (что необходимо для корректной работы библиотеки), выполняемого native-функциями и не всегда доступного в других программных средах, но и обширным использованием так называемых table lookup methods, которые, во-первых, требуют в общем случае большого количества памяти для хранения таблиц (см., например, [79]), а во-вторых, требуют наличия самих этих таблиц, точностью которых ограничена точность вычисления. Кроме того, возможности учета качества арифметики произвольной вычислительной системы в библиотеке filib-l- І отсутствуют.
Вычисление логарифма, экспоненты и арктангенса в библиотеке filib++ происходи і с использованием table-lookup methods
Еще одной распространенной библиотекой для гарантированного вычисления интервальными методами является созданная в Swiss Federal Institute of Technology {Lausanne, Switzerland) и разрабатываемая ныне в Laboiatoiie d Informatique de Nantes-Atlantique (France) библиотека интервальных функций GAOL [59] (в настоящий момент на сайте библиотеки доступна версия 2.0.2). Она позволяет проводить вычисления в интервальной манере и том числе полиномов и стандартных чрансценденіньгх функций.
Методы, реализующие общую для вгех функций логику вычисления
В базовом классе вычислителя BasicSeriesCalculator реализован общий алгоритм вычислений для рядов в соответствии с 1 3.1. Для корректной работы при вычислении различных функций необходимо создать наследующий класс и переопределить в нем некоторые функции в дополнение к тем. которые были переопределены для BasicCalculator. Класс содержит 3 метода. 1. public double calculatcSigmaO Вычисление а в соответствии с разделом 1.3.2. Метод должен быть переопределен наследующим классом в -72 соответствии со спецификой вычисляемой функции. Все необходимые данные берет из свойств класса. Возвращает а. 2. public double calculateSigmal О Вычисление а\ в сооївеїсівии с раз делом 1.3.2. Me і ода должен быть переопределен наследующим классом в соответствии со спецификой вычисляемой функции Все необходимые данные берет из свойсів класса. Возвращает crl. 3. protected double calculateXi_ 0 (double х) Вычисляет о по формуле в соответствии с разделом 1.3.2. Переопределять не следует, 3.5 Программная реализация вычисления экспоненты.
Вычислитель экспоненты реализован в классе Exponential на основе класса BasicCalculator. Для реализации были переопределены следующие методы класса BasicCalculator: 1. protected void enterBaszcInteryalAndCalculatelntermediateAgrumentidovble x) Вычисление осуществляется по следующему алгоритму: (a) Аргумент проверяется на неотрицательность. В случае, если аргумент отрицательный, далее рассматривается ею абсолютная величина и взводится флажок negativeargument. (b) Производится деление аргумента на base (= }ї) до і ex пор, пока он не будет меньше basicintervallo (в случае экспоненты 1). (c) Результат предыдущего шага записывается в поле и, степень основания системы счисления, необходимая для приведения к базовому промежуїку - в поле класса s. 2. protected Reallnterval exitBasirInterval(Reallnterval res) Выход из базового промежутка осуществляется по следующему алгоритму: -73 (a) В случае s 1 интервал res возводятся в степень base\ (b) Если установлен флажок negativeargument, выполняется вычисление 1/res. 3. protected Reallnterval calculateForPreciseKnownPomtO Возвраща ет значение en = 1. 4. protected double calculatePhi (double gamma) Вычисляет значение ап проксимирующей функции согласно (2.5). 5 double getGammal О Возвращает 1 (согласно (2.4)). 6. double getGamma2() Возвращает 3 (согласно (2.4)). 7. protected double ca,lculateNahveValue(double x) Методы вызывает стандартную функцию Java для экспоненты и возвращает полученный с ее помощью результат.
Численные результаты работы программы для экспоненты. Таблица 3.1 содержит результаты работы функции расчета экспоненты. Здесь х - аргумент, n, F(x) - нижняя граница результата, F(x) - верхняя граница результата, f(x)java - результат работы стандартной функции Java для расчета экспоненты.
Качесіво рабоїьі функции расчета экспоненш приведено на Рис. 3.1. На графике представлена зависимость ширины результирующего интервала от аргумента.
Как видно, интервал расширяется с увеличением аргумента, что вполне обьяснимо растущей погрешносіью при выходе из базового промежутка
Для реализации были переопределены следующие методы класса BasicCalculator. 1. protected void enterBasicIntervalAndCakulateIntennediateAqrum,ent(double x)
Вычисление осуществляется согласно 2.2.1 по следующему алгоритму: Определяются маніисса т и порядок р аргумента Если значение машиссы m y/Jl, взводится флажок reverseargument и в качестве аргумента для вычисления и используемся частное от деления основания внутренней системы счисления на мантиссу аргумента
Если значение мантиссы т y/JJ, оно используется в качестве аргумента для вычисления и (ситуация и 1 невозможна в соответствии с требованиями Стандарта [58]). (d) Вычисляется промежуточный аргумент согласно (2.6). 2. protected Reallnterval еіпШ(шс7пегуа/(Real Interval res) Выход из базового промежутка производится по одной из формул: (2 8) или (2.2.1} в зависимости от состояния флага reverseargument. При этом слагаемые также могут быть выведены отдельно, что позволяет во многих случаях сузить результирующий интервал. 3. protected String checkAгgument(double x) Проверяет положитель ность аргумента. 4. protected Reallnterval calculateForPreaseKnownPomtO Возвраща ет значение ln(l) = О 5. protected double calculate Phi (double gamma) Вьггаоляеі значение ап проксимирующей функции согласно (2.13). In + 1 6. double getGammalO Возвращает (согласно (2.11)). 7. double getGamma2() Возвращает 1 (согласно (2.12)). 8. protected double calculateAratwe Value (double x) Методы вызывает сіандартную функцию Java для натурального логарифма и возвращает полученный с ее помощью результат. Кроме того, были переопределены следующие методы класса Basic Calculator: 1. public double calculateSigmaQ Вычисление а в соответствии с (1.51). 2. public double calculateSigmal О Вычисление rl в соответствии с разделом (1.49).
