CC BY
33
УДК 519.677
А. Г. Бирюков1, А. И. Грипевич1'2
1 Московский физико-технический институт (государственный университет)
2000 «ГринМарк»
Метод оценки погрешностей округления решений задач
вычислительной математики в арифметике с плавающей запятой, основанный на сравнении решений с изменяемой длиной мантиссы машинного числа
Статья посвящена вопросам анализа погрешностей округления решений задач вычислительной математики на ЭВМ в арифметике с плавающей запятой и переменной длиной мантиссы машинного числа. Предложен метод оценки погрешностей округления, основанный на сравнении решений с различной длиной мантисс, сформулированы правила достижения требуемой точности.
Ключевые слова: погрешность округления, точность решения задач вычислительной математики, машинное число с переменной длиной мантиссы, бесконечноша-говый и конечношаговый алгоритмы, К-решение задачи вычислительной математики, гарантированная точность решений задач.
1. Введение
В работе авторов [1] предложен метод анализа погрешностей округления решения задач вычислительной математики (ВМ) в арифметике с плавающей запятой и переменной длиной мантиссы машинного числа (МЧ).
Настоящая статья является продолжением работы [1]; в ней предлагаются численные оценки погрешностей округления решений задач ВМ, гарантирующие достижение заданной точности, и дается их теоретическое обоснование.
Проблема анализа влияния погрешностей округления на решение задач ВМ актуальна со времени появления ЭВМ и остается таковой по сей день. Научные исследования над указанной темой ведутся в разных направлениях. Отметим классические работы по исследованию погрешностей округления при решении задач линейной алгебры [2,3]; исследование погрешностей округления в рамках интервального анализа [9,10]; статистический анализ погрешностей округления [6, 11]; исследование новых моделей по выработке машинного числа [7]; алгоритмы с автоматической коррекцией ошибок округления первого порядка — метод CENA [8].
В настоящей работе предлагается метод оценки погрешностей округления решений задач, отличный от рассмотренных в приведенной выше литературе. В соответствии с [1] под решением задачи понимается значение некоторой вектор-функции f (х) £ Rk, х £ Rn. Для оценки погрешностей решений вычисляется несколько их значений fmi (х) при различных длинах мантиссы МЧ: mi < Ш2 < ...mi < ...тЭти решения имеют погрешности A¿ = Wfmi (%) — f (ж)||) i £ [1, L]. Значение погрешности решения A¿, i £ [1,L — 1] оценивается значением погрешности A¿¿ = ||/mi (%) — fmL (ж)| i £ [1,L — 1], при этом предполагается, что Al << Ai. Здесь и далее под нормой ||^|| понимается евклидова норма
INI = \]S=i z £ ® разделе 2 рассматриваются свойства некоторого решения задачи fmL (х) при значении длины мантиссы т = m,L, названного контрольным или К-решением (КР). В разделе 3 вводится понятие гарантированной точности решений и рассматриваются правила оценок погрешностей округлений. В разделе 4 рассматриваются оценки погрешностей округления по правилу совпадения первых десятичных знаков (СПЗ)
решений с различной длиной мантиссы. В разделе 5 рассматриваются правила оценки погрешностей решений для класса бесконечношаговых алгоритмов (БША). В разделе 6 кратко рассматриваются вопросы эффективности метода К-решений. В разделе 7 на примерах решений систем линейных уравнений и задачи численного дифференцирования приведены результаты численного эксперимента, иллюстрирующие основные свойства метода КР.
Пусть / (х) - точное решение и / (х) - приближенное решение некоторой задачи ВМ.
Понятия абсолютной и относительной погрешностей А/ =
А/
||/(*)|| Рас~
сматриваются по отношению к / (х) - точному, но неизвестному значению решения. В ВМ
f (х) — f (х)
. Не рассмат-
используется также другое определение относительной погрешности: А// ривая специфики этих определений, в статье используется её первое значение. Незнание значения f (х) является принципиальным ограничением для получения оценок погрешностей А/ и А// ||/1|. Использование машинных чисел с возможностью изменения длины мантиссы при решении задач ВМ позволяет предложить вариант устранения этого принципиального ограничения.
2. К-решения задачи ВМ и их свойства
Определение 1. Совокупность L решений задачи fmi (х) при значениях длины мантиссы m,i,i £ [1,L\: т1 < т2 < ... < ть назовем итерационной последовательностью с переменной мантиссой (ИППМ) решения задачи. □
Определение 2. Числа ^ и ^о называются малыми по сравнению с 1, если 0 < ц ^ ^о ^ 0,1.
<<
ц << 1, щ << 1 □
Условие ^о << 1 задается Вычислителем (лицом, решающим задачу ВМ). Пусть а\ > 0, а2 > 0 ai = Щ2: V << 1) тогда а\ - малое число по сравнению с а2, т.е. а\ << а2. В вычислительной практике, в зависимости от требований задачи, число ^о ^ 0,1 может меняться в широком диапазоне значений. Например: ^о = 0,1; ^о = 0, 05; ^о = 10-fc, к = 1, 2,... и т.д. Определение 3. Пусть значения погрешностей решений равны Aj = Hfmi (х) — f (ж)! i £ [1,L]', Aij = ||/mi (x) — fmj (ж)||, j > i, i,j £ [1,L]. К-решением (КР) задачи ВМ называется значение fmL (х), если
А^ = AiL + iiLAiL,i £ [1, L — 1], (1)
где |&l| << 1, т.е. малое число по сравнению с 1. □
По смыслу К-решение означает «Контрольное решение», т.е. решение, позволяющее оценить значение погрешности решения fmi (х). В работе [1] (теорема 1) было получено значение погрешности функции f (х) при вычислении её с длиной мантиссы т. Для упрощения изложения представим его для частного случая порядка погрешности а = 1:
А = fm (х) — f (х) = СтЬ-т, (2)
где b - основание МЧ. Введем следующее
Определение 4. Обозначим gi = ^т1 = М1™г+1/^п11, г £ [1,L — 1]. Последовательность
решений fmi (х) назовем д-устойчивой, если д^ ^ до << 1, i £ [1,L — 1]. Число gi, г £ [1,L — 1] назовем коэффициентом уменьшения (КУ) погрешности на i-м шаге. □
Обобщением КУ д^ является чиело gij = г < j, г £ [1, L — 1]. Для него справедливо:
9ij = 9j-i9j-2...9i << 1,9%j < 9%,j £ [2, Ц , (3)
если ИППМ g-устойчива. Приведем достаточное условие ^-устойчивости. Лемма 1. Пусть для ИППМ выполнено условие (2) и погрешности: А' = Ш (х) — f (ж)|| = ЦСт, || Ъ-т' и А" = ||/т» (х) — f (ж)|| = ||Ст» || Ъ-т" удовлетворяют условию: п^Г"м ^ £о = const, Ут',т'' : m1 ^ т' < т" ^ ть- Тогда найдется
такая Am = min (m^+i — mi), г £ [1, L — 1]]; чт,о gi ^ до, до << 1, где до- заданное малое
число, т.е. ИППМ g-устойчива. Доказательство
Доказательство очевидно: Am ^ 1 — logb ^ , где LAI _ целая часть числа А □
Приведем без доказательства следующее утверждение.
Лемма 2. Пусть V1,V2 £ Rn и ||V2|| < тогда имеет место неравенство
||vi|| — HV2II < iiVi — V2||. □
Рассмотрим некоторые оценки погрешностей решений задач ВМ для ^-устойчивой ИППМ. Ввиду ограниченности объема, приведем следующую лемму без доказательства. Лемма 3. Пусть ИППМ g-устойчива. Тогда для значений погрешностей Ai; Aj, A^ имеют место двухсторонние оценки:
А < \ < А, < 1 - дг])Аг < (1+ дг])Аг, (4)
1 + 9í.j 1 - 9í.j 1 + дц 1 - дц '
где i,j £ [1,Ц- □
Замечание. Из (4) при малых значениях gi j: gij ^ до << 1 можно получить важные практические оценки:
Aij^ AirnAj = gijAij. (5)
Теорема 1. 1. Пусть в ИППМ выполнено условие ^ д0 << 1, j > i; i,j £ [1,L\. Для того, чтобы значение функции fmj (х) было К-решением,, необходимо и достаточно, чтобы ИППМ была g-устойчивой. 2. Пусть ИППМ g-устойчива. Тогда для любого е > 0 существует такое решение fmi (х), что Ai ^ е и A^ ^ (1 + g0) е, j > i. Доказательство.
1. Доказательство п. 1 теоремы не приводится, ввиду ограниченности объема статьи.
2. Т.к.. ИППМ д-устойчива, т.е. gt ^ д0, t = 1, 2, ..i, ..j-, д^ ^ д0ш Ai = дц,Ai = (gig2...gi) Ai, то Ai ^ д0Ai. Потребуем, чтобы Уе > 0 выполнялось неравенство Ai ^ д0Ai ^ е. Решая неравенство, найдем номер итерации, на которой гарантировано достижение требуемой точности решения: г = 1 + • Из (4) следует: A^ ^ (1 + д^) Ai ^ (1 + д0) е. □
Замечание. Определение (3) называет К-решением задачи ВМ fmL (х) - значение функции f при наибольшем заданном значении т = m¿. В доказанной теореме КР - значения fmj (х), j £ [2, L], оцениваемым решением является fmi (х), г £ [1,L — 1]. Из теоремы 1 следует важный практический вывод: если ИППМ д-устойчива, то достижима любая требуемая точность решения. На практике теорема 1 п.2 применима до максимального значения мантиссы, которую обеспечивает данная библиотека программ. В частности, для GNU MPFR mmax = 646 456 993 десятичных знаков.
Введем понятие значения КУ, в котором точные значения f (х) не используются. Определение 5. Пусть fmL (х) - КР задачи ВМ. Обозначим
9i = A/+1,¿ = IImT¿+1 (x) jmL ; i £ [1,L — 2]. Последовательность решений задачи
гЬ 11 Jmi\x)~JmL (ж)1
ВМ назовем квазиустойчивой (или ^-устойчивой по отношению к КР), если gL << 1. Число gL назовем КУ значений AiL. □
r\r г L L ДiL ufmi (x)—f™L(x)| • ^ • т/лг
Обобщением параметра gL является чиело gL = -г*— ^ tr^—т—гттг> J > i - Ку в
J ll J™i (x) — JmL (x)||
ИППМ для любой пары i,j £ [1,L — 1]. Очевидно, имеет место равенство
4 = gi_ igi_ 2 ...gÍL, (6)
Причем дL << 1, если g¡L << 1, i ^ t ^ j — 1.
Теперь необходимо сравнить значения gk и gij-
Лемма 4. Пусть ИППМ д-устойчива. Тогда для чисел gíj и gk справедливы оценки
< ^ < *> . <7>
1 + 9iL 1 - 9iL
Если же ^ 9о << 1, то существует такое число что | ^ д0 и
9^ = 9ц (1 + $) . (8)
Доказательство не приводится ввиду ограниченности объема статьи. □
Получим более простые приближенные двусторонние неравенства типа (7). Лемма 5. Пусть для мантисс выполнены, условия: тг < тз < т^ г, 3,1 € [1, Е — 1], т,1 — тз = т^ — тг = Дт и выполнено неравенство (7), в котором при Е = I значениям,и ди можно пренебречь. Тогда имеют место оценки
-1 + 1 - -
где ш = ^г. Если же << 1, то оценка (9) имеет вид
(1 — < дгз < д1гз (1 + ид{?) . (10)
□
Рассмотрим теперь оценки погрешностей округлений решений (х). Теорема 2. Пусть решение / (х), х € Кп, / € Кк оценивается значением /т1 (х), г € [1,Е — 1], ИППМ д-уст,ойчи ва, \\fmi (ж) || > \\Дг\\, ^ € [1,Е]. Тогда имеет место оценка,
Дг3 ^^Г < ®2 ,, Дг] , , (11)
(ж)Р ||/(ж)|р 2\\fmj (ж)\\'
где з >1,3 € [г + 1,Ь], Ч = Щ^Щ! ^ = ^ = (1-а^-)(!-«,) "
тирующие множители погрешности решений. Машинное значение числа м ^ ц имеет,
\\jmj И\\
отпноситпбльную позрсшносшь ^ —2—
ь1
которой можно пренебречь.
□
3. Правила оценки погрешностей округления решений задач ВМ.
Конечношаговый алгоритм (КША)
ИППМ можно рассматривать как метод решения задачи ВМ, позволяющий получать оценки погрешностей округления. Полезно следующее:
Определение 6. Пусть задано число е > 0 - требуемая точность решения задачи ВМ, то есть если точность решения задачи достижима, то имеет место неравенство А = || fm (х) — f (ж)|| ^ £ ми А/ ||/ (ж)|| ^ е. Будем говорить, что имеет место гарантированная точность решения (ГТР) задачи ВМ, если задача решается на ЭВМ методом, для которого известны оценки погрешностей решения, гарантирующие достижение указанной точности и значения оценок погрешностей определяются вместе с искомым решением. □
В теореме 1 доказано, что в ИППМ существует такое решение fmi (ж), которое при определенных условиях обеспечивает ГТР задачи ВМ.
В КША решение задачи ВМ получается за конечное число базовых (стандартных) вычислительных операций [1]. Оценки погрешности округления (ОПО) решения задачи ВМ проводятся путем сравнения значений fmi (ж) и fmj (ж), полученных при длине мантисс т = mi vi т = m,j. При этом возможны различные варианты построения ОПО. 3.1. Алгоритм последовательной оценки погрешности округлений решений. ИППМ решения задач ВМ считаем ^-устойчивой. Задается некоторое начальное значение длины мантиссы т = т\. Следующие значения длин мантиссы задаем по правилу
т%+1 = тг + Дтг+1, г = 1, 2,.... Особенностью настоящего алгоритма является то, что ОПО задачи ВМ проводится после каждого решения с мантиссой большей длины.
Пусть Ь = 1, т.е. найдено только одно решение задачи ВМ при длине мантиссы т = т-1, не исключающей одинарную, двойную и четверную точности или точности при большей длине мантиссы. Возможны ситуации, когда решаемая задача по мнению Вычислителя проста и нет необходимости находить её решение при других значениях длины мантиссы. В рассматриваемом случае за оценку погрешности решения ответственна интуиция Вычислителя, но указать величину погрешности решения в общем случае не представляется возможным.
при числе решений Ь ^ 2. При достаточно большом Дт ИППМ ¿»-устойчива (лемма 1) и по теореме 1 (х), ] > г, - К-решение по отношению к (х). Тогда по теореме 2 абсолютная Дг и относительная Дг/ ||/ (ж) || погрешности приближенного решения (х) имеют вид
Дг - ||т (ж) - { (ж)|| < 11 (Х1)_ (х)" - <
1 - г 1 - г
Дг ^_Дч_^ _
У (х)| " (1 - ^)(1 -«.) II ¡т, (х)|| ^ 0 || ¡т, (х)|| ,
где Дц и ц^^^ц ~ ОПО, ал = т—^) а0 = ^коэффициенты коррекции ОПО
для абсолютной и относительной погрешностей соответственно; д^ ^ до ^ 0,1 << 1, а.у ^ а0 ^ д0 ^ 0,1. Более точно, а^ << д^ и а0 << д0.
При решении задач ВМ возможны различные схемы получения ОПО. Выделим две
основные схемы. Пусть ¡т1 (х) и (х) - решение и К-решение задачи ВМ.
1. Пусть задана требуемая точность решения ел и для некоторых решений ¡т1 (х), /т. (х)
выполнено неравенство: алДц = Д1 ^ £л- Тогда полученная относительная погрешность 1
_ 1 ^ £а
= ^ < иг / -7 = £о. (13)
11 ^ (х)! (1 - а0) ^ ^ (х)! (1 - а0)
0
выполняется неравенство: ^—а ) ^ £0' Тогда Для абсолютной погрешности Д1
оаАц = А1 < е0 | | fmj (ж) || (1 - ао) = £а. (14)
Ниже рассматриваются два варианта алгоритма ОПО, в которых используются относительные погрешности решений (2-я схема); достижение требуемой точности решений гарантируется теоремой 1.
Пусть числа до ^ 0,1 и ао ^ до- Задавая различными способами Am^+i, г = 1, 2,..., будем иметь варианты алгоритма ОПО.
3.1.1. Вариант 1. Пусть ИППМ ¿-устойчива и ао = (i_go)1i_aoу Am^+1 = Am = const,
г = 1, 2,..., и для некотор ого г ^ 1 выполнено условие ао и \ц = £г ^ £о> a условия
£t ^ £о не выполнены для 1 ^ t ^ г — 1. Решением задачи является значение fmi (х), абсолютная погрешность его (14) не более £а, относительна я - ео', в (14) определено значение A1
лемм 4, 5.
3.1.2. Вариант 2. Пусть определены т1: т2 = т1 + Дт, решения /т_1 (ж) и /т_2 (ж) и
условие оо и ,А1.2 . и ^ £о не выполнено, число р1 ~ Д12ЬтК Оценим значение тз из условия \\/т2 О^у
Дз = РзЬ-тз = ее \\/„2 (ж)\\. (15)
Считаем, что рз = ХРъ гДе X можно брать равным % = 102, % = 103 и т.д. Из (15) получим %Д 12Ьт1 Ь-тз = ее \\/т2 (ж)\\; откуда Ьтз = Ь{ ^ • Если —^ 2Дт,
£0\/т2 (ж)\ £0\/т2 (х) \\
то т3 = т2 + Дт.
Если 1о& ^Йщ >2Д тз = т2 + 1о& ^ЙЬ
Находим /тз (ж) и проверяем
условие оо и . А3 ^ £е- Если это условие выполнено, то ¡т2 (ж) — решение задачи ВМ и
значение Д1 = ооД23- Если не выполнено, то далее ИППМ реализуется по Варианту 1 при тг+1 = тг + гДт, г ^ 3.
3.2. Табличный алгоритм оценки погрешностей округления решений. Пусть ИППМ решения задачи ВМ представлена в виде: (ж), г = [1,£], тг+1 = + гДт, г € [1,Ь — 1]. При выполнении Ь решений (ж), г € [1, Ь], становится возможным полученную информацию представить в виде таблиц. В частности, возможно построить таблицы значений Дг^, рг = \,Ст,11| = ДгзЪт% д^, £г^ и т.д. Совокупность всех указанных таблиц дает информацию о величине погрешностей решений (ж), г € [1,Ь — 1], и о ¿'-устойчивости ИППМ. '
3.3. Округление решений задач ВМ. Часто бывает так, что решение /т (ж) с требуемой точностью возможно представить при меньшем числе десятичных знаков. Такое представление реализуют различные процедуры округления чисел. Решением задачи ВМ в общем случае является вектор - набор чисел ¡т1 € Кк.
3.3.1. Рассмотрим сначала случай, когда к = 1, т.е. решением задачи является число. Далее, для упрощения записи обозначим А = / (ж) - точное решение, а = (ж) - приближенное решение - Шг-значное число. Следуя [12], напомним известное понятие.
а
числа удовлетворяет условию
\А — а\ < 110е-; (16)
являются верными в широком смысле, если выполнено условие
\А — а\ < 1 ■ 10е-, (17)
□
Округление решения проводится по известному правилу «по дополнению». Схема получения округленного решения следующая:
а) Пусть известна оценка решения \А — а\ ^ Д1 (см. пп. 3.1.1, 3.1.2).
б) Из решения неравенства Д1 ^ 210е-ь определим число верных знаков решения:
I = е + [— ^ю 2Д1] и е — 1 = — [— ^ 2Д1] . (18)
а1
погрешности
\А — а1\ < 1 ■ 10-^ 1о§ю2А^. (19)
а
А
а
Д1 : \А — а| ^ Д1) а1 - другое приближение числа а, такое что а = а1 + а, где \а| < \а|.
а1
IA — ai| = IA — a + а| < А1 + |а| = А2. (20)
Пусть мантисса числа а имеет т десятичных знаков. Представим числа a, ai, а в виде
(S* ioj
( Es * ■10-) ■106= ± ( iL Si ■10-) \i=1 / \i=i+1 /
а = ±ß ■ 10е = ± ( > ji ■ 10-i ) ■ 10е = ai + а,
(21)
ai = ± I > 's* ■ 10-i I ■ 10е,а = ± I > Si ■ 10-i I ■ 10е,
где e - порядок числа, 1 ^t^m.
Определение 8. Разбиение (21) числа а на ai и а назовем сечением числа по мантиссе. В (21) i-значное число ai - округленное значение числа am — i-значное число, а -погрешность округления числа ai, |а| - ошибка сечения числа a. Способ округления числа сечением по мантиссе назовем методом отбрасывания [17]. □
Для приближенного числа a оценка А2 = а1 + |а| состоит го численной т-значной оценки а1, полученной в процессе решения задачи вм и |а| - ошибки течения числа a. Число
А2
ai а
очевидно зависят от t, что можно представить как a*, а*. Значение чиела to, гарантирующее достижение точности решения еа, рационально определять из следующего условия:
t0 = min i, есл иА2 = А1 + |а*| ^ €а и 1 ^t ^ т. (22)
В свою очередь для удобства практического использования (экономичности записи) число А2 = А1 + |а*° | может быть округлено сверху и это значение А2 не должно превышать значение 8а — требуемой точности решения. □
3.3.2. В разделе 3.1 рассмотрены алгоритмы нахождения решения задачи ВМ, которое может быть как скаляром (к = 1), так и вектором (к ^ 2). Случай округления решения fmi (х) при к = 1 рассмотрен в п. 3.3.1. Правила округления компонент решения f^ (х), j £ [1, к], к ^ 2 могут быть различными. Они зависят от тех требований, которые предъявляет к решению Вычислитель. Рассмотрим два варианта:
1. Задача ВМ решается в соответствии с алгоритмами пп. 3.1 и 3.2 и полученные значения решений fmi (х) не округляются.
2. Задача ВМ решается в соответствии с алгоритмами п. 3.1 и, кроме выполненных требований точности, должны выполняться требования точности для компонент решения fmi (х) после их округления по правилу «по дополнению» или по методу сечения. Требования точности для абсолютной погрешности для компонент ej, j £ [1,к] должны удовлетворять условию J^jLi ^ £A>TRe £а ~ требуемая точность решения из п. 3.1. Более подробного изложения этого правила, а также других возможных вариантов в рамках этой статьи приводить не будем.
4. Оценка погрешности округления по совпадению первых десятичных знаков (СПЗ) решений с различной длиной мантиссы
В работе [4] автор предлагает способ достижения требуемой точности решения (перевод наш): «В тех случаях, когда основным источником погрешностей является округление, общий подход к оценке точности вычисления таков: пересчитать результат с помощью более точной арифметики и сравнить количество совпавших знаков в первом и втором случае. Интуитивно мы предполагаем, что требуемая точность результата может быть достигнута при вычислениях с достаточно точной арифметикой». Обоснование предлагаемого способа оценки точности решения в работе [4] не приводится, идея высказана на уровне интуитивного предположения. Приведем описание способа оценки точности решения, основанного на учете числа первых совпадающих десятичных знаков решений fmi (х) и fmj (х).
Введем соответствующее
Определение 9. Пусть дано I значений функции (j (ж), ] € [1,1},1 ^ 2 одинакового порядка, представление которых в десятичной системе имеет вид: (j (ж) = ± ^ • 10-;^ Т0е, ] € [1,1], г -натуральное число или г = го. Будем говорить, что у I функций одинакового знака совпадают £ первых десятичных знаков (СПЗ), если в 1 = в2 = ... = вг € [1, ¿], 1,2. □ Пусть £ (ж) € В,1, т.е. / - число, и получены два решения ¡т1 (ж) и (ж) при длинах мантисс т^ ту т, < ту Используя (2) при Ь = 10, запишем:
(ж) = / (ж) + Сгщ 10-т;, (ж) = / (ж) + Сгщ 10-т. (23)
В общем случае числа / (ж), Ст 10-т% Ст,. 10-т - бесконечнозначны, а (ж) (ж), ¡г (ж) кк (ж) - конечнозначны. Представим их в виде:
/ (ж) = ± ^ 5 ,10-; 10е; f (ж) = ? (ж) + к (ж);
£
=1
=1
/4(ж) = ± V • 10- • 10е; к (ж) = ±
Чг=¿+1
10-
10е
А,; = Ст. 10-
=
гп10~
10
,е-и
(24)
Аj = Ст3 10-т = ±
Е
10е- Ь-;
/т; (ж) = ±
(т; \
Е ^10-а]
=1
10е; и, (ж) = ±
10е
=1
где — Ь, и — ^ - порядки погрешностей чисел /т; (ж) и /т, (ж) Знаки + ми - у чисел /, /тц одинаковы. Знаки у чисел А,, Аj могут быть различными. Порядки погрешностей решений /т; (ж) и /т. (ж) при ¿-устойчивой ИППМ удовлетворяют условию: tj ^ ^ 1. Очевидно, что число первых совпадающих десятичных знаков может изменяться от 0 до и. Например, V чисел т (ж) = 0,4001111 и (ж) = 0,3999989, А, = 0,1111 • 10-3, Аj = —0,11 • 10-5, / (ж) = 0, 4, нет ни одного совпадающего десятичного знака. Однако
чается часто и этот вариант полезно использовать при оценке погрешностей решений. Итак, пусть у решений /т; (ж) и /т■ (ж) совпало í первых десятичных знаков. Тогда в качестве решения берется число ^ = ^ (ж) = ± га=1 ^а10-а) 10е. Погрешность числа (ж) оценивается по методу сечений, рассмотренном в п. 3.3.1 (теорема 3).
Совпадение Ь первых десятичных знаков у /т; (ж) и /т. (ж) ещё не означает, что совпадают £ первых знаков у решений / (ж) и /т; (ж) / (ж) ж /т. (ж). Например, пусть / (ж) = 0, 4, тогда у решений (ж) = 0, 3999888 и /т. (ж) = 0, 3999998 совпадают 4 первых знака и нет ни одного совпадающего знака с решением / (ж).
Рассмотрим условия, при которых у решений / (ж) /т; (ж) г ^ 1 < Ь совпадают £ первых десятичных знака, 1 ^ Ь < т,.
Теорема 4. Пусть Ыт = к1т (ж) = /т (ж) — /г(ж), кт (т — £)-значное число,
\к1п\ < 1 • 10'
IАI < 1 • 10
( ж)
т ( ж)
0 < шт — АI < 1 • 10'
(25)
причем если Ыт и А имеют разные знаки, то должно выполняться неравенство
\Кт — АI < 1 • 10'
если одинаковые, то \кт\ > \АI.
Погрешность решения /г (х) удовлетворяет условию I /1 (х) - /(х) I < 1 ■ 10е-
( х)
/ (х) = /4 (х) + к (х) = /* + к = /т (х) - Д = /4 + - Д. (26)
Из (26) следует, что к = Ыт — Д. Так как 0 ^ | к| < 1 ■ 10е-4, то условие (25), как следует из (26), эк вивалентно тому, что функции / (х) и /т (х) имеют £ СИЗ. Из неравенств | ккт | < 10е - и | Д | < 10е - следует, что условию | кгт — Д | < 1 ■ 10е - соответствует случай, когд а к1т и —Д имеют одинаковые знаки (т.е. \кт и Д имеют разные знаки), а усло|ию 0 ^ | ккт — Д | , при одинакотых знаках чисел кгтш Д, эквивалентно неравенство | к!т | ^ | Д | . Неравенство | / (х) — (х)| < 1 ■ 10е - следует из (25) и (26). □
Следствие 1. Пусть у функций / (х) и /т (х) найдено Ь СПЗ. Тогда для любого 1 ^ ^ ^ имеет место неравенство |/ (х) — /в (х)| < 1 ■ 10е-5. □
Для практического применения неравенств | кгт — Д | < 1-10е— и | кт| ^ | Д| необходимо знать оценки погрешностей Д, т.к. значение Д в общем случае неизвестно.
Но оценку погрешности Д^ некоторого решения /т (х) возможно получить, только зная К-решение /т^ (х), ] > г. Введем об означения: /т1 (х) = /г (х) + к\ (х) = /1 + Щ, /т (х) = /(х) + Дг, /т, (х)_ = /Чх) + к* (х) = /' + к,, /т, (х) = / (х) + Д„ Дг, = /т (х) — /т^ (х) = /Дг — Д,- Так как ИППМ ¿-устойчива, то /тj (х) - К-решение и выполнены равенства Дг = Д^ + ДД.^- = Д^ (1 + (1), Д, = ^-Ду, |^ . ^ д0 ^ 0,1 при д^ ^ 1/11 (теорема 1). Представим далее: пусть 1 ^ £ ^ ¿¿и
к = и* 10е -; к* = ^10е-; к, = ^ 10е-;
Д^ = 10е-*4 = щ, 10е-;,] >г,г < ¿¿;
| г/1 | < 1 | VI | < 1 — 10-т
^ 1 — 10 т\ | VI, | = ^г, 10* \ - мантисса числа | Дг,
округленная до т^ знаков, т.е. ^ ^ 1 — 10-т\ Теперь критерий для £ СПЗ (теорема 4)
□
Следствие 2. Для того чтобы решения / (х), /т1 (х), /т■ (х) имели £ СПЗ, необходимо и достаточно, чтобы
0 < |— (1 + ) V, | = | V, — ^ | < 1 — 10-т. (27)
Имеет место тождество: — (1 + ) = V* — □
Учитывая неравенство || ^ д0, из (27) получим достаточные условия £ СИЗ, которые можно использовать на практике:
| г!| + (1 + до) Щц < 1 — 10-т%ли + до^ц < 1 — 10-т, (28)
при и V, и имеющих разные знаки и
| V- | ^ (1 + до) ^ 10^или ^ до% 10^, (29)
при VI и VI,, и Vгj, имеющих одинаковые знаки. Для определения знака погрешности Д, = ДДнеобходимо найти решение /т1 (х), I > Если решение /т1 (х) не найдено, то
о
/т1 (х)-
5. Правило оценки погрешности решений задач ВМ. Бесконечношаговый алгоритм
В работе [1] введено понятие сходящегося и нормального бесконечношагового алгоритма решения задач ВМ (БША) и доказано, что решение задачи (х) после выполнения некоторых N базовых вычислительных операций представляется в виде
fN (x) = f(x) + C5 ! + 7, С51 = fZ (х) -Г (х), 7 = Г (х) — / (х), (30)
где fN (х) - точное значение решения после выполнения N операций, m - длина мантиссы МЧ, ¿1 = !&1-Z. По своей структуре БША являются итерационными алгоритмами, т.е. очередное приближение решения определяется после выполнения дополнительных Ns базовых вычислительных операций, где s - номер итерации, N = ^^=1 Ns, L - число решений N
N
условия теоремы 1 из [1], погрешность округления в (30) имеет порядок а : Сôi = СЬ-ат, 0 < а ^ 1 и выполнено уел овие ||С || ^ С Vx G G, G С Rn, У m ^ mmin, где mmin - минимальная длина мантиссы, при которой могут проводиться вычисления. Тогда существуют такой
N m
требуемая точность решения е, т.е. || fZ, (х) — /(х)|| ^ е. Доказательство
Из уравнения (30) имеем оценки:
|| fN (х) — f (х)|| = 7 + 1 = ¡7 + С&-атЦ < ||7 II + ||С || b-am < ||7 II +С Ъ-ат,
где Так как БША сходящийся, то для любого £i > 0 существует такой
N, что ||7|| = (х) — / (х) У ^ £1. Константа С не зависит от m, поэтому для любого £2 > 0 существует такая m, что Сb-am ^ £2. Последнее неравенство будет верно при m ^ |_1 — « logfe frj. Выберем £1 = а1£, £2 = а2£, где 0 < а1 ^ 1, 0 < а2 ^ 2- Для погрешности имеем оценку ||/Z (х) — f (х) || ^ Ц71| + Сb-am ^ £1 + £2 = а1£ + а2£ ^ £. □ Замечание 1. Существенными условиями для достижимости требуемой точности решений в БША являются его сходимость в точной арифметике и ограниченность нормы параметра погрешности: ||С|| ^ С, Vm ^ mmin аналогично требованию ||С1| ^ С, Vm ^ mmin для КША в [1] (определение 9, теорема 3). Именно ограниченность ЦС^ позволяет получить в БША и КША требуемую точность решения задачи ВМ.
2. Как показано в теореме 1, требуемое значение точности решения в КША = Ai+i = ||f4++i W-fwy = Ь+1+с+1 у << 1 > * р \х)= fN(i) (х)
9г = Ai = ||/4.(Х)-/(Х)Ц = Цъ+Oib-^i || ' 9г << 1 m%+1 > Jmi (х) - Jmi (х)
N (г) = S=1 Ns - число базовых операций, выполненных за г шагов (итераций) ИППМ. Очевидно, для ¿-устойчивости БША достаточно, чтобы р = Ya".^ и ^ = II|с^+||II b-(mi+1-Zi) были достаточно малы, т.е. Pi << 1 и qi << 1.
ведливы все результаты, теории, сформулированные в разделах 2 и 3, а пот,ом,у методика получения гарантированной точности решений для БША будет той же, что и для КША.
Представляет интерес изучение свойств погрешности А = Ц7 + С&-QZ|| для конкретных классов задач ВМ, решаемых БША. Учет специфики погрешности А для отдельных классов задач позволит повысить эффективность решения задач ВМ. К этим задачам относятся: нахождение суммы числового ряда; нахождение численного значения производной k-ro порядка; задача приближенного вычисления определенного интеграла; задача численного решения дифференциальных уравнений методом конечных разностей; численное решение систем нелинейных уравнений; численное решение экстремальных задач и т.д.
6. Об эффективности метода К-решений
Некоторый метод решения задачи ВМ назовем эффективным, если на данной ЭВМ решение задачи получено с заданной точностью за приемлемое время. Точность решения задается по-разному для разных классов решаемых задач и определяется некоторым признаком (условием) окончания решения задачи. Заданная точность решений часто отличается от гарантированной точности решений (ГТР), получаемой, например, в методах линейной алгебры [2]; методах решения задач, использующих интервальный анализ [10]; в ИППМ, использующей К-решения для оценки погрешностей округления.
ГТР - новое качество решений в отличие от многих методов решений, не обеспечивающих выполнение этого требования. Метод ИППМ обладает значительной универсальностью в решении задач ВМ, т.к. он не ориентирован на какие-либо классы решаемых задач. Таким образом, если метод ИППМ решает задачу за приемлемое время, а традиционный метод (ТМ), использующий «стандартное» программное обеспечение, решает, но не дает ГТР, то метод ИППМ можно считать высокоэффективным по сравнению с ТМ.
В методе ИППМ используется программное обеспечение (ПО), реализующее стандарт машинной арифметики IEEE 745 [14-16]. При выходе за диапазон длин мантисс «стандартной» арифметики (одинарной, двойной, четверной точности), наблюдается скачок увеличения времени вычислений от 10 до 100 раз в зависимости от сложности задачи и количества операций в ней. К примеру, решение СЛУ размерности К = 30 в «стандартной» арифметике с двойной точностью (га = 15 b = 10) находится, в среднем, за 5 • 10-4с, а решение той же системы при га = 16 b = 10, т.е. при выходе за пределы стандартной арифметики и «подключении» специального ПО, уже требует 3, 3 • 10-2 с. Рост времени решения при увеличении длины мантиссы - это плата за новое качество - ГТР.
7. Численный эксперимент
Решение системы линейных уравнений методом Гаусса
Рассмотрим задачу нахождения решения системы линейных уравнений (СЛУ):
Нг = с,
где Н - матрица Гильберта порядка К, т.е.
(31)
н = {hij}, i,j € [1,к] ,hij = т
г + j - 1
(32)
Алгоритм решения задачи 1 конечношаговый (КША), где функции / (х) соответствует $ = а аргументам ж соответствует матрица Н и вектор с, т.е. х = {Н, с}, г Е х Е Кп1 где п = к + т.к. матрица Н симметричная.
1
Рис. 1. Зависимость абсолютной погрешности решения Аш от длины мантиссы га
Рис. 1 представляет зависимость погрешности Аш решения СЛУ (31) от величины т при различных значениях К, Ат = \\zm — zУ, z — точное решение системы (31). Из данного графика видно, что, начиная с некоторого значения длины мантиссы ш, погрешность решения монотонно уменьшается при достаточно большом локальном увеличении длины мантиссы Аш.
а)
б)
Рис. 2. Зависимость р в области стабильности (га > 45) и области роста (га < 45) для СЛУ вида (31) от длины мантиссы т при К = 30 и шаге изменения длины мантиссы А га = 1
Как указывалось выше [11], погрешность округления носит случайный характер. На графике (рис. 2) представлено значение параметра погрешности р = \\гт — г|| Ьт при Ь =10 и известном точном решении г. Из графиков видно, что зависимость параметра погрешности р от т при т> 45 содержит элемент «случайности».
Приближенные значения мантисс чисел д^ для тех же значений %,] Е [1,10] и т^ = 110 отличались от значений мантисс д^ в 10-м или 11-м знаках, поэтому таблица для д^ не приводится.
Т а б л и ц а 1
Значение gjj'1 для СЛУ вида (31) от длин мантиссы т^ mj при К = 40
10 20 30 40 50 60 70 80 90 100
10 1 .ЗЗе+00 1.51е+00 3.09е+00 1.02е+00 2.18е-03 4.52е-13 1.99е-23 4.91е-33 4.55е-43
20 1.16е+00 2.47е+00 9.38е-01 1.63е-03 3.37е-13 1.49е-23 3.67е-33 3.40е-43
30 2.35е+00 9.73е-01 1.12е-03 2.32е-13 1.02е-23 2.52е-33 2.34е-43
40 1.19е+00 1.15е-03 2.38е-13 1.05е-23 2.59е-33 2.39е-43
50 3.23е-04 6.70е-14 2.95е-24 7.28е-34 6.75е-44
60 2.07е-10 9.14е-21 2.25е-30 2.09е-40
70 4.41е-11 1.09е-20 1.01е-30
80 2.47е-10 2.28е-20
90 9.26е-11
Таблица 1 представляет значения д3- , которые близки со значенпямп д^ при т{ ^ 60 (зона устойчивости) и отличаются от них при ^ 50.
Таблица2
Зависимость значений ^ и е^ от длины мантиссы т^ т^ при к = 40
£i 20 30 50 60 70 80 90 100
10 1.10945 1.75е+00 1.78е+00 7.36е+00 1.11е+00 1.11е+00 1.11е+00 1.11е+00 1.11е+00
20 1.48567 2.32е+00 7.98е+00 1.49е+00 1.49е+00 1.49е+00 1.49е+00 1.49е+00
30 2.1614 7.69е+00 2.16е+00 2.16е+00 2.16е+00 2.16е+00 2.16е+00
40 2.10893 6.31е+00 2.11е+00 2.11е+00 2.11е+00 2.11е+00 2.11е+00
50 7.48466 7.48е+00 7.48е+00 7.48е+00 7.48е+00 7.48е+00
60 0.00241895 2.42е-03 2.42е-03 2.42е-03 2.42е-03
70 5.01304е-13 5.01е-13 5.01е-13 5.01е-13
80 2.21032е-23 2.21е-23 2.21е-23
90 5.45175е-33 5.45е-33
Аг
Таблица 2 представляет значения относительной погрешности решений ei = jj^
.. ___А \ А 1 Гтлтт ™ . ^ ЯП т™.™™ П^^ГГ^Ч „^ГГ^^АТТТХП
^ = II/, (X
i G [1,10] и £ij = и jVsn, j > i. При mi ^ 60 имеет место хорошее совпадение их зна-
чении.
Задача нахождения производной 1-го порядка
Рассмотрим метод численного дифференцирования первого порядка:
^ ч ~ II \ Фт (х + К) - фт {х) 1 Ф(х) = Фт (х) = -К-, ф (х)
(33)
Метод численного дифференцирования относится к классу бесконечношаговых алгоритмов (БША) в том смысле, что для нахождения ф (х) с требуемой точностью надо
решить последовательность задач (33) для последовательности значений шага дифферен-
1 гл ■ // \ т фт-(x+hi (х) цирования пi ^ 0 г ^ те, а значение ф (х) = lim —1-h---ПРИ одновременном
h i ^ü
увеличении длины мантиссы т{ ^ те. Метод (33) превращается в конечношаговый при установлении зависимости между К и т. В [13] значение оптимального шага К рекомендуется брать пропорциональным \/Ё, где значения функции ф. Далее будем брать шаг = 2 • 10 2 . Задача (32) исследуется на примере функции
ф (х) = (sinx)c
х= 2/з.
(34)
Рис. 3. Зависимость параметра р(х) от m в предположении, что порядок погрешности а = 1/2
Рис. 3 представляет зависимость параметра погрешности р от т для значения порядка погрешности а = 0, 5, т.е. р = Из этого графика видно, что значение р ^ 2 при
10 < т ^ 300.
ТаблицаЗ
Значения решений iE [1,10]
(х)
10 8.8828000000е-01
20 8.88290852800000000000е-01
30 8.882908528554890000000000000000е-01
40 8.8829085285548970217500000000000000000000е-01
50 8.88290852855489702189867500000000000000000000000000е-01
60 8.882908528554897021898676192135000000000000000000000000000000е-01
70 8.8829085285548970218986761921499853000000000000000000. . . 0е-01
80 8.8829085285548970218986761921499854278385000000000000. . . 00000000000е-01
90 8.8829085285548970218986761921499854278396686450000000. . . 000000000000000000000е-01
100 8.8829085285548970218986761921499854278396686552550000. . . 0000000000000000000000000000000е-01
Таблица 3 представляет значения решений (х) = (х) при т\+1 = т\ + ¿Дт, т\ = 10 Дт = 10 1 £ [1, 9]. Как видно из таблицы, числа ^ СПЗ для решений (х) равны: ¿1 = 4, £2 = 10 ¿3 = 15 ¿4 = 19 Ь = 24, £6 = 29, £7 = 34 £8 = 39, ¿9 = 44 Для оценки числа ^ю необходимо вычислить значение фт11 (х). Значенпя ^ СПЗ выделены жирным шрифтом. Эта таблица иллюстрирует большую практическую важность правила СПЗ как
метода нахождения решения задачи ВМ (решение ftl (х) указано в явном виде), так и его погрешности, равной Aj = ||f (х) — fli (х)|| ^ 1 ■ 10е-ti.
8. Заключение
В настоящей работе предложен метод численного анализа погрешностей округления решения задач ВМ. Результаты, полученные в статье, сводятся к следующим положениям:
1) Введено понятие КР задачи, которое в оценках погрешностей решений с некото-
( х)
буемой гарантированной точностью и получены оценки погрешности, далее численно реализуемые в 111111 \!.
2) Предложены алгоритмы, позволяющие оптимизировать процесс решения задачи в 111111 \!. Рассмотрены методы округления полученных решений, причем округленное решение имеет гарантированную точность.
3) Доказана теорема об оценках погрешности метода (правила) округления решения по
шает значения е = 10e-i, где е - порядок числа.
4) Для бесконечношаговых алгоритмов (БША) решения задач ВМ доказана теорема о достижимости требуемой точности решения.
5) Предложенный метод КР оценки погрешностей округления обладает следующими свойствами:
б) Метод КР обладает универсальностью в том смысле, что он не ориентирован на решение конкретных классов задач ВМ.
6) Приведены результаты численного эксперимента, иллюстрирующие основные свойства метода КР оценки погрешности решений.
Актуальность предложенного метода КР в первую очередь заключается в возможности получения ГТР для различных классов задач ВМ. Метод имеет перспективы развития в том смысле, что определение границ его применимости и численной эффективности для различных классов задач открывает новую область исследований в вычислительной математике.
Литература
1. Бирюков А.Г., Грипевич А.И. О гарантированной точности решений задач вычислительной математики в арифметике с плавающей запятой и переменной длиной мантиссы // Труды МФТИ. - 2012. - Т. 4, № 3. - С. 171 ISO.
2. Годунов С.К., Антонов А.Г., Кирилюк О.П., Костин В.И. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. - Новосибирск: Наука. Сиб. Отд-ние, 1988. - 456 с. ISBN 5-02-028593-5.
3. Wilkinson J.H. Rounding Errors in algebraic processes. - Englewood Cliffs, N.J.: Prentice-Hall, 1963. ISBN 0-486-67999-3.
4. Higham N. J. Accuracy and stability of numerical algorithms. - Philadelphia : Society for Industrial and Applied Mathematics, 1996.
5. Воеводин В.В. Вычислительные основы линейной алгебры. - М.: Наука, 1977. - 304 с.
6. Henrici P. Elements of Numerical Analysis. - New York. - John Wiley к Sons Inc., 1964.
7. Clenshaw C. W. and Olver F. W. J. Beyond floating point //J. Assoc. Comput. Mach. -1984. - V 31. - P. 319-328.
8. Langlois P. A Revised Presentation of the CENA Method. - ARENAIRE - INRIA Grenoble Rhcjme-Alpes / LIP Laboratoire de l'lnformatique du Parallelisme.
9. Шокип Ю. И. Интервальный анализ. - Новосибирск: Сибирское отд. изд-ва «Наука», 1981.
10. Алефельд Г., Херцбергер Ю.. Введение в интервальные вычисления. - М.: Мир, 1987.
11. Воеводин В.В. Ошибки округления и устойчивость в прямых методах линейной алгебры. - М.: Изд-во МГУ, 1969. - 140 с.
12. Демидович Б.П., Марон И.А. Основы вычислительной математики. - СПб.: Лань, 2009. ISBN 978-5-8114-0695-1.
13. Бахвалов Н.С., Жидков Н.П., Кобельков P.M. Численные методы. - М.: Наука, 1987.
14. IEEE 754-2008: 754-2008 IEEE Standard for Floating-Point Arithmetic. - ISBN: 978-07381-5753-5.
15. GNU GMP: Multiple precision arithmetic library / http://gmplib.org/
16. GNU MPFR, http://www.mpfr.org/
17. Математическая энциклопедия Т. 4 - M.: Советская энциклопедия, 1984.
Поступила в редакцию 13.01.2013.
