Метод оценки погрешности округлений значений вычисляемой функции, основанный на варьировании длины мантиссы в арифметике с плавающей запятой тема диссертации и автореферата по ВАК РФ 01.01.07, кандидат физико-математических наук Гриневич, Алексей Иванович

• Гриневич, Алексей Иванович
• кандидат физико-математических наук
• 2013

Оглавление диссертации кандидат физико-математических наук Гриневич, Алексей Иванович

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

ГЛАВА 1. Математические модели машинных вычислений

1.1 Формат машинного числа в стандарте IEEE 754

1.2 Погрешность базисных математических операций

1.3 Переполнение и потеря значимости

1.4 Библиотеки программ использующие стандарт IEEE 754

1.4.1 GNUGMP

1.4.2 GNUMPFR

1.5 Другие подходы

1.5.1 Точная модель рациональных чисел

1.5.2 Арифметика с многоуровневой экспонентой

ГЛАВА 2. Элементарный алгоритм и его свойства

2.1 Основные определения

2.2 Погрешности округлений решений

2.3 О заданной точности решений

2.4 Выводы

ГЛАВА 3. Метод К-решевий оценки погрешностей округления

3.1 К-решения задачи и их свойства

3.2 Правила оценки погрешностей округления решений. Конечношаговый алгоритм (КША)

3.3 Оценка погрешности округления решения по совпадению первых десятичных знаков

3.4 Правило оценки точности решений. Бесконечношаговый

алгоритм

3.5 Об эффективности метода К-решений

3.6 Выводы

ГЛАВА 4. Численный эксперимент

4.1 Вычисление полинома

4.2 Решение системы линейных уравнений

4.2.1 Решение СЛУ методом Гаусса

4.2.2 Решение СЛУ методом сопряженных градиентов

4.3 Задача нахождения производной 1го порядка

4.4 Нахождение собственных чисел матрицы методом Данилевского

4.5 Задачи оптимизации

4.5.1 О погрешностях округления метода Ньютона

4.5.2 Тестовые задачи БМ

4.5.3 Тестовые СНУ

4.5.4 Параметры численного эксперимента

4.5.5 Правила варьирования мантисс

ЗАКЛЮЧЕНИЕ

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

ПРИЛОЖЕНИЕ 1. Метод штрафных функций

ПРИЛОЖЕНИЕ 2: Сеточные методы

Задача Коши

Уравнение теплопроводности

Сокращения

ВМ вычислительная математика

МЧ машинное число (число в арифметике с плавающей

запятой)

ЗТР заданная точность решения

ЭАВМ элементарный алгоритм вычислительной

математики

ИППМ итерационная последовательность с переменной

мантиссой

КУ коэффициент уменьшения (#)

g -устойчивость характеристика ИППМ

К-решение контрольное решение в ИППМ

ПВМ правило варьирования длины мантиссы

ПО погрешность округления

СЛУ система линейных уравнений

СНУ система нелинейных уравнений

Введение диссертации (часть автореферата) на тему «Метод оценки погрешности округлений значений вычисляемой функции, основанный на варьировании длины мантиссы в арифметике с плавающей запятой»

ВВЕДЕНИЕ

Повышение эффективности численных методов (ЧМ) решения задач вычислительной математики (ВМ) - проблема стоящая перед специалистами от самых начальных времен возникновения ВМ по настоящее время и останется актуальной в обозримом будущем, т.к. в прикладной математике постоянно возникают новые задачи, усложняются известные и растут требования к точности и надежности численных методов их решения. Существуют три основных фактора развития, которые повышают эффективность ЧМ:

]. Разработка новых и совершенствование известных методов решения задач ВМ, включая теоретический анализ их свойств сходимости, порядка аппроксимации, устойчивости и т.д.

2. Создание новых ЭВМ с лучшими техническими характеристиками: количеством оперативной и постоянной памяти, тактовой частотой и количеством процессоров.

3. Разработка программного обеспечения с новыми возможностями.

Для дальнейшего изложения высказанные выше общие положения требуется уточнить:

Численный метод решения некоторого класса задач ВМ будем считать эффективным, если при его использовании на данной ЭВМ решение задачи с заданной точностью находится за приемлемое время.

Однако, это определение достаточно расплывчато и оно также нуждается в уточнении. Чаще рассматривают сравнительную эффективность ЧМ, т.е.:

Данный метод эффективнее (предпочтительнее) по времени решения чем другой, если он решает задачу с заданной точностью за меньшее время.

Другой характеристикой эффективности численных методов является точность решения задач ВМ. В настоящей работе основное внимание будет уделено именно анализу погрешностей округления и достижению ЗТР.

Как будет определено ниже (Глава 2) решение некоторой задачи ВМ сводится к нахождению приближенного значения вектор-функции /(х)е Rk в точке хе R", где f(x)eRk - точное решение задачи. В дальнейшем вычисленное значение функции и найденное решение задачи мы будем восприниать как синонимы. Далее написание аргументов хе R" функций fix), /(*), там, где это возможно, будем опускать. Решением задачи ВМ

называется точка / е Rk из множества:

М = {/ е Rk : ф(7,/)< 4 где ф(7,/) - 0 1

неотрицательная функция определяющая заданную точность решения задачи при выполнении неравенства Ф(7,/)^£; е >0 - заданный параметр точности решения или точность решения. В процессе решения задачи неравенство Ф(7,/)^е является условием окончания метода решения. Функции ф(7, f) бывают различными и определяются методом решения. Например, Ф(-, ) - функция нескольких решений f: ф(7,...,/,,/) или ф(7,/)=ф,(7, ,7/) _ т-е- Ф](у) не зависит от значения точного решения. Приведем пример.

Пусть требуется решить невырожденную систему линейных уравнения (СЛУ) Az = c, z,ceRk; здесь f = z, х = (л,с). Пусть 01(z) = ||y4z-c||;

