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

005058312

На правах рукописи

ГРИНЕВИЧ Алексей Иванович

МЕТОД ОЦЕНКИ ПОГРЕШНОСТИ ОКРУГЛЕНИЙ ЗНАЧЕНИЙ ВЫЧИСЛЯЕМОЙ ФУНКЦИИ, ОСНОВАННЫЙ НА ВАРЬИРОВАНИИ ДЛИНЫ МАНТИССЫ В АРИФМЕТИКЕ С ПЛАВАЮЩЕЙ

ЗАПЯТОЙ

Специальность 01.01.07 - вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание учётной степени кандидата физико-математических наук

1 6 МАЙ 2013

МОСКВА-2013

005058312

Работа выполнена на кафедре математических основ управления Московского физико-технического института (государственного университета) Научный руководитель: кандидат физико-математических наук,

доцент Бирюков Александр Гаврилович Официальные оппоненты: доктор физико-математических наук,

профессор Лобанов Алексей Иванович, кафедра вычислительной математики Московского физико-технического института (государственного университета), профессор

доктор физико-математических наук, доцент Рогов Борис Вадимович, Институт прикладной математики имени М.В. Келдыша РАН, ведущий научный сотрудник

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

им. A.A. Дородницына РАН

Защита состоится iKQ&,_2013 г. в .Я ^асов на заседании

диссертационного совета Д 212.156.05 при Московском физико-техническом институте (государственном университете) по адресу 141700, Московская обл., г. Долгопрудный) Институтский пер., д.9 ауд. 903 КПМ.

I'

С диссертацией можно ознакомиться в библиотеке Московского физико-технического института (государственного университета). Автореферат разослан_jCSoiMps?-A _2013г.

Ученый секретарь диссертационного совета

Федько О. С.

Общая характеристика работы

Актуальность темы

В настоящее время большое внимание научного сообщества уделяется вопросам «вычислений с высокой точностью», т.е. таким вычислениям, при которых возможны изменения длины мантиссы машинного числа (МЧ) в широком диапазоне значений. Как отмечено в обзоре D. Bailey за 2012 г., можно выделить целый ряд направлений исследований, где стандартной арифметики оказывается недостаточно. Среди них есть как давно известные проблемы, так и относительно новые, активное изучение которых началось вместе с появлением достаточных вычислительных мощностей. Это моделирование солнечной системы за период в миллионы лет, вычисление рекуррентных соотношений, определение экспоненциально малых явлений в динамических системах, компьютерное исследование новых математических соотношений, моделирование явлений в сверхновых звездах и др. В частности, А. Фролов утверждает, что «сейчас мы научились рассматривать и решать задачу нескольких тел с ограничениями, о чем нельзя было и мечтать всего несколько лет назад». Для решения задачи им использовалась арифметика высокой точности (120 знаков) для нахождения собственных векторов почти вырожденных матриц размерами 5000x5000 в рамках решения задачи п тел.

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

1. Решение плохо обусловленных систем линейных уравнений (СЛУ).

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

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

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

Время вычислений существенно возрастает (на 1-2 порядка) при переходе на «высокоточную» арифметику с варьируемой длиной мантиссы. Поэтому подобные вычислительные проблемы считались слишком трудоемкими до недавнего времени, и ученые старались изыскивать специфические методы решения для конкретных задач или упростить задачу до «разрешимого» варианта. В последние годы ситуация начала меняться с появлением библиотек программ с интерфейсами высокого уровня. Здесь следует отметить и специальные подпрограммы для нахождения значений функций (ARPREC, GMP, MFPR, MPFR++, MPFUN90, QD), пакеты компьютерной алгебры, позволяющие находить решение символически без потери точности (Maple, MathCad, Mathematica) и реализации языков программирования, позволяющие задавать длину мантиссы (LISP, Python, Perl, Haskell, Ruby). При переходе на высокоточную арифметику, как правило, оказывается возможным не переводить программу целиком, а заменять лишь часть ключевых алгоритмов на более «точные» варианты. Это позволяет локализовать вычисления, требующие высокой точности и минимизировать влияние возрастающего времени выполнения до приемлемых значений. Таким образом, безусловно, вычисления в арифметике высокой точности не могут рассматриваться отдельно от других подходов по оптимизации вычислений. Поскольку набор практических задач, где необходимы высокая или заданная точность, продолжает расти, возникает вопрос о построении метода (методики) оценок погрешностей округлений при варьировании точности арифметики для получения заданной точности решения (ЗТР). Приобретают актуальность вопросы выработки обобщенных подходов и методов, применимых к широкому спектру задач в условиях арифметики высокой точности.

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

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

1. Построение алгоритма вычисления значения функции f{x)eRk, применимого для достаточно широкого класса задач ВМ. Для предложенного алгоритма - исследование оценок погрешностей округления приближённых значений /(х), зависящих от длины т мантиссы машинного числа (МЧ) и получение такого значения т0, которое обеспечивает достижение требуемой точности решения.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Для достижения требуемой точности решения рассматриваются оценки вида:

