Главная Промышленная автоматика.

1962, № 3. Бухнер К. X. «Подпверждение к алгоритму 60» «САСМ», 1962, № 5) и показал себя одним из удобных и надежных методов интегрирования. Поэтому имеет смысл указать, что алгоритм не вполне совершенен и что существует -значительный класс подынтегральных функций, для которых экстраполяционные значения являются худшими оценками интеграла, чем соответствующие суммы формулы трапеций.

Точность процедуры Ромберга зависит от возможности разложения погрешности метода трапеций по степеням /г где h - величина шага. Разложением такого типа является ряд Эйлера -Маклорена. Другое выражение может быть получено с помощью рядов Фурье. Коэффициенты при /гг в формуле Эйлера- Маклорена пропорциональны разности значений (2г+1)-й производной на двух концах интервала. Следовательно, любой интеграл, для которого нечетные производные подынтегральной функции на концах интервала либо равны нулю, либо равны между собой, не сходится при экстраполяции по Ромбергу. Простыми примерами таких интегралов являются интегралы периодических функций на интервале более длинном, чем период функции, и интегралы, для которых производные равны нулю на обоих концах интервала. Примером последнего типа служит интегральная аппроксимация модифицированной функции Ханкеля [341]

еКр {х) = J e*c-coftO cosh (pt) dt.

где L берется таким большим, что частью интеграла от L до оо можно пренебречь. Несколько других примеров даются под заголовком «Особые случаи» Бауэром, Рутисхаузером и Штифелем [151]. Эта статья - наиболее обширное исследование метода Ромберга на английском языке.

Алгоритм непригоден также и тогда, когда разложение остаточного члена содержит другие степени h или хотя бы одну такую степень. Рутисхаузер [161] рассматривает оценочные интегралы в форме

y{x)dx=f(<P(x)!Vx)

Если такие интегралы вычисляются по формуле трапеций в предположении, что/(0)=0, то погрешность может быть выражена в форме Ilckh+YhZdkh. Хотя стандартная экстраполяция Ромберга непригодна для этой последовательности оценок, Рутисхаузер дает модифицированную процедуру, которая является эффективной. Экстраполяция непригодна также тогда, когда подынтегральная функция разрывна, хотя с точки зрения вычислительной математики это исключение тривиально.

Было также указано (МсКеетап W. М. Личная переписка. Сент., 1963; Engeli М. Личная переписка. Янв., 1964), что процедура Ромберга может накапливать ошибки округления. Для большинства случаев применения потери ввиду их значительности предотвратить невозможно.



АЛГОРИТМ 616

Процедуры интервальной арифметики [А1]

Термин «интервальное число» был использован П. С. Двайером. .Машинные процедуры интервальной арифметики разработал примерно Б 1958 г. Р. Мор tli].

Если ахЬ и cyd, то процедура rangsum (сокращение от range - интервальная, sum -сум1ма) дает интервал [e,f] такой, что 4x+yf*. Из-за особенностей машинных операций (округления и усечения) .машинные суммы а+с и b+d могут не обеспечивать сохранности концевых точек выходного интервала. Поэтому для процедуры rangsum нужна нелокализованная вещественная процедура correcsum (сокращение от correction - поправка и sum - сумма), которая ko,m-пенсирует погрешности машинной арифметики. Тело процедуры correcsum будет зависеть от типа маШИиы, для которой эта процедура написана, поэтому она здесь не приводится (однако ниже дается пример такой процедуры). Предполагается, что процедура correcsum имеет вещественные параметры v и и и целый парзметр i и что ей сопутствует иелокализованная вещественная процедура correction, которая дает верхнюю границу величины погрешности, содержащейся в машинном представлении числа. Выходное значение функции correcsum дает левый конец интервала, выдаваемого процедурой rangsum, когда correcsum вызывается с параметром i=-1, и правый конец, когда i=l.