здесь и далее ||*| - евклидова норма вектора ие Rk: ||и|| = J^".2 • Множество

решений м удовлетворяющих условию заданной точности решения СЛУ очень часто представляется в виде: М - {?е В.к: ¡ЛЗТ-сЦ <

Известны другие оценки для погрешности решения СЛУ [58]. Например, если считать значение матрицы точным, то имеем оценку:

02

где достижение заданной точности для Ф{г,г)=\\гвозможно на практике только через оценку значения функции Ф,^)^ ЦуГ'Ц-Цл?-с||, т.е.

Ф(2>) <«&,(*)<£ 03

или \\1 -гЦ^Цл^'Ц-Цл?-с\<е, т.к. реально возможно вычислить только значение Ф,(г).

Таким образом, при численном решении задач ВМ желательно иметь оценки прямой погрешности решения:

0.4

и такие технические возможности вычислительного процесса, которые позволяют находить значения функции /(*) для заданных значений £>о, порядок которых г, где г = 1(Гг; г- натуральное число может изменяться в заданном диапазоне значений. Погрешность решения задачи (0.4) в общем случае может зависеть от параметров метода решения и длины т мантиссы МЧ, от которой зависит величина погрешностей округления [3],[16],[12]; значение д->о, если все составляющие погрешности ^ о.

Введем определение, уточняющее выше сказанное:

Определение 1: Пусть задано е > о - требуемая точность значения /, т.е. если точность решения достижима, то имеет место неравенство

д = ||/-/||<£ или Д/|/|<£г. Будем говорить, что имеет место заданная

точность решения (ЗТР), если задача решается на ЭВМ методом, для которого известны оценки Ф, (/],..,/г) погрешностей решения, обеспечивающие достижение указанной точности; значения оценок погрешностей определяются вместе с искомым решением и для них выполнено указанное условие достижимости точности (0.4). □

Данное определение ЗТР не учитывает специфики погрешностей вычисления функции, а потому его можно применять как для каждого вида погрешностей, так и для любой их комбинации. Настоящее исследование посвящено в первую очередь изучению погрешностей округления вычисленных значений функции при изменении длины т мантиссы МЧ.

В настоящее время большое внимание научного сообщества уделяется вопросам «вычислений с высокой точностью», т.е. таким вычислениям, при которых имеется возможность назначить необходимое (требуемое) значение длины мантиссы МЧ. Как отмечено в обзоре 2012 года [60], можно выделить целый ряд направлений исследований, где стандартной арифметики оказывается недостаточно. Среди них есть как давно известные проблемы, так и относительно новые, активное изучение которых началось вместе с появлением достаточных вычислительных мощностей. Это решение определенных видов ОДУ [65], вычисление рекуррентных соотношений [66], определение экспоненциально малых явлений в динамических системах [67], компьютерное исследование новых математических соотношений [68], моделирование сверхновых звезд [69]. В частности, А. Фролов утверждает, что «сейчас мы научились рассматривать и решать задачу нескольких тел с ограничениями, о чем нельзя было и мечтать всего несколько лет назад». Для решения задачи им использовалась арифметика высокой точности (120 знаков) для нахождения собственных векторов почти вырожденных матриц размерами 5000x5000 в рамках решения задачи п тел [71,72].

Как правило, вычислительные сложности возникают, когда решаемая проблема содержит одну из следующих вычислительных задач.

1. Задача решения плохо обусловленных СЛУ [2,3,16].

2. Распараллеливание вычислений. Задачи, разрешимые на одном процессоре приобретают новую степень сложности при распараллеливании. В частности, вычисление суммы ряда в условиях параллелизма, когда точный порядок суммирования неконтролируем, приводит к появлению дополнительных вычислительных погрешностей [62].

3. Моделирование долговременных физических процессов. Как правило, в условиях большого количества итераций возникают искажения, вносимые накоплением ошибок округления промежуточных вычислений [63].

4. Проблемы экспериментальной математики. Такие, как вычисление большого количества знаков в числах п и е для установления возможной «нормальности» [64], проверка гипотез о математических соотношениях и т.п.

Актуальность работы. Время вычислений существенно (рост времени выполнения на 1-2 порядка) возрастает при переходе на «высокоточную» арифметику с варьируемой длиной мантиссы. Поэтому подобные вычислительные проблемы считались слишком трудоемкими до недавнего времени, и ученые старались изыскивать специфические методы решения для конкретных задач или упростить задачу до «разрешимого» варианта. В последние годы ситуация начала меняться с появлением библиотек программ с интерфейсами высокого уровня. Здесь следует отметить и специальные подпрограммы для нахождения значений функций (АКРКЕС, вМР, МБРЯ, МРРЯ-Н-, МРРШ90, дБ), пакеты компьютерной

алгебры, позволяющие находить решение символически без потери точности (Maple, MathCad, Mathematica) и реализации языков программирования, позволяющие задавать длину мантиссы (LISP, Python, Perl, Haskell, Ruby). При переходе на высокоточную арифметику, как правило, оказывается возможным не переводить программу целиком, а заменять лишь часть ключевых алгоритмов на более «точные» варианты. Это позволяет локализовать вычисления, требующие высокой точности и минимизировать влияние возрастающего времени выполнения до приемлемых значений. Таким образом, безусловно, вычисления в арифметике высокой точности не могут рассматриваться отдельно от других подходов по оптимизации вычислений.

Отметим исследования, в которых получены оценки погрешности решений СЛУ и других задач линейной алгебры (JIA), зависящих от длины мантиссы МЧ, позволяющие получать ЗТР при заданном значении точности £>о [5], [16], [3], [2] и др. Однако использование теоретических оценок погрешности решений СЛУ и других задач ЛА в вычислительной практике затруднено из-за трудоемкости определения числа обусловленности матриц в методах, рассмотренных в [5], [16], [3], и из-за громоздкости оценок в [2].