д = ||7-/||<Ф1(/,...,/)<^, (1)

где f„i = \,l — приближенные значения решений. Предложен метод К-решений (Глава 2) для построения функций оценок погрешности Ф, Для

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

Определение 1. Пусть задано с > О — требуемая точность вычисления приближённого значения / функции /, т.е. если точность решения задачи достижима, то имеет место неравенство д = |/-/|<е или Л/||/|<с. Будем

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

Глава 1 диссертационной работы является вводной, носящей вспомогательный характер. Рассматриваются свойства машинного числа (МЧ) в формате IEEE 754; рассмотрены понятия машинного числа (МЧ), точности представления МЧ, функции округления МЧ и арифметических операций; рассмотрены ситуациии потери значимости и переполнения, возникающие при выполнении арифметических операций; рассматриваются характеристики библиотек програм GNU GMP и GNU MPFR, в которых реализуется вычисление в арифметике с плавающей запятой базовых (стандартных)

математических функций. Важнейшим свойством этих программ является то, что пользователь может задавать длину мантиссы МЧ в очень широком диапазоне от mm¡n =8 до mml> = 646456993 десятичных знаков. Появление в свободном доступе программного обеспечения с такими характеристиками расширяет границы возможностей для получения решений широкого круга задач ВМ с гарантированной точностью высокого порядка. Диссертационная работа посвящена изучению погрешностей округления решений, значения которых получены для варьируемых значений длины мантиссы.

Также рассмотрены другие подходы для достижения требуемой точности решений. В частности, рассмотрены свойства модели рациональных чисел и модель чисел с многоуровневой экспонентой. Отмечается, что в широкой практике они пока не нашли своего применения.

В главе 2 диссертации предложен алгоритм решения задачи ВМ, названный «элементарным» и исследованы некоторые его свойства.

В §2.1 вводятся определения понятий, на основе которых строится «элементарный алгоритм ВМ» (ЭАВМ).

Определение 2. Пусть <p\R" -> й' - некоторая функция, хеЛ" - вектор, хт- его машинное представление. Тогда: <р(х), <р(хт) - значения функции <р в точках х и Хт, <рт(хт) — машинное представление значения ip(xm); <рт(х) = <рт(xm). Значения функций <р{х) и <р(хт) в общем случае представляются бесконечным числом знаков; в этом случае будем говорить, что их значения получены в точной арифметике. □

Определение 3. Математические функции и операции, реализуемые в библиотеках программ, реализующих математическую библиотеку стандарта ANSI С99, назовем базовыми или стандартными. К базовым операциям относятся: округление чисел, арифметические операции, логические операции, операции вычисления математических функций, таких как sinx, cosx, tan*, и их обратные, ех, а', х", Iog0 х и т.д. □

Одной из библиотек, реализующих базовые функции стандарта С99, является GNU MPFR, в которой при заданной длине мантиссы m для значений базовых функции выполняется округление до последней значащей цифры, т.е.:

(*т ) = <PÍXт )(1 + и) , Н S 5,, <рт (хт ) Ф о, где (2)

<Р(хт) и <?,„(хт) - значение одной из стандартных (базовых) функций, вычисленное в точке хт е МЬ т р. Для случая, когда <рт (х1П) = 0,

к'ОО-РтО^ИЧ^тЦ^о- В (2) ¿>1 = "> ¿-основание системы исчисления, 30-погрешность нуля; она на много порядков меньше (см. Главу 1).

Определение 4. Задачей вычислительной математики будем называть совокупность понятий и условий /•"(/,*, т, г), определяющих возможность

вычисления значений вектор-функции /: Я" Як, где /(х) - решение данной задачи в точке хеС; е - точность решения, т - длина мантиссы. Условие достижения точности означает: найти значение вектор-функции /т(х) для которой выполнено одно из условий:

ф(/„, .-,/„, г или < г или |/. - /|| < г. □

Функция ф(») определяется методом решения задачи; /„,, Д, ..., /га — приближенные значения решений вычисляемые в точке хе6 для длин мантисс МЧ т,т,,...,т;.

Алгоритм решения задачи ВМ представляет собой композицию (объединение) базовых операций, т.е.

АЮ{Пх)) = <р\а1)ч<р1(а2)ч-чср''-1 (вдч) V(а„). (3)

