Демонстрационные возможности примера бабушкивитасека-прагера в точечных и интервальных расчетах Текст научной статьи по специальности «Философия, этика, религиоведение»

2005 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 10 Вып. 2

ИНФОРМАТИКА

УДК 519.6

Г. Г. Меньшиков

ДЕМОНСТРАЦИОННЫЕ ВОЗМОЖНОСТИ ПРИМЕРА БАБУШКИ-ВИТАСЕКА-ПРАГЕРА В ТОЧЕЧНЫХ И ИНТЕРВАЛЬНЫХ РАСЧЕТАХ

1. Введение. Точечное исполнение примера. Этой статьей мы обращаем внимание читателя на богатые демонстративные возможности примера, предложенного в известной книге [1]. Покажем, что сами по себе интервальные (даже доказательные, «локализующие») расчеты еще не несут в себе гарантии высокого качества вычислений. Все дело в устойчивости «внешнего алгоритма». Он негоден - все рушится. Рассмотрим пример нахождения чисел

1 Г1

1п = - хпех<1х, п = 0, 1,... . (1)

е ./о

Ясно, что 10 = I — е-1. Выразим 1п через 1п-\- Интегрируя по частям, получаем рекуррентное равенство

1п = 1 -п/п-ь (2)

Если проводить вычисление Д, /2,... по формуле (2) на реальном компьютере, т.е. с учетом округлений, то уже при п порядка 10-15 появляются заведомо неправдоподобные (отрицательные) значения 1п (см. столбцы 1, 2 табл. 1).

Таблица 1. Значения /т

1 2 3

п (1п) (1п)

0 .632 120 5 .632 120 558 828 557 7

1 ;367 879 5 .357 879 441 171 442 3

16 .-191 438 3 5.545 930 172 957 014Б-02

17 .325 445 3 5.719 187 059 730 757Б-02

18 -5.858 015Е+07 -2.945 367 075 153 627Б-02

19 1.113 023Е+09 1.559 619 744 279 189

20 -2.226 046Е-Ы6 -30.192 394 885 583 78

Выясним, в чем тут дело.

Очевидно, все 1п > 0. Они убывают, так как

1п = - [ х • хп~хехйх < - ( хп~гехс1х = /п_ь е Уо е к

Полученное двойное неравенство

0</п</п_1 • (3)

© Г. Г. Меньшиков, 2005

С^В-строки Комментарии

1200 в 2 Инициализация

1210 Б=0 Зануление п

1220 А=-1: С 86 -1 пара 0

1230 А=0: С 16 е~1 пара 0

1240 в 9 —е-1 =>• пара 0

1250 А=1: в 80 1 — е"1 =>■ пара 0

1260 РИШТ "п = ";Д" [1п] = "; Печать: п

1270 А=0: С 50.1 Печать: [1п]

1280 Б=Б+1 Увеличение п на 1

1290 А=Б: в 82 (п + 1) [7П] => пара 0

1300 в 9 -(п + 1) [/„] => пара 0

1310 А=1: С 80 [/п+1] = 1 - (п + 1) [1п] пара 0

1320 вОТО 1260 Переход к строке 1260

применим к формуле (2). Тогда имеем новое неравенство

0< 1-п/п_1 </п_ь (4)

которое решаем относительно 1п-\\

1 г 1

<1п-!<-• (5)

Теперь из (1) видно, что Тогда из (2)

. 1 ^ -УЬ— X ^

п + 1 п

/п - - 0 {п -> +оо). (6)

п

1п = 1 -п/п-1 - 1 - ^ (п-> +00). (7)

Возникает вычитание близких величин. Это объясняет кажущуюся аномалию с 1п.

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

2. Интервально-локализующее исполнение примера ВВП. Повторим тот же пример, на этот раз в интервально-локализующей манере. Следуя второй теореме о композициях [2], запишем правую часть рекуррентного соотношения в отрезках-локализаторах. Получим включение

/пЕ1-п[1п_1] (п = 0,1,...). . (8)

В нем под [/п_1] имеется в виду находимый машиной интервал-локализатор для 1п-1-Так как правая часть по упомянутой второй теореме является локализатором для /п, то можно было бы ее так и обозначить. Тогда, принимая во внимание мажоризацию (т.е. вспомогательное расширение локализатора при выполнении стандартных операций, связанных с появлением погрешностей), запишем дальнейшую разновидность соотношения (3):

[/п]51-п[/пч] (п = 0,1,...). (8а)

Соответствующая программа оформлена табл. 2. При этом взят язык С^Вавк: за основу интервального ассемблера (как в §17, 18 книги [2]).

Результаты же нескольких шагов интервального расчета приведены в столбцах 1-3 табл. 3.

С}В-строки Комментарии

1200 в 2 Инициализация

1210 Б=0 Зануление п

1220 А=-1: в 86 -1 => пара 0

1230 А=0: С 16 е~1 => пара 0

1240 С 9 —е-1 пара 0

1250 А=1: в 80 1 — е~1 пара 0

1260 РИШТ "п = ";£)," [/»! = »; Печать: п