Большой цикл работ посвящен анализу погрешностей округления в рамках интервального анализа (ИА) [41,42,44,46]. В этих работах оценивается величина интервала, которому принадлежит значение погрешности округления на всех шагах алгоритма решения задачи ВМ. Интервальная арифметика являет собой принципиально иной подход к анализу погрешностей вычислений. Интервальная арифметика оперирует не числами, а интервалами возможных значений, заключающими в себе искомое решение. Отметим, что по мере роста количества вычислений в композиции, трудоемкость оценки интервала возрастает, причем относится к классу NP-трудных задач. В реальности это означает, что при работе с

итерационными вычислениями и задачами большой размерности построение интервальной оценки становится затруднительным за приемлемое время.

В монографии [43] рассмотрено применение методов интервального анализа для решения задач глобальной оптимизации, в которой дано описание программного обеспечения для решения указанных задач.

Существуют библиотеки, реализующие интервальные вычисления [13]. Также ведётся разработка стандарта для интервальной арифметики рабочей группой IEEE [14]. Особенностью оценок погрешностей округления в рамках ИА является то, что они строятся для отдельных классов задач ВМ и для получения численных значений необходимо иметь специальное программное обеспечение [40].

Известный алгоритм оценки погрешностей вычислений с автоматической коррекцией ошибок округления первого порядка - метод CENA применяется для узкого класса функций [19] и в широкой вычислительной практике не нашел применения.

Перспективны возможные алгоритмы учета вычислительных погрешностей на основе специальных моделей представления машинного числа [25], [40]. Практическим ограничением таких алгоритмов на настоящий момент является недостаточная база библиотек, обеспечивающих реализацию стандартных вычислительных процедур и функций на основании этих представлений чисел.

Уже указанные работы по анализу методов оценок погрешностей округления, которые позволяют находить значения функции с заданной точностью, говорят об актуальности поиска таких научно-технических решений. В рассмотренных работах длина мантиссы рассматривается либо как ограничение, либо как как параметр, длина которого определяется заранее на основании известных оценок. Но современные технические

средства позволяют производить вычисления, используя мантиссу как динамически изменяемый параметр. Поэтому приобретает актуальность вопрос построения методов оценивания погрешностей округлений в условиях, когда длины мантиссы является варьируемым параметром. Желательно, чтобы методы анализа погрешностей округления были применимы для возможно более широкого класса задач. По нашему мнению построение методов (методики) оценок погрешности решений задач ВМ обладающих выше указанными свойствами является актуальной научной проблемой.

Приведем некоторые данные о возможностях современного программного обеспечения, на основе использования которого будут строиться методы оценивания погрешностей округления рашений задач ВМ. В настоящее время в бесплатном доступе получила распространение библиотека программ GNU GMP [9], реализующая стандарт IEEE 754 [1,38], в которой длина мантиссы в арифметике с плавающей запятой варьируется в широком диапазоне значений. Библиотека GNU GMP позволяет оперировать числами с длиной мантиссы от mmin = 24 вплоть до mmax = 231 = 2 147 483 648 двоичных знаков, чему соответствует 8 и 646 456 993 десятичных знаков. Верхнее значение длины мантиссы МЧ т^ невообразимо огромно. В указанной библиотеке также реализована возможность динамического изменения длины мантиссы m в различных сегментах программы от т^,, до т,^. Появление в свободном доступе программного обеспечения с такими возможностями расширяет границы для получения решений широкого круга задач ВМ с гарантированной точностью высокого порядка. Гигантское значение мантиссы т^ в библиотеке GNU GMP в широкой вычислительной практике вряд ли когда будет востребовано. Поэтому были созданы пакеты математических программ с меньшим значением т^. Например в пакете Maple

ramax = 65535 десятичных знаков, которое также намного больше значения мантиссы для четверной точности.

Библиотека GNU GMP стала основой реализации вычислений во многих математических программах, языках программирования и библиотеках.

Исторически сложились два основных подхода к уменьшению погрешности округлений решений задач ВМ:

1. Разработка алгоритмов, учитывающих специфику решаемой задачи, особенности машинной арифметики,

2. Повышение точности МЧ, т.е. увеличение длины мантиссы.

До настоящего времени основной упор делался на первый подход. Это объясняется тем, что возможности второго подхода исторически были ограничены возможностями аппаратной части и ограничениями оперативной памяти.

Можно сказать, что в последнее время основные ограничения на 2-й подход сняты, т.к. оперативная память и доступные вычислительные ресурсы достаточны для решения широкого класса задач ВМ.

Подход с варьированием мантиссы для получения более точного решения задач ВМ известен давно. До недавнего времени длина мантиссы машинного числа рассматривалась как параметр с 2-3 возможными значениями (одинарная, двойная, возможно, четверная точность). Кроме того, известно, что ошибка вычислений, вносимая погрешностями огругления, не является монотонной функцией длины мантиссы [47].

Цели диссертационной работы

1. Построение алгоритма вычисления значения функции /С*), применимого для достаточно широкого класса задач ВМ. Для предложенного алгоритма - исследование оценок погрешностей

округления значений /(х), зависящих от длины т мантиссы машинного числа (МЧ) и получение такого значения т0, которое обеспечивает достижение требуемой точности решения.

2. Разработка численных методов оценки погрешностей округления значений /(*), гарантирующих достижение заданной точности (ЗТР), и исследование их свойств.

3. Численное исследование предложенных методов оценок погрешностей округления для некоторых классов задач ВМ.

Научная новизна

1. Предложен численный алгоритм решения задачи ВМ, названный «Элементарный алгоритм вычислительной математики» (ЭАВМ), представляющий собой естественное объединение стандартных (базовых) вычислительных операций из какой-либо библиотеки программ. ЭАВМ обладает тем свойством, что при его реализации в точной арифметике будет получено точное значение вычисляемой функции.