АШ(Г(х)) = <р\а1)ч<р1{а2)ч-ч<р,''-\а^)ч1р''' К), N -> со (4)

где символ V означает объединение операций. Первый алгоритм вычисления /(*) является конечношаговым (КША), второй - бесконечношаговым (БША). Оба алгоритма (3) и (4) рассматриваются в точной арифметике. Соответствующие КША и БША алгоритмы для вычисления /„(*) представлены следующими формулами:

А1С{Гт(х)) = <р'Ла^1р1{а1)у-^Ч>:-\а^<р:(а,). (5)

ЛЮ(/т (*)) = <р\ (а,) V <р1 (аг) V • • ■ V рУ (ал, ) V (ол. ), (6)

- число шагов БША, которое определяется по некоторому правилу его окончания.

Определение 5. Алгоритм вычисления функции / е Л' (5) в точке хеЯ" при длине мантиссы ш будет называться элементарным алгоритмом для решения задач вычислительной математики (ЭАВМ), если вычисленное значение функции в точной арифметике по алгоритму (5), в котором логические операции не выполняются, дает точное значение функции /(х), т.е. имеет место:

л¿G(/W) = pl(o1)vf)2(a2)v••■ví)''-,к_l)v^>iVЮ □ (7)

В §2.2 изучаются погрешности округлений, возникающие в итерационном процессе ЭАВМ.

Лемма 1. Пусть <р, - базовые вычислительные операции (кроме логических) из некоторой библиотеки программ, ге[1,7У6] - номер базовой операции. Тогда значение <р'т(а,) можно представить в виде

< (а,) = Ч>' (о,) + а А> гДе (8)

И6 - число базовых операций библиотеки, ]а11<(а,)[, <5, = а, еЯ', а, =ат)

- вектор машинных чисел. □

Лемма 2. Пусть для чисел у, г и их приближенных значений с!т, ут, гт выполнены условия Ас/ = с!т - </ = , Ау = ут - у = ссу, Аг = -г = ¡)г \ ут= цуЪ', 2т = Ь1, [лу и рг - мантиссы чисел ут и гт, г - порядок числа; - цг= г\т_ч Ь"', 1 < <7 < т, >/,„_„ - мантисса числа //у - ц., где ш- <? её длина; | = , а = 53,, р = /35{.

Тогда погрешность д = —^---— = ь' х.5\> гДе X = —~— % -——— -о

у-2 ЬЧ у-2 )

Лемма 3. Пусть функция <р: Я" -> Д1 удовлетворяет условию Липшица:

+ Д х)— < ¿|Д , х,х + Ахей, где (9)

(7 с= /Г — компакт. Тогда существуют такие числа /, е [- Ь, Ь\, что

^(х + Д х) - ^(х) = ^ /,. Д х,.. □

ы

Теорема 1. Пусть функция /(х)еЯ1, хеСс/Г; базовые функции (р'(а:) (кроме логических) либо являются функциями округления числа, либо непрерывны по Липшицу, т.е. в некоторой окрестности Ц(д) точки я, е Л'удовлетворяют условию: (р\а1 )-<р'(Ь1 )<Ц 10,-6,1, /6[1,Л'], а алгоритм (3) вычисления функции /(х) является элементарным для х еО. Тогда существует такой вектор Сей', что

/т(х)-/(*) = СЯ1, где ¡5, =-1ь'"т, т - длина мантиссы. □ (Ю)

Определение 6. Пусть в бесконечношаговом алгоритме выполнено N первых базовых операций и для /(х) получено приближение значения вычисляемой

функции (х). БША вычисления функции / называется элементарным, если для всех указанных N алгоритм вычисления значения /„,"(*) будет элементарным. □

Определение 7. Бесконечношаговый алгоритм называется сходящимся, если /(*) = lim /"(*), где /v (х) - точное решение задачи после N операций БША. □

Л'—»ас

Из доказательства теоремы 1 следует, что для ЭАВМ решение задачи представляется как:

/.{*)= Л*)+Сб„ (11)

где С - вектор параметров погрешностей округления, [|с||<со, S, т -

длина мантиссы МЧ. Анализ формулы, приведенный в Следствии к теореме 1 показывает, что решение (11) возможно представить в другом виде:

/Лх)~ f(x)+Cb~am, где 0 < а < 1 и с . (12)

Если БША решения задачи является сходящимся т.е. /(*)= Нт/"М> гДе /Л(х)

v

— точное решение задачи после выполнения N операций его определяющих, то результат теоремы 1 уточняется.

Теорема 2. Пусть для функции f"(x) БША является элементарным, сходящимся и выполнены условия теоремы 1. Тогда существуют векторы CeR1, yeRk, такие, что Vi > О: < е имеет место представление

/mNW = /W+C<5,+y □ (13)

В §2.3 получены оценки длины мантиссы, гарантирующей достижение требуемой точности ЭАВМ решения задачи ВМ для КША и БША. Определение 8. Будем говорить, что метод (алгоритм) вычисления значения функции /ей* называется корректным (КМ), если для любого е>0 найдется такой размер мантиссы ш, что:

А = \\т- /„ «IH 6, и ли < £ , при ||/( х)\\ * 0 . (14)

Величину е в (14) будем называть требуемой точностью. □

Определение 9. Будем говорить, что значение функции /т (х) имеет погрешность порядка а, 0<а<1, относительно погрешности мантиссы ¿>;,, если существуют константы С и С, такие, что \fxeGcR", |С|<С, С/|/(х||<С0, Vm > mmin, где mmin - минимальное значение длины мантиссы при которой могут проводиться вычисления и

Л = ||/(— fm(х)|| ä CS" , для абсолютной погрешности; (15)

= |. А „ < С0<Г при II/(х)|| * 0, для относительной

||/м|| ii/wi 0 " р ^

погрешности,

где 5fl = b"m . □

Векторы С в формуле (10) (или С в (12)) назовем параметром погрешности (ПП) значения функции

Теорема 3. Пусть погрешности Дзначения функции f(x) или имеют

порядок а. Тогда для любого е > 0 данный метод вычисления функции будет

корректным при m >

а С

или ш>

, 1 i £

1--'оёьТГ

а С„

, где W - целая часть

числа А . □

Таким образом теорема 3 (и 4) дают условия, при которых имеет место заданная точность решения (ЗТР). Конечно, верхнее значение длины мантиссы, необходимое для обеспечения ЗТР, не может превышать максимальных значений, которые имеет используемая библиотека программ. В частности для пакета Maple mniax =65535 десятичных знаков, а для GNU GMP mmax = 646456993 десятичных знаков.

Теорема 4. Пусть БША является сходящимся и элементарным, для каждого N выполнены условия теоремы 2, погрешность округления для каждого N в (13) имеет порядок а\ т.е. СЗ^=СЬ~а", 0<а< 1 и выполнено условие ||с|<С VxeG, Ge«", Vm > mmin, где mmm - минимальная длина мантиссы при которой могут проводиться вычисления. Тогда существуют такое число базовых операций N и такая длина мантиссы ш, при которых достигается требуемая точность решения с, т.е. ||/ту(л-)-/(д:|<£. □

Глава 3 диссертации посвящена построению вычислимых оценок погрешностей округлений значений fm(x).

В §3.1 третьей главы вводится определение итерационной последовательности с переменной мантиссой (ИППМ).

Определение 10. Совокупность L решений задачи fm (*) при значениях длины мантиссы m¡,/e [l.X.]: m,<m2<...<mt назовем итерационной

последовательностью с переменной мантиссой (ИППМ) решения задачи ВМ.

□

Определение 11. Числа rj и ?j0 называются малыми по сравнению с 1, если О < 77 < 770 < 0,1. Условие малости чисел по сравнению с 1 обозначается символом ««»-много меньше: 77 «1, 70 «1. □

Определение 12. Пусть значения погрешностей решений равны Д/=||Л>)-/М||. |' = [U]; У>'. Ue[lL], К-решением

(Контрольным решением, KP) задачи ВМ называется значение (*), если Д,.=Д,.1+<?,ХД,., «б[1,1-1], (16)

где | ¿¡„ l« 1, т.е. | | малое число по сравнению с 1. □

Определение 13. Обозначим g, = ^^ = 1/""'.^ } .1, /e[l,Z,-ll. - --А,

Последовательность решений задачи назовем g -устойчивой, если g, < g„«1, / е [l,£-l]. Число g,, ; е [l,I-l] назовем коэффициентом уменьшения (КУ) погрешности на ¡-м шаге. □

Для упрощения изложения представим значение погрешности для частного случая его порядка а = 1:

л = /„(*)-(17) где Ь - основание МЧ. Тогда очевидно, что при достаточно большом значении т, и |с„| < С, где С не зависит от т, ||д|| может быть малым числом.

Лемма 4. Пусть для ИППМ выполнено условие (17) и погрешности: Д'=|/т,(х)-/(х| = |Ст.||ь~"'' и А"= ||/„-(*)-Д*| = Цс».-!6""" удовлетворяют условию:

Ы^-соМ, Vrn',m":m, <т'<т" <mL. Тогда найдется такая Am = min(mj+l - m,.),

Ir»'!

iе[i,z,-i]], что g, <g0, «1, где g„- заданное малое число, т.е. ИППМ g-устойчива. □

Приведем достаточное условие g -устойчивости ИППМ. Рассмотрим некоторые оценки погрешностей решений задач ВМ для g -устойчивой ИППМ. Лемма 5. Пусть ИППМ^- устойчива. Тогда для значений погрешностей д(, Дj, Д,j имеют место двухсторонние оценки:

Л Л !> Л у Л , , , .

—— <Д. < ——; ——— < Д < ——— \ -г.к<д. <1+г.К, П81

где _/ е [1,1]. □

Теорема 5. 1. Пусть в ИППМ выполнено условие —— < «1, / >;; е [|Д].

Для того, чтобы значение функции /т (х) было К-решением, необходимо и достаточно, чтобы ИППМ была £-устойчивой. 2. Пусть ИППМ ^-устойчива. Тогда для любого е > 0 существует такое решение /„, (х), что д, < £ и

Из П.2 теоремы следует, что существует такой номер г ИППМ, для которого имеет место ЗТР Х/е > 0. Номеру соответствует длина мантиссы т,, которая, конечно, не должна превышать значения тп]Ш библиотеки программ.

Препятствием к практическому применению полученных оценок является значение , в котором присутствует «истинное» решение /(х), в общем случае неизвестное при вычислениях.

Для практического применения указанных оценок вводятся новые понятия. Определение 14. Пусть /„ (х) - КР задачи ВМ. Обозначим

Д- ( ) ( 11 ^ _ ни _ I ■".-! -"■< 41 ^ , е П 2]. Последовательность решений задачи ВМ

назовем квазиустойчивой (или ^-устойчивой по отношению к КР), если g¡■ «1. Число g¡■ назовем КУ значений Д(1. □

В леммах 6 и 7 приводятся оценки, связывающие теоретические значения для коэффициента уменьшения g. и вычисляемые на практике величины ^.

В теореме 6 рассматриваются оценки погрешностей округления для g-устойчивой ИППМ, которые можно применять при численном решении. Теорема 6. Пусть решение/(х), хеД", /еД* оценивается значением /„ (х), / е [1, £ -1], ИППМ ^-устойчива, |/„. (х| > ||д,|, /е [1,1]. Тогда имеет место оценка

оА ^ „ а

шШг^ШГ (19)

где ;>,, М,+1.ь], а,-^; "

корректирующие множители погрешности решений. Машинное значение числа

<4 к + 6,,.„ |—имеет относительную погрешность я—^—Ь', которой можно

пренебречь. □

В §3.2 рассматриваются подходы к практической оценке погрешностей, получаемых в процессе построения ИППМ. Исследованы алгоритмы построения оценок погрешностей округлений.

В разделе 3.2.1 рассматривается алгоритм последовательной оценки погрешности округлений решений. Задается некоторое начальное значение длины мантиссы т = тх. Следующие значения длин мантиссы задаем по правилу т„,=т;+Дтм, 7 = 1,2,.... Особенностью настоящего алгоритма является то, что оценка погрешности округления (ОПО) задачи ВМ проводится после каждого решения с мантиссой большей длины. При решении задач ВМ возможны различные схемы получения ОПО. Выделим две основные схемы. Пусть /„ (д:) и /т (х) - решение и К-решение задачи ВМ.

Введем числа <тл = —1— и ег0 = - коэффициент коррекции погрешности,

где g¡j <g0 и а0 « .

1. Пусть задана требуемая точность решения сл и для некоторых решений (х), /„,(х) выполнено неравенство: д,,<<гА&и =Д' <гл. Тогда полученная относительная погрешность решения е' удовлетворяет условию:

е.

- = £' <

¡ДМЦО"««) ||/т>|(1-«о)"

(20)

2. Пусть задана требуемая точность решения £а для относительной

Д. сг0Аа

погрешности, т.е. выполняется неравенство: .. . < ^—-^-¡т < . 1огда для абсолютной погрешности решения д' выполняется условие:

а; Д'<4/„ (*)||(1(21)

В оценках (20) и (21) значению функций Ф,(у) из (1) соответствует функция Ф,(/ ,/ )=(т/<Д;. для абсолютной погрешности и функция Ф2(/„ )=-¡р^пт —

' ' ' 1КИ1

для относительной погрешности.

П1 (Вариант 1) В этом варианте<т0 = --i-г, дтм = Am = const, / = 1,2,..., и для

(l-g0Xl-e0)

А--.и

некоторого i > 1 выполнено условие <т0 „—^^-тг = е,< сй, а условия е, < е0 не

выполнены для 1</</-1. Решением задачи является значение /„, (дг), абсолютная погрешность его не более еа: Д, <A,i+l/(l-g0) относительная -¿-„.

е Д.

Функции оценки ЗТР (1) Ф,(-,) здесь равны: Ф ,(•,•)=,, 0 для относительной

К. И

погрешности и Ф, (•,■)= А"4' -для абсолютной. 1-го

П2 (Вариант 2) Пусть определены т,, т2=т1+&т, решения / (*) и /mj (*) и условие а0 А'2 < е0 не выполнено, число р, = ЛР6'"'. Оценим значение т3 из

IM*!

условия Д3 = рф~щ = г0|/т . Считаем, что р, = да, где ^ можно брать равным

7д

.2Г = Ю2, J = 103 И Т.Д. Если log, 12 ^,<2Ат, ТО т, = т,+Ат.

уД

Если log,, —¡jVbli > 2Am' то тз=т2 +

Со||ДМ||

1 Х^ р log

. Находим / (л:) и проверяем

условие <т0.. < г0. Если это условие выполнено, то / (*) - решение задачи

К-И

ВМ и значение д' =д23/(1-^0). Если не выполнено, то далее ИППМ реализуется ПО Варианту 1 при /и,, = /и + 1Ат, ;>3.

В разделе 3.2.2 (П2) рассматривается табличный алгоритм оценки погрешностей округления решений. Пусть ИППМ решений представлена в виде: /т{х), / = [1д], тм=тк+\Ат, /е[1,X, — 1 ]. Значения т,, Ат, число решений ь задает Вычислитель (лицо, решающее задачу ВМ). Построим таблицу погрешностей Д., у > /, г 6[1,1-1], /е[2,/.]. Таблица значений дгД.Г -треугольник значений р., на / строке которого находятся приближенные значения р., у > /. Строится также таблица чисел '¿^ и таблица оценок

относительной погрешности решения: е^= <т0 й—¡г > )>и е [1,/,], где

<т0 = 1 + , /0 < 0,1 и т.п. Совокупность всех указанных таблиц дает информацию о величине погрешностей решений / (*), / е [1,1,-1], и об % -устойчивости ИППМ.

В разделе 3.2.3 рассматриваются методы округления решений. Полученное значение функции представляется как /„.(*) т,-значное, т.е. представление решения в 10-ичной системе имеет т, 10-ичных знаков. Значение /„, (*) округляется до t знаков, t < т,-

В §3.2 рассматривается естественное обобщение сформулированных выше правил округления на случай к > 2, т.е. на случай вектор-функций.

В §3.3 исследуется метод оценки погрешности округления скалярных решений по правилу «совпадения t первых десятичных знаков».

В §3.4 показано, что при g-устойчивости БША для них справедливы все результаты теории, сформулированные в разделах 2 и 3, а потому методика получения заданной точности решений для БША будет той же, что и для КША.

В §3.5 рассматривается вопрос эффективности предложенного подхода. Метод КР обладает значительной универсальностью в определении погрешностей округления значений функции, т.к. он не ориентирован на какие-либо классы решаемых задач. Таким образом, если метод КР решает задачу за приемлемое время, а традиционный метод (ТМ), использующий «стандартное» программное обеспечение решает, но не дает ЗТР, то метод КР можно считать высокоэффективным по сравнению с ТМ.

В главе 4 анализируются результаты численных экспериментов, проделанных для проверки предложенных теоретических результатов. Для практических целей использовался математический пакет Maple 13, в котором регулировка точности вычислений производится в соответствии со стандартом IEEE 854, который является расширением стандарта IEEE 754 на случай ь = ю. Ввиду ограничений на размер автореферата, результаты приводятся для наиболее значимых экспериментов.

В §4.1 рассматривается задача вычисления полинома. Рассмотрена зависимость погрешности от т и точки х.

В §4.2 рассматриваются различные варианты методов решения СЛУ: метод Гаусса с выбором главного элемента. Рассмотрим задачу нахождения решения системы линейных уравнений (СЛУ):

Hz = с (22)

где н - матрица Гильберта порядка кт.е.:

й=К).иб[и]Л~. (23)

К=10 К=20 К=30 К=40

-А

.Л.

V

100 150 200 250

Рис. 1. Левый график: Зависимость абсолютной погрешности решения СЛУ (23) Д,„ от длины мантиссы т. Правый график: Зависимость / - времени решения от длины мантиссы т и размерности СЛУ (К = 10,20,30,40)

Из первого графика (Рис. 1) видно, что зависимость погрешности для плохо обусловленной матрицы имеет зону неопределенности, после чего выходит на участок эксспоненциального убывания вместе с ростом длины мантиссы. На этом участке проявляется свойство g -устойчивости.

45 50 55

Рис.2. Зависимость р в области стабильности (ш > 60) и области роста (т <60)для СЛУ вида ((23) от длины мантиссы шири К =40 и шаге изменения

длины мантиссы дт=1.

На графике (Рис. 2) представлено значение параметра погрешности р = \г„-г\ьт при 6 = 10 и известном точном решении г. Из графиков видно, что зависимость параметра погрешности р от т при ,я>45 можно трактовать как постоянное значение, на которое наложено «случайное» возмущение. Рассмотрим различные правила варьирования длины мантиссы МЧ (ПВМ) для для нахождения решения задачи ((23) с заданной точностью гг0 =ю 20 в условиях ИППМ:

Для всех правил т, =100. Критерий остановки варьирования: Дм м, < 0,1 и д,у+1 < " это грубые условия, характеризующие выход из «зоны

неопределенности», g!^ - <go = o;qi - это оценка g -устойчивости данного ИППМ;

Н

ПВМ №1 clog: Правило, на основании П1.2 Вариант 2, описанного в п. 3.2.

т. , = ш + Дот, ти. = от, + min

Ю,-

log4

, 7 = 103 - выбран на основе анализа

графиков для р, рассмотренных выше. ПВМ №2 НпЮ: ПВМ №2. тм = от, + Дот,Дот = 10 ПВМ №3 Нп20: тм = т. + Дот,Дот = 20 ПВМ №4 2х: ПВМ №4. тм = 2т,

ПВМ №5 ор1ш: т. = т - единственное вычисление при «необходимой» длине мантиссы.

ПВМ №6 тЮОО: т, = 1000 - правило с достаточно большой длиной мантиссы и единственным вычислением.

Эксперимент проведён для матриц Гильберта порядка N=100,200,300.

Подход N /, сек mi mL £ A,x i

clog 100 30.563 175 254 7.48402e-103 2.53773e-25 6

linlO 55.234 180 200 9.86545e-49 6.83909e-29 11

lin20 36 180 220 1.85175e-68 6.83909e-29 7

2х 32.954 200 800 1.34404e-648 9.86545e-49 4

optm 4.922 175 175 2.53773e-25 1

ml 000 18.203 1000 1000 1.47395e-848 1

clog 200 542.266 350 499 2.25154e-194 1.46094e-46 9

linlO 908.375 330 350 1,46094e-46 2.87443e-27 16

lin20 579.719 340 380 2.96797e-75 1.53653e-36 10

2x 589.953 400 1600 2.33753e-1296 4.24806e-96 4

optm 57.063 350 350 1.46094e-46 1

ml 000 168.859 1000 1000 4.28177e-696 1

clog 300 3651.16 500 638 5.88136e-181 8.38046e-42 11

linlO 6545.81 480 500 8.38046e-42 1.19533e-21 21

lin20 3819.86 480 520 1.94937e-61 1.19533e-21 12

2x 3349.88 600 2400 1.76402e-1941 9.43380e-142 4

optm 248.844 480 480 8.38046e-42 1

шЮОО 517.891 1000 1000 8.60282е-542 1

Таблица 1. Зависимость времени решения ? для СЛУ вида ((23) от ПВМ при N = 100,200,300, где ш, - достаточная длина мантиссы, е - истинная погрешность полученного решения, - практическая оценка погрешности, / - число шагов

и ИППМ.

Из результатов эксперимента (Таблица 1) видно, что подход с варьированием мантиссы позволяет получать решение с требуемой точностью вне зависимости от способа варьирования. Далее в работе для решения СЛУ с матрицей Гильберта рассмотрен метод сопряженных градиентов. Произведены оценки коэффициентов gfj и их сравнение с теоретическим результатом. На основании этого сделан вывод о #-устойчивости данных методов. Рассмотрено поведение параметра погрешности с(х,т).

В §4.3 рассматривается задача численного нахождения производной. Показано, что параметр погрешности а = 1/2 при соответствующем выборе шага дифференцирования И ~ ^¿Г. Изучен вопрос &-устойчивости ИППМ для этой задачи.

В §4.4 задача нахождения собственных чисел методом Данилевского. На эксперименте показано, что ИППМ для метода Данилевского является g-устойчивой начиная с некоторой длины мантиссы т .

В §4.5 рассмотрены задачи безусловной минимизации и решения систем нелинейных уравнений. Для метода Ньютона приведены формулы для предсказания Ат при варьировании мантиссы на основании правил изложенных в главе 3. Рассмотрены другие правила к варьированию мантиссы и показана эффективность предложенного метода предсказания. В частности, произведен следующий численный эксперимент.

Задача БМ Б (Функция Розенброка 1)

/« = ¿(1 -х,У+а

, где (24)

а> 0, Р = 1,2,3.... Для функции (24) известно точное решение — х' =1,/ = 1,«, а значение /(*') = 0, причем решение (точка *') у этой функции - единственное. При численном эксперименте принимались следующие общие условия'. а) Алгоритм заканчивал работу, если Ах1[ = \\х„рк\<е, и где

е1 =£2=10"20 - требуемая точность. При такой величине требуемой точности, стандартной арифметики двойной точности уже недостаточно и варьирование мантиссы становится необходимым.

б) Шаг X, находится из условия =а^тт/(дг4_|+Лрк). Одномерная

минимизация происходит с помощью метода золотого сечения. Причем

критерий остановки одномерной минимизации - погрешность АЛ<Ь 1-3-1

Для нахождения требуемой длины мантиссы использовались различные ПВМ. Общий подход к реализации ИППМ состоит в использовании последовательности длин мантисс т8 =т, <т2 <т3..., где ^ = 1,2,... - группы вычислений с указанным значением мантиссы. Для каждой группы вычислений определены следующие параметры: 1. Требуемая точность для данной группы:

10 ,

. 2. Параметр перехода. Способ определения точности для

каждой группы: = Цт,,*,,/7). Для первой группы т, = 7] (*„,/•',£,). 3. Шаг

численного дифференцирования в зависимости от т : Ъ= 10

При вычислениях в каждой группе используются следующие условия окончания: Условие окончания 1: Если ||.г4то происходит увеличение мантиссы и вычисления продолжаются для следующей группы начиная с достигнутого хк. Условие окончания 2: Если количество шагов в данной группе превосходит л = лг, то происходит переход к следующей группе. Условие окончания 3: Если на любом шаге группы срабатывает более сильное условие полного окончания а) то алгоритм считается завершенным.

Рассмотрены ПВМ, аналогичные приведенным ранее для СЛУ:

ПВМ №1: (с1ой): тм = + шах

0,—

, т, =15-

ПВМ №2: (НпЮ):тм = т, +10, т,=15 ПВМ №3: (1т20):тм = т, + 20, „,, = 15 ПВМ №5: (2х):

т,ч| = 2 ■ «7;, /и, =15

ПВМ №5: (т*): использует т - наименьшую из достаточных точностей полученную на одном из предыдущих шагов. = т1, т, = т . ПВМ №6 (шЮО): В этом подходе фиксированная достаточно большая длина мантиссы 100, т.е. тм=т,,т- 100.

Подход N т /, сек 1у/Ы[ л-/' к г

с1ОБ 5 52 2.628 9.28103е-26 4.64229е-21 9.80325е-55 55 31

НпЮ 75 5.016 4.68859е-36 3.41326е-28 4.23552е-76 33 7

Ип20 135 5.19 4.68869е-36 3.41303е-28 4.23540е-76 33 7

2х 960 5.72 4.68877е-36 3.41264е-28 4.23510е-76 33 7

ор!т 52 6.781 4.55641е-29 1.04900е-33 1.02147е-61 33 7

тЮО 100 8.327 6.98179е-54 1.34388е-33 1.35510е-111 33 7

с1с^ 10 43 13.727 4.17797е-23 8.68153е-21 1.56645е-51 94 50

1т10 105 26.046 2.52116е-34 2.61099е-23 9.24818е-74 55 10

Пп20 195 28.048 2.52112е-34 2.61112е-23 9.24846е-74 55 10

2х 7680 78.735 2.52108е-34 2.61126е-23 9.24875е-74 55 10

ор^п 43 28.563 9.12954е-22 2.66973е-23 4.99117е-48 54 10

тЮО 100 33.469 3.96826е-41 2.69291е-23 2.72010е-87 55 10

Таблица 2. Время решения задачи (24) с ЗТР е, = е2 = 10~20 с помощью различных

ПВМ, градиент вычисляется численно. - погрешность иолучепного

решения.

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

В заключении приведены основные результаты диссертации.

В приложении 1 рассмотрены различные подходы к варьированию мантисс применительно к методу штрафных функций.

В приложении 2 рассмотрены сеточные методы для решения задачи Коши и уравнения теплопроводности. В частности, исследован вопрос о возможности связывания шага сетки А с длиной мантиссы т по аналогии с задачей о численном дифференцировании.

Основные результаты диссертации

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

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

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

Автор выражает глубокую благодарность своему научному руководителю доценту Бирюкову А.Г. за постановку задачи, ценные советы, терпение и постоянное внимание к работе.

Список публикаций автора по теме диссертации

1. Бирюков А.Г., Гриневич А.И. О гарантированной точности решений задач вычислительной математики в арифметике с плавающей запятой и переменной длиной мантиссы // Труды МФТИ. — 2012. - Т.4, №3, С. 171180.

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

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

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

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

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

Гриневич Алексей Иванович

МЕТОД ОЦЕНКИ ПОГРЕШНОСТИ ОКРУГЛЕНИЙ ЗНАЧЕНИЙ ВЫЧИСЛЯЕМОЙ ФУНКЦИИ, ОСНОВАННЫЙ НА ВАРЬИРОВАНИИ ДЛИНЫ МАНТИССЫ В АРИФМЕТИКЕ С ПЛАВАЮЩЕЙ ЗАПЯТОЙ

Автореферат

Подписано в печать 12.04.2013. Формат 60 х 84 Усл. печ. л. 1,25. Тираж 100 экз. Заказ № 120 Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Московский физико-технический институт(государственный университет)» Отдел оперативной полиграфии «Физтех-полиграф» 141700, Московская обл., г. Долгопрудный, Институтский пер., 9

Московский физико-технический институт (государственный университет) Кафедра математических основ управления

На правах рукописи

04201356199

МЕТОД ОЦЕНКИ ПОГРЕШШ

ГРИНЕВИЧ Алексей Иванович

К 519.8

ОКРУГЛЕНИИ ЗНАЧЕНИИ ВЫЧИСЛЯЕМОЙ ФУНКЦИИ, ОСНОВАННЫЙ НА ВАРЬИРОВАНИИ ДЛИНЫ МАНТИССЫ В АРИФМЕТИКЕ С ПЛАВАЮЩЕЙ ЗАПЯТОЙ

Специальность 01.01.07 «Вычислительная математика»

Диссертация на соискание учётной степени кандидата физико-математических наук

Научный руководитель кандидат физико-математических наук, доцент Бирюков А.Г.

МОСКВА-2013

СОДЕРЖАНИЕ

ВВЕДЕНИЕ................................................................................................5

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

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

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

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

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

1.4.1 GNUGMP...............................................................................37

1.4.2 GNUMPFR.............................................................................38

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

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

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

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

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

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

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

2.4 Выводы.........................................................................................60

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

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

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

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

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

алгоритм........................................................................................................90

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

3.6 Выводы.........................................................................................93

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

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

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

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

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

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

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

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

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

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

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

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

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

ЗАКЛЮЧЕНИЕ.....................................................................................133

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

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

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

Задача Коши.......................................................................................146

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

Сокращения

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

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

запятой)

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

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

математики

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

мантиссой

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

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 стала осново