Быстрый алгоритм расчета интегралов специального вида Текст научной статьи по специальности «Математика»

быстрый алгоритм расчета интегралов специального вида

В.А. Куделькин, Ю.Л. Ратис Консорциум «Интегра-С» Самарский государственный аэрокосмический университет Институт систем обработки изображений РАН

Аннотация

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

Введение

В работах [1, 2] было показано, что проблема нахождения функции отклика оптикоэлектронного датчика с учетом дифракционных поправок сводится к задаче вычисления интегралов специального вида Л, (р, д) и Лс (р, д).

Л, (р, д) = [—Бт( р( 2)соБ(д/), t

о 1

Лс (р, д) = [—со8^ 2^т(дО. t

(1)

(2)

Как показал последующий анализ, интегралы (1), (2) в общем случае не могут быть представлены в виде комбинации конечного числа элементарных и/или специальных функций. С формальной точки зрения численные значения интегралов специального вида Л, (р, д) и Лс (р, д) могут быть найдены на

основе общеизвестных квадратурных формул. Однако фактически, вычисление этих интегралов по квадратурным формулам Симпсона, Гаусса и т. п. оправдано только при небольших значениях аргументов р и д .

В тех случаях, когда выполняются условия р »1 и/или д »1, подынтегральные выражения в

формулах (1), (2) становятся быстроосциллирую-щими, и использование квадратурных формул для вычисления величин Л, (р, д) и Лс (р, д) становится некорректным, т. к. погрешность вычисления быстро растет с увеличением р и/или д .

В конечном счете, задача корректного расчета величин Л, (р, д) и Лс (р, д) численными методами сводится к необходимости избавиться от сложения большого числа знакопеременных слагаемых, имеющих одинаковый порядок величины. С точки зрения теории аналитических функций, общий подход к преодолению указанной проблемы основан на использовании идеи метода перевала [3]. Однако технически эта идея может быть реализована несколькими способами. А для использования датчиковой аппаратуры в реальном масштабе времени необходимы быстрые и устойчивые алгоритмы обработки сигналов, т. е., функции отклика.

В настоящей работе проделано сопоставление вычислительной эффективности различных алго-

ритмов расчета интегралов (1), (2) с целью выявления наиболее быстрого и устойчивого из них.

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

1. Случай р ^ 1 и д ^ 1

В соответствии с общей идеологией численных методов в этом случае алгоритм расчета интегралов специального вида строится на основе степенных рядов.

Очевидно, что подынтегральные выражения в формулах (1) и (2) имеют хорошие аналитические свойства, и в соответствии с этим величина Л, (р, д) может быть представлена в виде:

Л, (р, д) = [—sin(pt2)cos(gt) = t

=ё

(-1)У

(2п)!

БШ( pt .

(3)

Произведем замену переменных х = t2. В результате приходим к выражению:

4 (р, д) = р) +

1 ™ (_ 1)ид2и 1 +1Ё((12) д 15т(рх)хп_хйх,

(4)

2 И=1 (2п)! о

где 81(х) - интегральный синус, а интеграл в фор муле (4) является табличным

1 (_1)кр2к+1

| 8Ш( рх) х" Хйх = Ё

к=0( п + 2 к +1)(2 к +1)!

(5)

Из (3) и (5) следует, что интеграл Л, (р, д) может быть представлен в виде двойного степенного ряда:

"=0

А (р, д) = = 2 ^ (-1)пд2

-х-

(-1)кр2к+1

и=0 (2п)! к=0(и + 2к + 1)(2к +1)!

= р) + +1 х (-1)"д2

(6)

-х-

(-1)кр2к+1

2 (2п)! к=0(и + 2к + 1)(2к +1)! Совершенно аналогично вычисляется интеграл

А (р, д): Ас (р,д) =

(-1)кр2к =

= хх (-1)"д2я+1 х

п=0 (2п +1)! к=0(4к + 2п + 1)(2к)! Г п,—л

= д,

х

• с

(-1)пд2

(7)

х

(-1)кр2к

=1 (2п +1)! й(4к + 2п +1)(2к)!

В формулах (6) и (7) в явном виде выделены главные (нулевые) члены разложения, для которых имеются стандартные алгоритмы расчета.

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

2. Случай р ~5 и д ~5

В случае умеренных значений аргументов (р, д ~ 5) степенные ряды (6) и (7) сходятся существенно медленнее, чем при р « 1, д « 1. Поэтому для указанного диапазона аргументов вычисление интегралов (1) и (2) производилось двумя способами 1) разложением в степенной ряд; 2) по квадратурным формулам. Прямым расчетом установлено, что при р, д ~ 5 наиболее эффективной является процедура

вычисления интегралов Л! (р, д) и Лс (р, д) по квадратурным формулам Гаусса по 20 х к точкам, где к = тах {[р], [д]} . Относительная погрешность вычислений е при этом не превосходит е < 10-5, что вполне приемлемо для практических расчетов. Контроль вычислений с помощью квадратурной формулы Симпсона по 500 точкам подтвердил этот вывод.

3. Случай умеренных р < 5 и произвольных д > 5

В соответствии изложенной выше общей идеологией построение эффективного алгоритма расчета интегралов (1) и (2) сводится к задаче их представ-

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

Представим выражение для интеграла Лх (р, д) в следующем виде:

А (р, д) = Яе

[—Бт( рг2) ехр(/'дг) t

= Яе

XI" (2п +1) ]„ (д)фп (р)

(8)

При переходе от (1) к (3) мы воспользовались соотношениями [4]:

соз(уг) = Яе (ехр(/>/)) . (9)

ехр(/у0 = X ? (2" +1)]п (у)Р„ (г). (10)

п=0

Из (8) легко видеть, что

да

а ( р, д) = X (-1)п (4 п+1) кп (д) ф2п (р), (11)

где

1 ¿г 2

Ф2п ( р) = | — ®1П( 2) Рр2п Ц ) = 0 *

да (-Пк р2к+1 1

=х( 1) р г ¿г • г4к+' • р2 (*). (2к +1)! Г 2п

(12)

Содержащийся в формуле (12) интеграл является табличным [4, 5]

Г ¿г • гАк+1 • Р2п (г) =

0

= (- 1)п Г(п - 2к -1/ 2)Г(2к +1) = 2Г(-2к - 1/2)Г(п + 2к + 2) ;

(13)

откуда вытекает окончательное выражение для коэффициентов разложения (12):

ф2 (р)=А X(-1)^р2к+1 • (4к+1)п .

2п 2 2к +1 (п + 2к +1)!

(14)

где (а)п - символ Похгаммера (см. Приложение). Совершенно аналогично вычисляется интеграл

лс (р, д):

да

Ас (р, д) = X (-1)п (4п+3) • у2я+1 (д) • Ф2„+, (р), (15)

где

1 ¿г

Ф2п+1(р) - Г — С0«(рг2) • р2п+1(г): -1 г

( 1\к „2к 1 = х 1 ¿г • г4к-1 • р!п+1 (г).

(16)

(2к)!

и

п=0

п=0

Коэффициент ф2и+1(p) также выражается через табличный интеграл [4, 5]:

( - а

\dt • tAk-1 • R

i(t) =

2 (2k +1/2 )n

(17)

Таким образом, окончательное выражение для коэффициентов ф2п+1(р) имеет вид:

Wp) 2 ^ (2^ (^k ч-1 / 2))n

(18)

Соотношения (11), (14), (15) и (18) полностью решают задачу вычисления интегралов (1), (2) в случае умеренных значений аргумента p < 5 при произвольных значениях аргумента q. В дополнение к сказанному отметим, что ряды (11) и (15) сходятся быстро, а сферические функции Бесселя вычисляются с помощью вычислительно эффективного и устойчивого алгоритма, построенного и исследованного в работах [6, 7, 8].

4. Случай умеренных q < 5 и произвольных p > 5

Описанный в предыдущем разделе алгоритм представляется весьма изящным, однако ряды (11) и (15) имеют достаточно сложную структуру коэффициентов разложения и содержат специальные функции математической физики (сферические функции Бесселя). Все вышесказанное ориентирует на поиск гибридных алгоритмов, в рамках которых изначально выделены «главные члены» разложения интегралов специального вида As (p, q) и Ac (p, q) по той или иной комбинации параметров.

Из соотношения (4) следует, что

Ах (p, q) = 2si(p) +

1 ™

+2 ^

z n=1

(-1)V

(2n)l

-Im J exp(ipx) xn 1dx.

(19)

Выражение (19) удобно представить в виде:

Ах (p, q) = |si(p) +

1 » (—1)nq2" дn-1

+—Z -——— Im-

2 ti (2n)l d(ip)n

1J exp(ipx)dx

(20)

Откуда немедленно следует, что

As (p, q) = |si(p) - ^ >

dm ^

(—1)nq2

дn

о (2n + 2)! d(ip)n

exp(ip) -1 ip

(21)

Производная n - го порядка в (21) вычисляется элементарно:

дn

d(ip)"

дn

exp(ip) -1 ip

[exp(ip) • (ip)-1 ]-

дт

откуда следует, что:

дn (ip)'1 д(ip)n

(22)

дn

дт-

exp(ip) -1

ip

= exp(ip)Z

n )дr (ip)-1 дn (ip)-

(23)

=0 ^r ) д^У д(/p)и и, таким образом:

дn

д(ip)n

exp(ip) -1

ip

= (-1)nnl(ip)

-(n+1) .

exp(ip)X

(-1)k (ip)k k l

-1

(24)

Конечная сумма, фигурирующая в квадратных скобках в формуле (24) есть не что иное, как частичная сумма ряда Тейлора для экспоненты (см. Приложение):

дn

exp(ip) -1

/p .

д(/p)и

= (-1)nnl(ip)-(n41) [exp(ip) • е (-ip) -1].

(25)

Учитывая это обстоятельство, преобразуем выражение для Л, (р, д) к следующему виду:

1 L-, ч q cosp As (P, q)=2 •|si( p)++

+ Im [ S (p, q) - exp(ip) x

Z nlq2n+2(ip)-(n41) ( . ) xZ-(2 + 2)l--en (-/p)

n=1 (2n + 2)l

(26)

где

Z nlq2n+2(ip)-(n41)

s ( p, q) -Z q(2 +p))l n=0 (2n + 2)l

(27)

Для того чтобы выразить мнимую часть степенного ряда в формуле (27)

Sta(p, q) = Im {S(p, q)} ,

(28)

через специальные функции математической физики, поступим следующим образом. Перепишем (28) в виде:

SIm (P, q) = -Z

(2k)l(-1)k I q

to (4k + 2)l ^ p

(29)

Введем переменную z = q • p 1 2 и исследуем свойства функции:

S (1z) -Jt (2k)l(-1)kz«+ 2

SIm(1, Z) = Z (4k + 2)l Z .

Из (30) видно, что

¿2

= -л/Л^,

=-х

(2к )!(-1)к м (4к +1)!

-2 =

(31)

где Пу (2, - функция Ломмеля двух переменных (см. Приложение). Таким образом, получаем:

^1т(1, 2 ) = - ^Х

<Ит

(-1)к Г 22

X ,

к=04к + 2 | 2£

J 2

Д)

(32)

Существует и другой способ нахождения функции £ (р, д). Для того чтобы просуммировать исследуемый ряд (27), воспользуемся соотношением:

х2 = 1 + ^ х ехр егг I-к=0(2к)! 2 Н 4 1 12,

(33)

и проинтегрируем его в пределах от 0 до 2 . В результате, мы приходим к соотношению:

(X — X 2кс1х = [ к=0(2к)!

= 1

1 х ехр (т 1ег£" 1|

(34)

• ¿х,

из которого вытекает, что:

X-

к!

к=0(2к +1)! и, таким образом: к!

= ^ехр| хр!^I |

X"

22к+2 =л/П| ¿х ехр I — | етГ

(35)

(36)

к=0(2к + 2)!

Выражение (36) можно также представить, как к!

К ■ 2к+2

-2 =

)!

2/2 , (37)

X

к=0(2к + 2)!

2/2 :

= егГ1 (2/2)• егГ(2/2)-4 | ¿у• ^(у)

где ^ (2) - интеграл Досона (см. Приложение).

Вычислим интеграл Лс (р, д) для исследуемого случая умеренных д < 5 и произвольных р > 5 . С учетом проделанных выше вычислений легко показать, что:

Лс (р,д) =

= 1ке 2

=0 (2п +1)!I /р

I У(п +1/2,ip)

(38)

Воспользуемся известным свойством гамма-функций:

(2п +1)! = Г(2(п + 1)) = = (2л)-1/2 22п+3/2 Г(п + 1)Г(п + 3/2). Отсюда следует, что Лс (р, д) - 1

(39)

х Яе

X

^л/п

(-1)пу(п +1/2,1р) Г д

(40)

=0 п!Г(п + 3/2) 14р

Из теории гамма-функций известно также, что у (п + 1/2, 2 2)

Г(п +1/2)

—^ехр(-2 л/п к =1

= е/ (2) -

Н 2к-1( 2)

(41)

'(п - к )!(2к)!

Принимая во внимание соотношение (41), мы приходим к выражению:

Лс (р, д) =

1

2\/Л

Яе

X

(-1)п

0 п !(п +1/2) 1 4/р

(42)

ГегГ(2)-Ае-22 X Я2к-'(2) ^

>/Л к=1(п - к )!(2к)!

Соотношение (42) допускает существенное упрощение. Для этого представим его в следующем виде:

Ас (р, д) =■

1

гЯе

2л/Л

--Яе [е 2( р, д)е-р ]

[ЕД р, д ^(^/р)]-

(43)

где через Е1 обозначена следующая сумма (см. [5]):

Е1( p, д) = X

(-1)п

0 п !(п +1/2) 14/р

= У

1 £.

2' 4/р

(44)

Вторая сумма, обозначенная через Е2, входящая в формулу (43):

\ п+1/2

е2( p, д) = X

(-1)п

0 п!(п +1/2) 14/р ху Нк-1 ^л//р)

^п - к )!(2к)!,

(45)

вычисляется несколько сложнее.

Хорошо известно, что гамма-функция целого отрицательного аргумента обращается в бесконечность. Учитывая это обстоятельство, преобразуем (45):

е2( p, д) =

(-1)п

X

н2к-1^) (46)

1 п!(п +1/2) 1 4/р 1 к"1(п - к )!(2к)!

п=

п=

В результате указанного преобразования мы можем изменить порядок суммирования по индексам к и п :

2к-Дл/^),

Е,( р, а) = V-

к (2к)!

(-1)п

хУ-

П=1 п !(п - к)!(п +1/2) 1 Мр Из (47) вытекает, что

^2(р, а) =

= ^ (-1)*Н2*) (\

(47)

к=1 (2к)! 1 4/р )

(48)

(-1)"

о т!(т + к)!(т + к +1/2) 14/р

Внутренняя сумма в формуле (48) выражается через интеграл от функции Бесселя. В результате получается выражение для Е2:

а

^ р, а) = £ (-1)ки2 к) 2к (2к)! |

гк/к (г)Сг . (49)

В заключение отметим, что сумма Е2 может быть представлена в виде:

^2 (р,а) =

= %/Лу (-1)к г(к+1/2)я 2к-1^Л//р) х (50)

2 (2к)! Х ' ^ ^

х {[./к (^)Нк-, (г) - /-1 (^)Нк (г)] },

Фр

где Нк (2) - функция Струве.

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

5. Случай больших значений аргументов

а » р »1

В случае очень больших значений обоих аргументов а » р » 1 оценка исследуемых интегралов становится тривиальной. При этом наибольший интерес представляет случай экстремально больших (р > 50) значений аргумента р, поскольку описанный выше алгоритм для случая умеренных р < 5 и произвольных а > 5 на практике пригоден вплоть до значений аргумента р ~50.

В рассматриваемом случае а » р »1 зависимость величины интеграла от верхнего предела интегрирования становится несущественной. Другими словами, мы можем распространить пределы интегрирования в (1) и (2) до бесконечности. Тогда, со-

гласно [5], искомые интегралы вычисляются точно. Согласно [5] интеграл Лл. (р, д) равен

а ( р, а)'

1+с (—1] 2 + I + £ (£ Л]

[ 2 ( 4 р | [ 2 ( 4 р J_

(51)

В рассматриваемом приближении интеграл Лс(р, а) также сводится к табличному:

4 (р, а) * -и®

ей"

( п А

/1 а

24р ,

(52)

Аналитические свойства интеграла ошибок позволяют свести выражение (52) к линейной комбинации интегралов Френеля:

( п Л

/ 7 а

вг/

V (

2у[р

= с

Л ( + £

72этр ) 2пр

+/

с

- £

(53)

Откуда немедленно получается окончательное выражение для интеграла Лс (р, а) в приближении

а » р » 1 :

Лс (р, а) = - •

с

Л ( + £

л/2Пр ) 2пр

(54)

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

6. Случай больших значений аргументов

р » а »1

Если для аргументов р ид выполняется соотношение р ^ д ^ 1, то интегралы л! (р, а) и Ас (р, а) можно вычислить методом перевала. Покажем, что это действительно так. Для этого заметим, что

а ( р, а) =

= -21т } си [е'(р,г+"')-ы'+е'(р,г)-1п' ]. (55)

2 0

В соответствии с основной идеей метода перевала, разложим показатели экспонент в формуле (55) в ряд Тейлора в окрестности г0:

/(рг2 ±аг)- 1пг * /(р^ ±аг0)-

1

- 1п г0 +— 0 2

2'р + -2

(г - г,)2.

(56)

Стационарная точка г0 удовлетворяет уравнению

/(2р*0 ± д) --1 = О

»с

имеющему очевидное решение

г0(±д)=)-которое мы будем записывать как

г>(± д) = +"Т" ^л/р ехр(/ф / 2),

4 р

где

(57)

(58)

(59)

(60)

I 8 р

Ф = -агСап! ——

У д ,

В рамках поставленной задачи выполняются условия р >0 и д >0. В этом случае стационарная

точка г0 расположена далеко от исходного контура интегрирования г е [0,1], проходящего по вещественной оси. Поэтому первое слагаемое в формуле (55) с «перевальной» степенью точности не дает вклада в величину Л, (р, д), а интеграл от второго слагаемого вычисляется элементарно.

Г Л е(

0

Т2Л

(рг - дг )-1п г)

(р'2)

^0

2/р + -2

0 1

(61)

По аналогичной причине из двух корней уравнения (58) в «перевальном» выражении для искомой величины Л! (р, д) следует удержать лишь тот, для которого выполняется условие Яе(г0) > 0 . Тогда в рамках метода перевала мы приходим к следующему выражению

1 1

Л, (р, д) ~ 21тГ С ехр (/(р/2 - дг) - 1пг) и

-[-(2/р+¿02(-д))]-

<I р'0 (-Я)-Я»0(-Я) I

(62)

г0(-д)

и в результате очевидных преобразований мы получаем

Л (р,д) =

= /л I ехр[/ (2(-д) - д^-д) )] = ±12

(63)

л/2/р?02(-д) +1

При извлечении корня из (-1) при переходе от формулы (62) к формуле (63) лист римановой поверхности (знак + или -) выбирается из условия,

что в рассматриваемом случае р » д »1 стационарная точка г0 и 0 , и, следовательно, Л, (р, д) > 0.

Учитывая, что согласно (58) дг0 = 2р»0 + /, приведем выражение (63) к виду

[ехр (-/рг02(-д))

п (64)

Л (р, д) = ±<Ъ • е • Яе \ ,-2-

V2 2/р?2 (-д) +1 ]

откуда получается окончательная формула для вычисления интеграла Л, (р, д) методом перевала

Л, (р, д) = + е • Яе {

,-г"о(-д)-'ф1/2 I

(65)

где

ргК-д) = р |д+4Р ехр(/ф /2) 11 =

= А+д л/р ехр(/ф /2) + р ехр(/ф), 16 р 2

(66)

а модуль и фаза подкоренного выражения в знаменателе формулы (64) определяются соотношением

Р1 ехр(/ф!) = 2/р'02(-Я) +1. (67)

Перейдем к рассмотрению интеграла Лс (р, д) в приближении перевала:

11

Лс (р, д) = --1т | с [

е> ( р» + я' ) - е' (р»2 - яг)

] г-1. (68)

Легко видеть, что в этом приближении выполняется соотношение

Лс(р,я) и-Л,(р,д). (69)

Другими словами, в рассматриваемом случае интеграл Л, (р, д) набирается на первой осцилляции, в то время как интеграл Лс (р, д) набирается на второй осцилляции.

7. Случай больших значений аргументов

р и я » 1

В случае, когда параметры р и д велики и соизмеримы, метод перевала может внести в расчеты заметную погрешность, поскольку при р и д стационарная точка расположена достаточно к нижнему пределу интегрирования. Для того чтобы построить искомый алгоритм, представим интеграл Л, (р, д) в виде:

1 р ¿х

Л, (р, д) = — Г—эт( х) со5(дл/х / р).

2 о х

(70)

Разложим в ряд Тейлора косинус, стоящий под знаком интеграла

>., ч 1 рлх • , ч^ (-1)п Г я2х |

Л (л я) = 2®1П(x)X72пг! — I .

^ х и=0 (2п)! 1 р I

20 х

(71)

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

4 (р, а) = ^к р) -

-1V (-1)п

2 ;П=0(2п + 2)! 1 р

— I | хп Б1п хСх.

(72)

Интеграл, фигурирующий в формуле (72), вычисляется элементарно. В итоге мы приходим к следующему выражению

1

4 (р, а) = ^ р) +

+1V (-1)

,-! V к!

2 п=0 (2п + 2)! 1 р ) 1 к

<[рп-к СОБ

(73)

(р + к п /2 )-5кпп!соБ (пп /2)].

В результате достаточно громоздких, но несложных преобразований мы приходим к следующему выражению для Л! (р, а):

4 (р, а) = р) +

+Т 1 £ 1П"с,п,2|( р) -

2 П=0(2п + 2)!I р ) 1 1

б1п р V (-1)" п! (а

1 (2п + 2)! 1 р ) [(п-1)/21

(74)

|( р) -

-1 V (-1)п (2п)! ( о2 2 ¿0 (4п + 2)! [ р

где [а| - целая часть числа а, а sn (2) (сп (2)) — частичная сумма ряда Тейлора для синуса (косинуса). Полученное выражение весьма похоже на соотношение (26). Однако, в отличие от формулы (26), оно содержит только функции вещественного аргумента.

Перейдем к рассмотрению интеграла Лс (р, а). Проделаем выкладки, полностью аналогичные, тем, что были выполнены для интеграла Л3 (р, д). Для этого преобразуем выражение (2) к виду

Лс (р, а) =1 [—соб( х)$1п(^х/Х7р).

с

20 х

(75)

Разложим подынтегральное выражение в ряд Тейлора. В результате получаем:

лс ( р,а) =

(-1)п

-1 I" — СОБ( х)У 2 о х п=0 (2п +1)!

(

х а

р

2

(76)

Поскольку под интегралом стоят аналитические функции, операции суммирования и интегрирования коммутируют. Следовательно

Лс (р, а) =

=1 (а! Лп+1/2 Ь->/2со5( х)Сх. (77)

2 п=0(2п +1)! 1 р ) 0

Интеграл в формуле (77) удобно выразить через неполную гамма-функцию:

| хп-1/2СОБ( х)Сх = 0

= Ие/п+1/2у(п +1/2, -/р). Тогда

/ о \ п+1/2

=1 V (-1)п (I

(78)

Лс(p, а) = - V (2 +1)! I х

2 п=0 (2п +1)!1 р )

(79)

х Ие ['п+1/2у (п +1/2, -/р)].

Воспользуемся тем, что у(п +1/2, -/р) = = Г(п + 1/2) -Г(п + 1/2, -/р),

а также тем, что для неполной гамма-функции имеет место асимптотическое разложение:

Г(п +1/2, -/р) ж

(80)

ж (-/р)п-1/2 ехр(/р)

п-1/2 1 + /-+...

(81)

Следовательно, для случая р »1 получаем

\ п+1/2

., ч 1 ^ (-1)п | а-Лс (л а) 1

2 п=0 (2п +1)! 1 р

• Ие

^ /—(и+1/2) 1/2 ■ г т

в 2 Г(п +1/2) -/р е,р [1 +...]

(82)

откуда получается выражение, удобное для проведения численных расчетов

Лс ( р, а) * +

2р

+ 1 ^(-1) п г(п+1/2) (а! г+1/2 х (83)

2^2 п=0 (2п +1)! I р )

пп ) . ( пп СОБ| - I-Б1ПI -

2 ) I 2

Заключение

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

1. В работе построены простые аналитические выражения для расчета интегралов специального вида.

2. Проведен численный анализ полученных аналитических выражений.

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

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

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

Приложение

Определения и обозначения, используемые в работе

В настоящей работе используются следующие стандартные определения, обозначения и соотношения, подробно описанные в справочниках [4, 5]:

Н (х) =

¿п г полином Эрмита

= (-1)п ехр( х2)—ехр(-х2)

1 Cп

Рп (2) =--(22 - 1)п - полином Лежандра

п 2пп! й2п

да

Г(2) = Г ехр(-г)г2-1С - гамма-функция

Г( z + n) ^

(z) =--символ Похгаммера

Г( z)

Y (a, z) = J exp(—t)ta ldt - неполная гамма-функция

0

erf (z) = J exp(-t2 )dt - интеграл вероятностей Vn о

dz

inerfc(z) = -in 'erfc(z) - кратный интеграл веро-

ятностей, причем по определению

_ 2 i_1erfc(z) = —;= exp(_z2), Vn

i0erfc( z) = erfc( z) = 1 _ erf (z) 2 z

erfi(z)=—exp(t2)dt - интеграл вероятностей N/л 0

мнимого аргумента .

F(y) = exp(-y2)J exp(t2)dt - интеграл Досона.

0

C(z) = J cos t2 ] dt - интеграл Френеля. S(z) = J sin 12 j dt - интеграл Френеля.

/ \v + 2k

Uv (z, O = £ (-1)k |-Z] Jv+2t © - функция Лом-

меля двух аргументов.

" (z)k

en (z) = V--частичная сумма ряда Тейлора для

к=0 к!

экспоненты (lim en (x) = exp(x)).

n

(— 1)k

sn (x) = ^-x2k+1 - частичная сумма ряда Тей-

k=о(2k + Г)!

лора для синуса (lim sn (x) = sin x).

n (—1)k

cn (x) = V-x2k - частичная сумма ряда Тейлора

k"0 (2k)!

для косинуса (lim cn (x) = cos x).

Литература

1. Ратис Ю.Л., Леонович Г.И. Дифракция светового потока на чувствительных элементах электронно-оптических и оптикоэлектронных датчиков механических перемещений // Компьютерная оптика. - Самара. 1996. - № 16. - С. 74-77.

2. Ратис Ю.Л., Леонович Г.И., Курушина С.Е., Мельников А.Ю. Нелинейные дифракционные искажения оптической функции отклика в кодирующих сопряжениях оптикоэлектронных датчиков // Компьютерная оптика, - Самара-Москва, 1998. - Т. 18. - С. 61-71.

3. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. - М.: Наука, 1987. - 688 с.

4. Абрамовитц М., Стиган И. Справочник по специальным функциям, Москва, Наука, 1979, 832 с.

5. Прудников А.П., Брычков Ю.А., Маричев О.И., Интегралы и ряды, Элементарные функции. - М.: Наука, 1981. - 800 с.

6. Ratis Yu.L., de Cordoba P.F. A code to calculate (high order) Bessel functions based on the continued fractions method. Computer physics communications, 76 (1993) 381.

7. Ratis Yu.L., Segura J., de Cordoba P.F., A code to evaluate Modified Bessel function based on the continued fractions method, Computer Physics Communications, 1997. V. 105. - Р. 263-272.

8. Ratis Yu.L., Bastardo J.L., Abraham Ibrahim S., de Cordoba P.F., Urchueguia J.F.S., Evaluation of Fresnel integrals based on the continued fractions method, Applied Mathematics Letters, 2005. - № 18. - P. 23-28/

k=0

d