2. Для ЭАВМ с конечным и бесконечным числом шагов (КША и БША) получены теоретические оценки погрешностей значений функции, зависящие от её аргументов и длины мантиссы т МЧ, гарантирующие достижение заданной точности решений.

3. Предложен метод К-решений (КР) - численной оценки погрешностей округления значения вычисляемой функции. Введены понятия итерационной последовательности значений функции с переменной мантиссой (ИППМ) и её £-устойчивости. Для g-устойчивой ИППМ доказана возможность достижения заданной точности решения.

4. Предложены алгоритмы, позволяющие оптимизировать процесс достижения заданной точности решения.

Практическая ценность. 1. Предложенный метод оценок погрешностей округления применим для достаточно широкого класса

14

вычисляемых функций. 2. Его использование не требует изменения используемых вычислительных алгоритмов. Для получения значений погрешностей необходимо вычислить значения функции / (х) при разных

длинах мантисс mi МЧ и сравнить их. 3. Предложен вариант метода варьирования длины мантиссы для получения заданной точности решения при решении задачи безусловной минимизации методом Ньютона.

Апробация. Основные положения диссертационной работы докладывались, обсуждались и получили одобрение специалистов на научных семинарах кафедр высшей математики и математических основ управления Московского физико-технического института (государственного университета) (2007-2013), научном семинаре отдела прикладных проблем оптимизации в Вычислительном центре им. A.A. Дородницына РАН (2013).

Публикации. По теме диссертации опубликованы 6 печатных работ из них 4 [48][49][52][53] из списка, рекомендованного ВАК РФ.

Личный вклад автора. Как содержание диссертации, так и основные положения, выносимые на защиту, отражают личный вклад автора в опубликованные по теме диссертации работы. Все представленные в диссертации результаты получены лично автором.

Структура. Диссертация состоит из введения, четырех глав, заключения и приложений. Список использованной литературы содержит 72 наименования.

Содержание работы.

Заключение диссертации по теме «Вычислительная математика», Гриневич, Алексей Иванович

3.6 Выводы

Результаты полученные в главе сводятся к следующим положениям:

1. Введено понятие КР задачи, которое в оценках погрешностей решений с некоторой точностью заменяет истинное решение задачи f(x). Исследованы свойства g-устойчивости последовательности решений, в том числе доказана теорема о достижимости решения с требуемой гарантированной точностью и получены оценки погрешности, далее численно реализуемые в ИППМ.

2. Предложены алгоритмы, позволяющие оптимизировать процесс решения задачи в ИППМ. Рассмотрены методы округления полученных решений, причем округленное решение обеспечивает заданную точность.

3. Предложен метод (правило) округления решения по совпадению t первых десятичных знаков (СПЗ), погрешность которого не превышает значения е = юе', где е-порядок числа.

4. Предложенный метод КР оценки погрешностей округления обладает следующими свойствами: a. В g-устойчивой ИППМ обеспечивается ЗТР, b. Метод КР обладает универсальностью в том смысле, что он не ориентирован на конкретные алгоритмы решений задач.

4 ГЛАВА 4. Численный эксперимент

Для проведения численного эксперимента использовался математический пакет Maple 13 [7], в котором регулировка точности вычислений производится в соответствии со стандартом IEEE 854 [27], который является расширением стандарта IEEE 754 на случай ь = ю . Непосредственно в среде выполнения Maple изменение m производится с помощью присвоения значения глобальной переменной Digits. Например, если необходимо получить 25 значащих десятичных цифр потребуется следующая команда:

Digits:=25;

Ранее было показано, что погрешность округления значения функции представляется в виде:

А,=|ДМ-/И=Р11^, 41 или д, = сг д,, где сг - некоторый коэффициент коррекции погрешности, который при достаточно большом Am = m,+i - т, близок к 1.

Представим эти равенства в виде:

Р\ =1оёкА =log1o||c,||-fl»«lj

С,! 4.2

7,

Графики функций и

КША полагаем а = 1; для БША значение а может быть неизвестным и тогда его надо определить из численного эксперимента или, если имеется гипотеза о его значении, то это значение надо сравнить с экспериментальным.

Если в результате численного эксперимента установлено, что значение а соответствуют требуемой гипотезе, то это означает что достижима ЗТР, т.е. при достаточно большом Ат = т1+1 - т1 значение С можно оценить с приемлемой точностью: С = А^Ь"";

Значение погрешности решения для заданного т находится как Д( = СЬ~ат' = Д, ,+1, а требуемое значение мантиссы, обеспечивающее ЗТР определяется из условия: т

1 1 1 £

1--1оё7г а С

Но часто в реальном численном эксперименте найти достаточно большое количество решений может оказаться затруднительно, поэтому придется довольствоваться для оценки погрешностей округления ограниченным числом решений. В этом случае необходимо воспользоваться оценками погрешности решения Д, < <тД, где Д, = Д, ;+1, а значение относительной погрешности по формуле , < <т, рассмотренными ранее. Для гарантии выполнения условий сг;Д, =е1

4.1 Вычисление полинома

Рассмотрим задачу вычисления значения для следующего полинома по схеме «прямого» суммирования:

100 ;=1

Ожидаемая погрешность значений р(х) существенно зависит самого значения х и растет с ростом х, поскольку растет значение полинома 100

Р(|х|) = 2>|х|;. Для данного многочленна известна формула суммы: 1 ч п{— х)"+1 - (« + 1)(— х)" + 1 р (х) = -х—-—-——г^—--, которую можно использовать для х + 1) получения истинной погрешности.

Численный эксперимент для х = ^ показывает следующую зависимость:

1и I 1 1 1 1

1П15 I I 1 1 1 1

1 1 1

Ю10 ---- ----- 1 1 1 -1 ю5 1 1 1 1 1

1 1 1

10° 1 1 X 1 1 1 1

1 1 1 ю-5 1 1 1 1 1 1

1 1

Ю"10 1 1