1270 А=0: в 50.1 Печать: [7П]

1280 Увеличение п на 1

1290 А=П: в 82 (п + 1) [/п] пара 0

1300 С 9 —(п + 1) [1п] пара 0

1310 А=1: в 80 [/п+х] = 1 - (п + 1) [1п] пара 0

1320 вОТО 1260 Переход к строке 1260

применим к формуле (2). Тогда имеем новое неравенство

0< \-nIn-x <1п-и (4)

которое решаем относительно 1п-\:

1 г 1

</»-!<-. (5)

п + 1 п

Теперь из (1) видно, что

. /п ~ - 0 (п +оо). (6)

Тогда из (2)

1п = 1 -п/п-1 - 1 - тт1 (п-> +оо). (7)

1/п

Возникает вычитание близких величин. Это объясняет кажущуюся аномалию с 1п.

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

2. Интервально-локализующее исполнение примера ВВП. Повторим тот же пример, на этот раз в интервально-локализующей манере. Следуя второй теореме о композициях [2], запишем правую часть рекуррентного соотношения в отрезках-локализаторах. Получим включение

1п € 1-п[1„_1] (п = 0,1,...). (8)

В нем под [1п-1] имеется в виду находимый машиной интервал-локализатор для /п_1. Так как правая часть по упомянутой второй теореме является локализатором для /п, то можно было бы ее так и обозначить. Тогда, принимая во внимание мажоризацию (т.е. вспомогательное расширение локализатора при выполнении стандартных операций, связанных с появлением погрешностей), запишем дальнейшую разновидность соотношения (3):

[1п]2 1-п[/„_1] (П = 0,1,...). (8а)

Соответствующая программа оформлена табл. 2. При этом взят язык С^Вавк за основу интервального ассемблера (как в §17, 18 книги [2]).

Результаты же нескольких шагов интервального расчета приведены в столбцах 1-3 табл. 3.

1 2 3 4 5

п ['»] ['»] Ш([/п])

0 .632 120 3 .632 120 8 4.77Е-07 .632 120 3 .632 120 8 4.77Е-07

1 .367 878 9 .367 88 1.01Е-06 .367 878 9 .367 88 1.01Е-06

9 -.169 765 5 .359 630 4 5.29Е-01 9.090 908Е-02 .1 9.09Е-03

15 -941 903 965 832.4 1.91Е+06 5.822 351Е-02 6.250 001Е-02 3.62Е-03

Далее, при возрастании п ширина локализатора растет, причем с ускорением. Наконец, уже при п = 9 можно констатировать, что дальнейший расчет беспредметен из-за резкого падения точности: отношение ширины [1п] (обозначаемой как и) (/п)) к /п имеет порядок 0,1. Тем более велика и) ([/15]) •

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

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

Для экспериментального снятия мажоризации можно положить в базовой интервальной программной системе С — 0. Тогда результаты работы той же программы будут в форме вырожденных отрезков (табл. 4). Поэтому каждый результат даем не в форме отрезка, а одним числом. Как видно, качественно получается тот же эффект, что в точечном исполнении, и можно полагать, что интервальное, но не обязательно локализующее исполнение примера не принесет существенного сужения ширины локализатора. Ему также может быть свойственна низкая точность.

4. Локализующее усовершенствование примера

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

С одной стороны, по интервальному аналогу рекуррентной формулы

[/п]31-п[/п-1] (9)

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