Процедуры rangsub, rangpro и rangdiv (сокращение от subtraction- вычитание, product - произведение, division - шлепш) предназначены для остальных основных операций интервальной арифметики. Процедура rangsqr (сокращение от square - квадрат) дает интервал, внутри которого должен лежать квадрат интервального числа.

Процедуры rangsumc, rangsubc, rangproc и rangdivc предусмотре-ны для операций с комплексньши интервальными аргументами, т. е. вещественные и мнимые части этих аргументов являются интервальными Числами.

procedure rangsum (a,b,c,d) result: (e,f);

value a,b,c,d; real a,b,c,d,e,f; begin e: = correcsum(a,c,- 1); f:=correcsum(b,d,l) end rangsum;

Процедура rangsub вычисляет разность интервальных чисел {с,&}и {c,d}, т. е. находит интервал [e,f] такой, что ех-yf для всех х и у, удовлетворяющих условиям ахЬ и cyd. В теле процедуры rangsub процедура rangsum не локализована.

procedure rangsub (a,b,c,d) result: (e,f);

value a,b,c,d; real a,b,c,d,e,f; rangsum (a,b,-d,-c,e,f);

Процедура rangpro находит произведение {e,f} интервальных чисел {a,b} и {c,d}, т .е. находит интервал [e,f] такой, что exXyf для всех X и у, удовлетворяющих условиям ахЬ и C:yd. В теле прОЦеду-ры rangpro имеется глобальная процедура correcpro, аналогичная вышеописанной процедуре correcsum (см. нижеследующий пример).

* Другими словами это можно выразить так, что процедура rangsum находит сумму {е. f} двух интервальных чисел {а, Ь} и {с, d). {Прим. ред.)



procedure rangpro (a,b,c,d) result: (e,f);

value a,b,c,d; real a,b,c,d,e,f; begin real r;

e: = i: = aXc;

for r: = aXd,bXc,bXd do if r<e then e: = r else if r>f then f:=r;

e: = correcpro (e,-1);

f:=correcpro(f,l) end rangpro;

Процедура rangdiv находит частное от деления интервального числа {а,Ь} на интервальное число {c,d}, т е находит интервал [e,f] такой, что enx/yf для всех хну, удовлетворяющих условиям ахЬ и cyd. В теле процедуры rangdiv используется глобальная вещественная процедура correcdiv, аналогичная (возможно идентичная) процедуре correcpro. Если интервальный делитель содержит нуль, то осуществляется выход из программы « нелокализованной метке signalGl, procedure rangdiv (a,b,c,d) result: (e,f);

value a,b,c,d; real a,b,c,d,e,f; begin if cO j\dO then go to signaI61; e: = f: = a/c; for r: = a/d,b/c,b/d do if r<e then e:=r else

if r>f then i: = r; e:=correcdiv (e,-1); i: = correcdiv(f, 1) end rangdiv;

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

procedure rangdiv2 (a,b,c,d) result: (e,f);

value a,b,c,d; real a,b,c,d,e,f; begin if cO Д dO then go to signal61;

rangpro (a,b, 1 /d, 1 /c,e,f);

e: = correcdiv (e,-1); f: = correcdiv (f,l) end rangdiv2;

Процедура rangsqr находит квадрат интервального числа {а,6}, Т. е. находит интервал [е,/], такой что exf для всех х, удовлетворяющих условию ах,Ь. В теле этой процедуры нспользуется нелокализованная процедура correcpro.

procedure rangsqr (a,b) result: (e,f);

value a,b; real a,b,e,f; begin e: = aXa; f:=bxb; if a<0 then

begin e:= if b<0 then f else 0;

if-a>b then f: = aXa end;

e:=correcpro (e,- 1); f: = correcpro (f, 1) end rangsqr;

Процедура rangsumc находит сумму комплексных интервальных чисел {а\,а2,Ь\,Ь2} и {c\,c2,d\,d2}, т. е. находит интервалы [el,e2] и





0 1 2 3 4 5 [6] 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

0.0017