1 1 ю-15 1 1 1

1 1 ю-20 1 х 1 1 ю-25 1

10-30 1 1

ЧГ.-35 1 I 1

10 15 20 25 30 35 40 45 50 55 60 т

Рисунок 1 Зависимость ошибки е в зависимости от т при вычислении полинома . . 5 р(х) = 2^(-\)Ч-х' в точке Х = —

1 з

Величина ошибки е{т) получена сравнением с известным результатом для суммы рп{х), вычисленным при т-100. Таким образом видно (Рисунок 1), что зависимость е(т) линейна. Сравнивая значения разности порядков е и соответствующей разности мантисс, получаем, что а-1. Далее рассмотрим зависимость погрешности от величины * для того же полинома. Будем брать значения х в диапазоне от 1 до 5 с шагом I

15'

Возьмём т = 35 и т=55. Вот как зависит точность получаемых значений от выбранной точки: ю41

10й

10" т=35

10"

10

10"

10"'

10 I иг *!-,'й .I. : II ' ч I) - г - - - г

УуШ ШШ

Ш-ЦЦ+н! и '> им г -.г "Ят ш

I I

7 9 ' I

1.5 2 2.5 3 3.5 4 4.5 5 X

1 1.5 2 2.5 3 3.5 4 4.5 5 X

100

Рисунок 2: Ошибка вычисления полинома р(х) = / ^\ (— 1)'/ • х' в зависимости от X

По графикам можно видеть, что алгоритм сложения точен при х = —.

Это связанно с тем, что исходные данные при основании ь = ю представляются точно и все задействованные операции дают точный результат.

Построим зависимость параметра погрешности С(х) от т при х = ~

С(х)

10"

10 15 20 25 30 35 40 45 50 55 60 т

100

Рисунок 3 с(х) = с(х,т) при вычислении полинома р{х) = ^ (— 1)' / • х'

1=1

Откуда видно, что величина С(х) является «почти» постоянной. Степень отклонения от постоянной величины можно видеть из таблицы: ш е С(х) Ш £ С(х) гп £ С(х)

10 1.93е+16 1.93е+26 30 0.000192 1.92е+26 50 1.93е-24 1.93е+26

11 1.93е+15 1.93е+26 31 1.93е-05 1.93е+26 51 1.93е-25 1.93е+26

12 1.93е+14 1.93е+26 32 1.92е-06 1.92е+26 52 1.93е-26 1.93е+26

13 1.93е+13 1.93е+26 33 1.93е-07 1.93е+26 53 1.93е-27 1.93е+26

14 1.93е+12 1.93е+26 34 1.92е-08 1.92е+26 54 1.92е-28 1.92е+26

15 1.93е+11 1.93е+26 35 1.93е-09 1.93е+26 55 1.93е-29 1.93е+26

16 1.92е+10 1.92е+26 36 1.93е-10 1.93е+26 56 1.92е-30 1.92е+26

17 1.93е+09 1.93е+26 37 1.93е-11 1.93е+26 57 1.93е-31 1.93е+26

18 1.92е+08 1.92е+26 38 1.93е-12 1.93е+26 58 1.93е-32 1.93е+26

19 1.93е+07 1.93е+26 39 1.93е-13 1.93е+26 59 1.93е-33 1.93е+26

20 1.93е+06 1.93е+26 40 1.93е-14 1.93е+26

21 192 1.92е+26 41 1.93е-15 1.93е+26

22 19300 1.93е+26 42 1.93е-16 1.93е+26

23 1930 1.93е+26 43 1.93е-17 1.93е+26

24 192 1.92е+26 44 1.93е-18 1.93е+26

25 19.2 1.92е+26 45 1.92е-19 1.92е+26

26 1.93 1.93е+26 46 1.93е-20 1.93е+26

27 0.193 1.93е+26 47 1.92е-21 1.92е+26

28 0.0193 1.93е+26 48 1.93е-22 1.93е+26

29 0.00193 1.93е+26 49 1.93е-23 1.93е+26

ЗАКЛЮЧЕНИЕ

В диссертационной работе получены следующие результаты:

1. Введены понятия элементарного алгоритма решений задач вычислительной математики (ЭАВМ), конечношагового (КША) и бесконечношагового (БША) алгоритмов. Для ЭАВМ получены оценки погрешности решений как для КША, так и для БША, зависящие от длины мантиссы и аргументов алгоритмов. Для КША получены оценки длины мантиссы, обеспечивающие достижение требуемой точности значения функции.

2. Введены понятия контрольного решения (КР) задачи, которое в оценках погрешностей решений с некоторой точностью заменяет истинное значение функции и g -устойчивой последовательности решений. Исследованы свойства g -устойчивости, в том числе доказана теорема о достижимости заданной точности значения функции; получены оценки погрешности, далее численно реализуемые в итерационной последовательности с переменной мантиссой (ИППМ).

3. Предложены алгоритмы, позволяющие оптимизировать процесс решения задачи в ИППМ. Рассмотрены методы округления полученных решений, причем округленное решение имеет заданную точность.

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

Список литературы диссертационного исследования кандидат физико-математических наук Гриневич, Алексей Иванович, 2013 год

СПИСОК ИСПОЛЬЗОВАННЫХ источников

1. IEEE 754-2008: 754-2008 IEEE Standard for Floating-Point Arithmetic -IEEE, NY, USA, ISBN: 978-0-7381-5753-5

2. Годунов С.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах - Новосибирск: Наука. Сиб. Отд-ние, 1988. - 456с. ISBN 5-02-028593-5

3. Higham N. J. Accuracy and stability of numerical algorithms -Philadelphia : Society for Industrial and Applied Mathematics, 1996.

4. Иванов В.Д., Косарев В.И., Лобанов А.И., и др. Лабораторный практикум "основы вычислительной математики" - Москва, Мз пресс, 2003.