С .другой стороны, теоретическим путем в [2] доказана ограниченность {1пу - полу-

Таблица 4• Результаты расчета в форме вырожденных отрезков

п

0 .632 120 5

1 .367 879 5

14 -797. 597 3

15 -119 64.96

чено неравенство

1 г 1

-—г </„<-• п + 1 п

Соответствующий отрезок - не что иное как еще одна система локализаторов:

1 1

п + 2' п+ 1

С их помощью можно «унять» катастрофический рост ширины. В самом деле, пересечение интервальных локализаторов есть интервальный локализатор:

1 1 '

п + 2' п + 1

Как видно (в идеальной модели интервальных вычислений [3]),

С (1-п [/„_!]) П

(10)

(п + 1)(п + 2)'

(П)

Отсюда следует, что ю (1п) ->• 0.

Посредством этой усовершенствованной интервальной алгоритмики проведем еще одну серию численных экспериментов (табл. 5). Затем проделаем вычисления при тех же п, что ранее. Результаты содержатся в 1, 4, 5 столбцах табл. 3. На этот раз ширина не только не возрастает, но, более того, имеет тенденцию к относительно сколь угодно малым значениям.

Более того, относительная ширина также является бесконечно малой:

*>{%)) < и>{[1п]) <

1

0.

(12)

. Итак, усовершенствование интервального подхода позволяет экспериментально исследовать точность и управлять ею.

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

Практическая ценность предложенного интервально-локализующего усовершенствования состоит в том, что оно не только устраняет катастрофический исход вычислений (это делают и две модификации примера БВП [см. 1, с. 31 и 32]), но и приносит новое качество вычислительной работы - гарантию точности расчетов и, более того, позволяет управлять точностью за счет подбора (посредством аналитических рассуждений) более узкого локализатора, которым путем пересечения формируется еще более узкий. Такой гарантирующий механизм отсутствует в упомянутых модификациях примере БВП, поскольку они используют традиционный, неинтервальный подход к организации вычислений.

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

Заметим, что в табл. 2 и 5 аббревиатура в обозначает оператор ООЭиВ.

Вообще, использование столь «древнего» языка, как С^Вавк;, может вызвать недоумение. Причина нашего внимания к (^Вавк^'у в том, что до сих пор не обнаружены арифметические аномалии, в этом языке, т. е. случаи нарушения гипотез о машинной

QB-строки

1200 G 2

1210 D=0

1220 А— -1: G 86

1230 А=0: G 16

1240 G 9

1250 А=1: G 80 1 - е-1 =»

1260 PRINT "п = ";Д" [7П]

1270 А=0: G 50.1

А= —D: G 40

1280 D=D+1

А= —D: G 82

А=1: G 80

А=1: G 40

А=1: G 86

A=D=2: G 83

А=2: G 40

А=1: G 86

A=D+1: G 83

А=2: G 61

А=1: G 34

Комментарии Инициализация Зануление г

1=> пара О е-1 => пара О е-1 пара О пара О Печать: п Печать: [1п] [/»] =>пара 1 Увеличение п на 1 —(п -1-1) [1п] => пара О [7п+1] = 1 - (п + 1) [7П] пара О 1 — (п + 1) [1п] =>• пара 1 [1, 1] =» пара О

* пара О

п + 2 1

[1

п + 2 1] = 1

71 + 1

пара 2 пара О пара О

1 1

п + 2' п+1 / • //

GOTO 1260

1 ^ J [п + 2' 72 + 1 _

[Jn+l] = [Jn+l]' П [Jn+l]" Переход к строке 1260

пара 0 пара 0

арифметике, перечисленных в §11 книги [2] и обеспечивающих доказуемость результатов интервально-локализуюгцих вычислений. Такие случаи не были обнаружены, несмотря на то, что в составе интервального ассемблера этот программный продукт (QBasic) используется в студенческих лабораторных работах еще с 1997 года. Таким образом, в упомянутом смысле это - высоконадежный язык.

Summary

Men'shikov G. G. Show oportunities of Babushka-Prager-Vitasek example in the point and interval computation.

We show on this well-known example that the interval computations (even the proving, localizing ones) themselves does not yet guarantee a high quality of the computational work. The matter is the stability of the "outer" algorithm. If it is unfit then the whole fails.

Литература

1. Бабушка ИВитасек Э., Прагер М. Численные процессы решения дифференциальных уравнений / Пер. с англ. В. JI. Каткова; Под ред. Г. И. Марчука. М., 1969. 368 с.

2. Меньшиков Г. Г. Локализующие вычисления: Конспект лекций. Вып. 1. Введение в интервально-локализующую организацию вычислений. СПб., 2003. 89 с.

3. Меньшиков Г. Г. Локализующие вычисления: Конспект лекций. Вып. 2. Задачи композиционного расчета и проблема грубости их интервально-локализующего решения. СПб., 2003. 59 с.

Статья поступила в редакцию 21 апреля 2005 г.