5. Wilkinson J.H. Rounding Errors in algebraic processes - Englewood Cliffs, N.J.: Prentice-Hall, 1963, ISBN 0-486-67999-3

6. Федоренко P.П. Введение в вычислительную физику: Учеб. Пособие: Для вузов. - М.: Интеллект, 2008. - 504 с. ISBN 978-9-91559-011-2

7. Maples oft Maple User Manual - Maplesoft, 2011, ISBN 1-926902-07-6

8. Henrici P. Elements of Numerical Analysis - John Wiley & Sons Inc., New York, 1964.

9. The GNU Multiple Precision Arithmetic Library [HTML] (www.gmplib.org)

10.THE MPFR LIBRARY: ALGORITHMS AND PROOFS [HTML] (http://www.mpfr.org/algo.htm)

11.Ulrich Kulisch: Wissenschaftliches Rechnen mit Ergebnisverifikation. Eine Einführung, Vieweg-Verlag, Wiesbaden 1989, ISBN 3-528-08943-1

\2.Бахвалов H.C., Жидков Н.П., Кобельков Г.М. Численные методы

13.Software for Interval Computations collected by Vladik Kreinovich, University of Texas at El Paso [HTML] (http ://www. с s .utep. edu/interval-comp/main.html)

14.IEEE Interval Standard Working Group - P1788

15.iRRAM - Exact Arithmetic in С++ [HTML] (http://www.informatik.uni-trier.de/iRRAM/)

16.Воеводин В.В. Вычислительные основы линейной алгебры. - М.: Наука, 1977.-304с.

17 .Хемминг Р.В. Численные методы - Москва: Наука, главное отделение, 1972.

18.Рябенький B.C. Введение в вычислительную математику - Москва: Физматлит 2008

19.Langlois P. A Revised Presentation of the CENA Method, ARENAIRE -INRIA Grenoble Rhône-Alpes / LIP Laboratoire de l'Informatique du Parallélisme - INRIA - CNRS : UMR5668 - Université Claude Bernard -Lyon I - École Normale Supérieure de Lyon.

20.Demmel J. Effects of Underflow on Solving Linear Systems // EECS Department University of California, Berkeley Technical Report No. UCB/CSD-83-128 August 1983.

21.Demmel J. On Error Analysis in Arithmetic with Varying Relative Precision // IEEE 8th Symposium on Computer Arithmetic (ARITH), 1987 Page(s) 148-152.

22.Кнут Д. Искусство программирования. Том 1. Основные алгоритмы. 3-е издание, исправленное и дополненное - Вильяме, 2002 г. ISBN 5-8459-0080-8, 0-201-89683-4.

23.Кнут Д., Искусство программирования. Том 2. Получисленные алгоритмы. 3 - е издание - Вильяме, 2000 г. ISBN 0-201-89684-2.

24.Shouichi Matsui, Masao Iri An overflow/underflow-free floating-point representation of numbers. J. Inform. Process., 4(3):123 133, 1981.

25.Clenshaw С. W., Olver F. W. J., Turner P. R.. Level-index arithmetic: An introductory survey. In Numerical Analysis and Parallel Processing, Lancaster 1987, Peter R. Turner, editor, volume 1397 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, 1989, pages 95 168.

26. Wilkinson J. H. The Algebraic Eigenvalue Problem. // Oxford University Press, 1965. xviii+662 pp. ISBN 0-19-853403-5.

27.A Proposed Radix- and Word-Length-Independent Standard for FloatingPoint Arithmetic. // Reprinted in IEEE Micro, (August 1984): 86-100.

28.Марков А.А., Нагорный H.M. Теория алгорифмов - M.: Наука, 1984.

29.Петров И.Б., Лобанов A.M. Лекции по вычислительной математике: Учебное пособие - М.: Интернет-Университет Информационных Технологий; БИНОМ. Лаборатория знаний, 2010. ISBN 978-5-94774542-9.

30.Enright W.H. The Relative Efficiency of Alternative Defect Control Schemes for High Order Continuous Runge-Kutta Formulas - Technical Report 252/91, Dept. of Computer Science, University of Toronto, June, 1991.

31 .Verne J.H. Explicit Runge-Kutta Methods with Estimates of the Local Truncation Error. - SIAM Journal of Numerical Analysis, Aug. 1978.

32.Stone P. Jim Verner's "Maple" (dverk78) 13 stage combined order 7 and 8 Runge-Kutta scheme [PDF]

(http://www.peterstone.name/Maplepgs/Maple/nmthds/RKcoeff/Runge_K utta schemes/RK8/RKcoeff8c 2.pdf)

33.Curtis A.R. High-order explicit Runge-Kutta formulae, their uses, and limitations // J. Inst. Maths. Applies., V. 16 (1975) P. 33-55.

34.Фаддеев Д.К., Фаддеева B.H.. Вычислительные методы линейной алгебры,- Изд. 2-е.- М.: «Наука», 1963.- 656с.

35.Hammarling S.J. Latent Roots and Latent Vectors - The University of Toronto Press, 1970. ISBN 0 8020 1709 6

ЪЬ.Фиакко А., Мак-Кормик Г. Нелинейное программирование. Методы последовательной безусловной минимизации -М.: Мир, 1972.

37.Иванов В.В. Методы вычислений на ЭВМ: Справочное пособие. -Киев: Наукова думка, 1986. - 584с.

38. IEEE standard for radix-independent floating-point arithmetic: ANSI/IEEE Std 854-1987, 1987. [HTML] (http://grouper.ieee.org/groups/754/)

39.ISO/IEC 9899:1999 Standard for the С programming language (C99) -1999.

40.Clenshaw C. W. and Olver F. W. J. Beyond floating point. - J. Assoc. Comput. Mach., 31(2):319-328, 1984.

41 .Шокин Ю. И. Интервальный анализ. - Новосибирск: Сибирское отделение изд-ва «Наука», 1981.

42.Шарый С. 77. Конечномерный интервальный анализ. - М.: 2007.

43. Хансен Э,, Уолстер Дж. У. Глобальная оптимизация с помощью методов интервального анализа. - Москва-Ижевск: R&C Динамика, 2012.

АА.Жолен Л., Кифер М., Дидри О., Вальтер Э. Прикладной интервальный анализ (2-ое изд., испр.).- М.-Ижевск: Институт компьютерных исследований, 2007. - 468 с. ISBN 5-93972-585-6.

45 .Chaitin-Chatelin F., Fraysse V. Lectures on Finite Precision Computations. - Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1996.

46.Алефельд Г., Херцбергер Ю. Введение в интервальные вычисления. М.: Мир, 1987. 356 с.

47 .Воеводин В.В. Ошибки округления и устойчивость в прямых методах линейной алгебры. М., Изд-во МГУ, 1969, 140 с.

4%.Бирюков А.Г., Гриневич А.И. О гарантированной точности решений задач вычислительной математики в арифметике с плавающей

запятой и переменной длиной мантиссы // Труды МФТИ. - 2012. -Т.4, №3. - Р. 171-180.

49 .Бирюков А.Г., Гриневич А.И. Метод оценки погрешностей округления решений задач вычислительной математики в арифметике с плавающей запятой, основанный на сравнении решений с изменяемой длиной мантиссы машинного числа // Труды МФТИ. - Том 5, № 2(18), 2013, С. 160-171.

50.Гриневич А.И. Математическая модель погрешностей округления при вычислениях в арифметике с плавающей запятой и переменной длиной мантиссы. // Труды 55-й научной конференции / Управление и прикладная математика. Том 1. - М.: МФТИ, 2012, С. 15-16.

51 .Бирюков А.Г., Гриневич А.И. Итерационный процесс с переменной длиной мантиссы для решения задач вычислительной математики с заданной точностью // Труды 55-й научной конференции / Управление и прикладная математика. Том 1. - М.: МФТИ, 2012, С. 86-87.

52.Бирюков А.Г., Гриневич А.И. Анализ погрешностей некоторых алгоритмов БМ. Труды Института системного анализа Российской Академии Наук. Динамика неоднородных систем. Том 42(1), 2009, С. 106-111.

5Ъ.Бирюков А.Г., Гриневич А.И. Об эффективности и устойчивости численных методов решения систем нелинейных уравнений и задач безусловной минимизации // Труды Института системного анализа Российской Академии Наук. Динамика линейных и нелинейных систем. Том 25(1), М.: КомКнига, 2006, С. 111-123.

54.Демидович Б.П., Марон И.А. Основы вычислительной математики. -СПб: Лань, 2009, ISBN 978-5-8114-0695-1

55.Ulam S., Metropolis N. The Monte Carlo Method, — J. Amer, statistical assoc. 1949 44 № 247 335—341.

56Жиглявский А.А., Жилинскас А.Г. Методы поиска глобального экстремума. - М.: Наука, 1991

57.Zhigljavsky A., Zilinskas A. Stochastic global optimization. Berlin: Springer, 2008

58.Крылов В.И., Бобков B.B., Монастырный П.И. Вычислительные методы.Том I - М., Наука, 1976.

59. Математическая энциклопедия / Под ред. И.М. Виноградова - М.: «Советская энциклопедия», 1980.

60.Bailey D. Н., Barrioy R. Borweinz J. М. High-Precision Computation: Mathematical Physics and Dynamics. March 27, 2012.

61 .Robey R. W., Robey J. M., Aulwes R. In search of numerical consistency in parallel programming, Parallel Computing, vol. 37 (2011), 217-219.

62.He Y., Ding C. Using accurate arithmetics to improve numerical reproducibility and stability in parallel applications, J. Supercomputing, vol. 18 (Mar 2001), 259-277.

63 .Hayes W. Is the outer solar system chaotic? 11 Nature Physics, vol. 3 (2007), 689-691.

64.Borwein J. M., Bailey D. H. Mathematics by Experiment: Plausible Reasoning in the 21st Century // A.K. Peters, Natick, MA, second edition, 2008.

65.Simo C. Global dynamics and fast indicators // Global Analysis of Dynamical Systems, 373-389, Inst. Phys., Bristol, 2001.

ee.Jiang H., Barrio R., Li H., Liao X., Cheng L., Su F. Accurate evaluation of a polynomial in Chebyshev form // Applied Mathematics and Computation, vol. 217 (2011), 9702-9716.

61 .Lamb J. S. W. Reversing symmetries in dynamical systems // J. Phys. A: Math. Gen., vol. 25 (1992), 925-937.

68 .Bailey D. H., Crandall R. E. On the random character of fundamental constant expansionsfo // Exp. Mathematics, vol. 10 (2001), 175-190.

69.Hauschildt P., H. Baron E. The numerical solution of the expanding Stellar atmosphere problem // J. Comp. and Applied Math., vol. 109 (1999), 41-63.

10.He Y., Ding C. Using accurate arithmetics to improve numerical reproducibility and stability in parallel applications // J. Supercomputing, vol. 18 (Mar 2001), 259-277.

11.Frolov A. M., Bailey D. H. Highly accurate evaluation of the few-body auxiliary functions and four-body integrals. // J. Physics B, vol. 36 (2003), 1857-1867.

ll.Bailey D. H., Frolov A. M. Universal variational expansion for high-precision bound-state calculations in three-body systems. Applications to weakly-bound, adiabatic and two-shell cluster systems // J. Physics B, vol. 35 (2002), 42870-4298.

ПРИЛОЖЕНИЕ 1. Метод штрафных функций

Рассмотрим задачу нелинейного программирования:

X* =аг§пип/(х), где /: Я"-> Я'

при ограничениях

5.1

(р,(х) = 0,/ = 1,/, ф,(х) < 0,г = 1 + \,т

Суть семейства методов, к которым относится рассматриваемый метод штрафных функций [36], состоит в преобразовании данной задачи с ограничениями в последовательность задач минимизации без ограничений. Это преобразование выполняется с помощью соответствующей вспомогательной функции, которая выражается через функции задачи и определяет новую целевую функцию, имеющую безусловные минимумы в некоторой области. Постепенно изменяя параметр и тем самым увеличивая влияние ограничений на вспомогательную функцию, строят последовательность или семейство задач без ограничений, решения которых сходятся к решению исходной задачи.

Одной из реализаций описанного подхода является следующий метод:

х* = Нт(аг§тт.р(г,х)),

г—>°о

Р(г,х) = /(х) + г£р,2(х) + г ¿тах(0,?>, (*))2 5.2

1=1 ;=/+1

Рассмотрим следующие задачи:

Задача НП1:

Задача НП2:

1 т

fix) = ~x Hx-bx

n

JYx, =n2

I=]

1 J

f(x) = —x Hx-bx

и2-Ji'-x, <0

i=i

Конкретные реализация метода (5.2) выглядит так:

1.

е, = max

J 0 L 3 J cSrad

2. h = 10

-l"J

3. mi+1 = m, + 20,m, = 15

4. r, = 10,r, = 10/;.,/' = 104

Подход N m* /, сек к g

Clog 5 41 2.72 1.5883 le-14 1 4

linlO 35 2.483 1.04965e-14 6 3

lin20 35 2.172 2.42719e-13 6 2

2x 30 2.188 8.4681 Ое-12 6 2

Optm 30 5.344 4.05892e-12 9 4

mlOO 100 4.109 6.52007e-19 6 2

Clog 15 101 590.811 2.29593e-18 116 9

linlO 105 734.749 2.36129e-18 134 10

lin20 115 339.547 1.29264e-18 69 6

2x 240 218.578 7.46890e-15 51 5

Optm 101 28.656 4.16227e-16 4 2

mlOO 100 28.891 2.19711e-16 4 2

Clog 25 141 7210.58 2.78881e-12 336 14

linlO 175 9393.64 6.48244e-12 403 17

Ип20 195 5261.3 1.18178e-14 233 10

2х 480 1894.81 1.58310e-24 106 6

Optm 141 107.672 5.38002e-12 4 2

ml 00 100 1752.95 1.90953e-12 70 13

Таблица 20: Время решения задачи НП1 с помощью различных подходов.

Подход N m* t, сек II^MI к g

Clog 5 32 6.422 6.39319e-13 13 4

linlO 35 5.734 5.33276e-14 12 3

lin20 55 5.828 7.13687e-15 12 3

2х 60 5.858 7.13687e-15 12 3

Optm 32 5.406 4.43984e-12 7 2

mlOO 100 6.125 7.26702e-21 7 2

Clog 15 101 947.062 6.20323e-13 134 10

linlO 105 1047.2 6.26742e-12 133 10

lin20 95 477.251 7.13779e-12 64 5

2x 120 318.406 7.13779e-12 46 4

Optm 95 66.39 4.16227e-16 6 2

mlOO 100 68.125 2.31956e-17 6 2

Clog 25 104 6184.7 1.95185e-12 179 9

linlO 105 8158.19 3.56855e-12 213 10

lin20 115 4225.24 9.93039e-12 125 6

2x 480 3306.61 1.58310e-24 106 6

Optm

mlOO 100 2592.1 9.81700e-12 54 20

Таблица 21: Время решения задачи НП2 с помощью различных подходов.

Задача НПЗ:

/=1

п

Ег'"х< = *2

1=1

i=i

Подход N m* t, сек к g

clog 5 35 1.594 6.43935e-14 12 3

linlO 35 1.579 6.43935e-14 12 3

lin20 55 1.796 8.20816e-18 13 3

2х 60 1.814 8.20816e-18 13 3

Optm 35 3.048 3.61050e-16 16 4

ml 00 100 3.657 3.89910e-24 16 4

clog 15 42 13.454 5.29275e-15 17 3

linlO 35 9.688 2.03123e-15 12 3

lin20 35 12.188 3.11462e-14 16 2

2x 60 13 3.52487e-17 16 3

optm 35 20.719 8.00299e-15 15 2

ml 00 100 24.093 3.53841e-22 15 2

clog 25 41 39.171 1.93933e-14 15 3

linlO 35 37.484 9.60224e-13 14 3

lin20 35 37.813 1.43567e-13 15 2

2x 60 40.406 3.85343e-18 15 3

optm 35 64.156 3.27794e-14 15 2

mlOO 100 75.734 1.53715e-20 15 2

Таблица 22: Время решения задачи НПЗ с помощью различных подходов.

Задача НП4:

п i=i

«-¿(-1)4 < о

Подход N m* t, сек ||VF(x,) к g

clog 5 40 2.829 6.60066e-18 17 4

linlO 45 2.861 3.67842e-19 17 4

lin20 75 3.218 3.67842e-19 18 4

2x 120 3.031 3.67842e-19 17 4

optm 40 3.796 1.38507e-17 16 4

mlOO 100 4.564 4.56187e-25 16 4

clog 15 44 14.219 1.37324e-15 15 3

linlO 35 14.218 2.46418e-14 15 3

lin20 35 12.609 1.11779e-13 14 2

2x 60 14.796 1.50413e-17 15 3

optm 35 29.14 2.62535e-15 15 2

ml 00 100 35.157 3.53841e-22 15 2

clog 25 41 52.797 1.91131e-14 16 3

linlO 35 48.454 3.91658e-13 14 3

Нп20 35 49.422 4.13485е-14 15 2

2х 60 58.531 1.85659е-17 16 3

ор1:т 35 127.453 6.45038е-14 15 2

тЮО 100 166.359 1.53715е-20 15 2

Таблица 23: Время решения задачи НП4 с помощью различных подходов.