CC BY
41
УДК 535.42
ОПИСАНИЕ СПОСОБА ВЫЧИСЛЕНИЯ ОДНОГО КЛАССА СПИРАЛЬНЫХ ФУНКЦИЙ
© 2013 А.В. Устинов
Институт систем обработки изображений РАН, г. Самара
Поступила в редакцию 30.10 2013
Статья посвящена вопросу вычисления значений неэлементарных функций, являющихся обобщением интегральных синуса и косинуса и функций Френеля. Отмечена возможность задавать комплексное значение параметра и/или аргумента функции. Также мы обсуждаем использование этих функций в практических расчётах и аспекты их вычисления.
Ключевые слова: осциллирующая функция; спираль; степенной ряд; асимптотический ряд.
ВВЕДЕНИЕ
В ряде случаев [1-3] при вычислениях появляются выражения, содержащие под знаком интеграла осциллирующую функцию. Достаточно часто их можно преобразовать к комбинации элементарных и спиральных функций.
В данной работе произвольной спиральной функцией будем называть функцию вида
fc ,s (t,a) = jV.fC? & )) 2 (u, a)du , (i)
u 0 I J
где a - вектор параметров, а функции g 1, g2 таковы, что при t — ^ интеграл сходится, при этом g 1 не имеет конечного предела при t —У ^ . Если эти условия выполнены, то при параметрическом задании x(t) = fc(t) y(t) = fs(t) получим кривую с формой закручивающейся спирали. Стандартными примерами могут являться sici-спираль или спираль Корню. Если условие отсутствия конечного предела не выполнено, то кривая не будет похожа на спираль - с ростом аргумента кривая стремится к предельной точке без закручивания вокруг неё. Если при бесконечном аргументе интеграл расходится, то получим раскручивающуюся спираль. (Очевидно, имеются целые семейства спиральных функций, не описываемых формулой (1), так как косинус и/или синус можно заменить другими осциллирующими функциями, например, функциями Бесселя.) Далее мы ограничимся лишь одним классом из семейства (1).
1. ОБСУЖДАЕМЫЙ КЛАСС СПИРАЛЬНЫХ ФУНКЦИЙ
Мы рассмотрим следующий класс спиральных функций
Устинов Андрей Владимирович, аспирант, ведущий программист лаборатории лазерных измерений. E-mail: andr@smr.ru
cos I sin I
(t, k) = J
cos | du (u
(2)
где u0 = 0 при 0 < k < 1, и u0 = при k > 1. Если k < 0 (при этом спираль является раскручивающейся), то интегрированием по частям мы можем выразить требуемую спиральную функцию в виде комбинации элементарных функций и спиральной функции с параметром 0 < k < 1. Назовём эти функции обобщёнными интегральным синусом и косинусом. При k = 1 имеем стандартные интегральный синус и косинус, а при k = 1/2 с точностью до постоянного множителя - функции Френеля.
Функции (2) встречаются при вычислении электромагнитного поля, формируемого оптическим элементом в виде обобщённой параболической линзы. Кроме того, спираль, соответствующая функциям (2), удобна для тестирования некоторых алгоритмов обработки изображений, в которых вычисляется кривизна кривой. Дело в том, что для такой спирали получаются простые формулы для кривизны (в текущей точке) и, что бывает намного реже, для длины дуги от некоторой начальной точки до текущей. Более того, аналитически получаются обратные формулы - при каком значении параметра будет заданная кривизна и длина дуги. Примеры спиралей приведены на рис. 1. Отметим, что при k < 1 полная длина спирали бесконечна, а при k > 1 конечна.
В ряде работ похожие функции уже рассматривались. Так, в [4] упомянуто однопараметри-ческое семейство функций, которые при частном значении параметра выражаются через cos I(t,3 / 2) и sin I(t,3 /2). В [5] рассмотрен намного более близкий случай. Там введены обобщён-
t Г—s|
С 1/иа р
ные интегралы Френеля N зр ^ с р > 1,
получены приближённые формулы для большого значения аргумента, на основе которых вычис-
u
а)
б)
в) г)
Рис. 1. Примеры спиралей для различных значений параметра. Начало координат взято в предельной точке - центре спирали.
а) к=0,2; Ь в диапазоне 0...25; масштаб 100; центр: х=0,360 у=1,107
б) к=0,8; Ь в диапазоне 1.25; масштаб 200; центр: х=4,366 у=1,419
в) к=1,3; Ь в диапазоне 1.25; масштаб 300; центр: х=0 у=0
г) к=2; Ь в диапазоне 2.20; масштаб 1000; центр: х=0 у=0
лены некоторые фрактальные характеристики соответствующей спирали, называемой p-клото-идой. Эти интегралы с точностью до постоянного множителя сводятся к cos I (t,1 -1/p) и sin I(t,1 -1/ p), таким образом, покрывая наш диапазон 0 < k < 1. Однако в [5] не упомянуты вычислительные аспекты получения значений этих интегралов. Кроме того, наш класс шире - имеется область параметров k > 1. В принципе в этой области, если параметр не равен 1, можно привести функцию к комбинации элементарных и функции с параметром, меньшим 1. Но такое преобразование не всегда оправдано с точки зрения вычислительной сложности (см. ниже) и не очень удобно для программной реализации. Наконец, главное расширение состоит в возможности да-
вать для г и к комплексные значения, естественно, с определёнными ограничениями, о которых далее будет сказано. Комплексное значение г получается, например, при приближённо-аналитическом вычислении амплитуды электромагнитного поля модифицированным методом стационарной фазы [6] при освещении оптического элемента гауссовым пучком.
2. ПРИБЛИЖЁННО-АНАЛИТИЧЕСКИЕ ВЫРАЖЕНИЯ ДЛЯ ДАННОГО КЛАССА ФУНКЦИЙ
Выражение будем записывать отдельно для обобщённого интегрального синуса/косинуса в обоих диапазонах: о < к < 1 и к > 1.
2.1. Обобщённый интегральный косинус при 0 < к < 1
Представление в виде степенного ряда:
5I (t, k) = J
cos u 1-k -du = t
(-1) nt2 n
,(2n)!<2n +1 - k)
. (3)
Этот ряд имеет бесконечный радиус сходимости, но при большом значении аргумента сходится медленно, поэтому в этом случае более целесообразно использовать другое выражение.
Асимптотическое разложение при большом аргументе t.
Представим искомую величину в виде разно-
7 cos u 7 cosu сти cos I(t, к) = J —^"u - J —^"u . Первое сла-
0 u t u гаемое вычисляется точно [7], а втором производится интегрирование по частям (тригонометрическая функция является интегрируемым множителем) некоторое число раз. В результате после преобразований получим представление:
cos I (t, k) = r(1-k)sm(fa/2)^ 1-4-; .1 + 4+1 -J-All cost( 1 k 1 1 k 1 1 k 1 1
——I k ■—Ak+1 ■—+a;+-1 ■— - 4-1 ■—+
tk I t t3 k+4 t5 k+6 t7
(4)
Здесь А^ - число размещений из и по т. Вопрос о выборе между формулами (3) и (4) и требуемом количестве слагаемых обсуждается ниже в отдельном параграфе.
2.2. Обобщённый интегральный синус при 0 < к < 1
Здесь всё аналогично пункту 2.1. Представление в виде степенного ряда:
11 (t, k) = J
sin u , 2-k -du =t2 k
(-1) nt2n
n~0 (2n + 1)!-(2n + 2- k)
.(5)
Асимптотическое разложение при большом аргументе t .
Представив искомую величину в виде разно. г, , ч 7 sin u , 7 sin u , сти sin I (t, k) = J——du- J——du , аналогично
0 u t u cos I(t, k) (с использованием формулы из [7]), получим представление:
sinI (t, k) =Г(1 -k)cos(fcrt /2)-
1-Akk+1 ■ -+Ak+3 ■ --Akk++5 ■ -+...I-
, 1 .k-1 1 .k-1 1 .k-1 1 , k-Г Ak+2 —Г + Ak+4 ^^-Ak+6 ^^ + ■ t \ t t3 t t'
sint
.(6)
не нулю, то проще получается разложение при большом аргументе г. А именно, сохраняют силу формулы (4) и (6) без свободного слагаемого, так как при выводе разложения не использовалось условие о < к <1.
Степенной ряд при не очень большом аргументе г получить несколько сложнее. Вначале рассмотрим нецелые значения параметра к .
Пусть т = [к] - целая часть числа к. Тогда т-кратным интегрированием по частям (тригонометрическая функция является дифференцируемым множителем) мы сведём задачу к случаю, когда степень знаменателя меньше единицы, который уже рассмотрен.
Рассмотрим обобщённый интегральный косинус.
cosI (t, k) =
cost
sint
cost
(1- k)tk-1 (1- k)(2 - k)tk-2 (1- k)(2 - k)(3 - k)tk-3
sin t
(1 -k)(2 - k)(3 - k)(4 -k)tk
+ (сумма2). (7 а)
Далее в первой сумме чередование знаков у cos t и sint периодически продолжается, число слагаемых в первой сумме равно m. Подобно структуре первой сумме вторая сумма имеет различный вид в зависимости от величины остатка от деления m на 4.
(sin I(+7, k - m) - sin I (t, k - m)); m(mod4) = 1 (cos I (+7, k - m) - cos I (t, k - m)); m(mod4) = 2 (sin I (+7, k - m) - sin I (t, k - m)); m(mod4) = 3
(сумма2) =
(1-k) ■■■ (m - k) 1
+-*
(1 - k) ■■■ (m - k) 1
+-*
(1 - k) ■■■ (m - k)
1
(1-k) ■■■ (m - k)
:(cos I(+7, k - m) - cos I(t, k - m)); m(mod4) = 0
(7б)
Представление (7а), (7б) верно для любого значения аргумента, поэтому, в принципе, мы всегда можем перейти к параметру меньшему единицы, но при большом значении аргумента это делать нецелесообразно. Степенной ряд для функций с параметром к - т даётся формулой (3) или (5), а значения с параметром к -т и бесконечным аргументом - это свободные слагаемые в (4) или (6). В силу громоздкости явный вид степенного ряда, получаемого по (7а), (7б), выписывать не будем.
Аналогично для обобщённого интегрального синуса.
2.3. Обобщённые интегральные синус и косинус при к > 1
В данном диапазоне параметра удобнее рассматривать обе функции в одном пункте. Так как здесь постоянный предел равен бесконечности, а
sin I (t, k) =
sint
cost
sint
(1 - k)tk-1 (1- k)(2-k)tk-2 (1- k)(2-k)(3-k)tk-3
cos t
(1 - k )(2 - k )(3 - k )(4 - k )tk
- +...
+ (сумма2). (8 а)
Число слагаемых в первой сумме равно m.
u
0
+
*
u
k
+
(сумма 2)
(1 - к) ••• (m - к) 1
(1 - к) ••• (m - к)
+--* (cos I, k - m) - cos I(t, k - m)); m (mod 4) = 1
(1 -k) •• • (m -k)
+-1-* (sin I, k - m) - sin I(t, k - m)); m(mod 4) = 2
(1 -k) •• • (m -k)
* (cos I, k - m) - cos I(t, k - m)); m(mod 4) = 3
* (sin I, k - m) - sin I(t, k - m)); m (mod 4) = 0
(8б)
Пусть теперь параметр k является целым. К сожалению, мы не можем придти к нулевому значению параметра, но можем его сделать равным единице - все формулы (7), (8) сохраняют силу и для целого k > 1, если положить m = k -1. Осталось рассмотреть только случай k = 1.
Для sin I(t,1) ряд легко получается переходом к нижнему пределу, равному нулю, как при параметре, меньшем единицы. (На самом деле, для обобщённого интегрального синуса можно было бы сохранить нижний предел, равный нулю, для значений параметра вплоть до двух.)
„ rsmu rsinu , rsinu , rsinu , ,„/<-> \ ni (/,1) = J-du = J-du-J-du = J-du-n /2 .(9 а)
• sinu u
sinu u
Разложение интеграла в ряд производится по формуле (5), если подставить к = 1.
Получить таким же способом ряд для cos I(t,1) невозможно, так как оба интеграла в разности, аналогичной имеющейся в (9)расходятся. Собственно, поэтому мы и поменяли нижний предел на бесконечный. Однако cos I(t,1) совпадает со стандартным интегральным косинусом, для которого в [8] есть готовая формула:
I (t ,1) = J du = ln t + C + У J u
(-1) nt2 n
"1 (2n)!-2n
Здесь C = 0,5772... - постоянная Эйлера. Ряд для sin I (/,1) тоже можно было взять в готовом виде. Заметим, что в (9б) слагаемые степеней второй и выше такие же, как в формуле (3) при подстановке k = 1. К сожалению, в [8] нет хотя бы краткого (нестрогого) доказательства формулы (9б). Правда, логарифмическое слагаемое и члены с положительными степенями получаются тривиально, вопрос лишь в доказательстве того, что слагаемое нулевой степени равно постоянной Эйлера. Отсутствие доказательства помешает нам ниже, когда готовой формулы нет, но доказательство можно было бы сделать по аналогии с (9б).
3. НЕКОТОРЫЕ ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ
В этом параграфе мы опишем некоторые подробности численной реализации аналитических представлений из предыдущего параграфа. Далее в параграфе символом е обозначена допустимая абсолютная погрешность. Аргумент t везде предполагается положительным.
3.1 . Формула (3)
Перепишем её правую часть в виде tк • ^ ап.
п = 0
Ряд является знакочередующимся, поэтому для достижения желаемой точности достаточно, что-
бы было \an \ < —— (=е1). Для ускорения вычис-t
лений члены ряда вычисляются рекуррентно:
1 an
2n +1 - к
1 - к an (2n + 1)(2n + 2) 2n + 3 - к
. (10)
3.2. Формула (4)
Перепишем её правую часть в виде
sint
3 -
cost
^2
2Г(к)^(кгс/2) ^ 1 ^
Так как <1, <1, поэтому достаточно найти суммы £2 с точностью до е1 = £• ^ . Так как ряды, определяющие эти суммы, являются не сходящимися, а асимптотическими, то обрывать их следует не позже, чем когда слагаемые перестанут убывать. В обозначениях, таких же, как в (10) имеем.
Для :
a 0 =1
a1
-n+L =---(к + 2п)(к + 2n +1)
an t2
Слагаемое an+1 не берём, если
■ 1. Ми-
(9б) нимальное значение слагаемого an примерно
■ i i (t0 -2)! 1Г равно min an ~-, где 10 = J/[ - целое сверху,
например ]2,1[= 3 ]1[= 2 .
Для S2:
к an+
t an
1
(к + 2n + 1)(к + 2n + 2)
Слагаемое ап+1 не берём, если
нимальное значение слагаемого а для суммы ¿1.
> 1. Ми-
Итак, если
(t0 - 2)!
(t 0)t0 - 2
< e1, то применяется фор-
мула (4), иначе формула (3). Естественно, проверка имеет смысл только при 10 > 3 , при меньшем значении всегда применяется (3).
3.3. Формула (5)
Здесь аналогично пункту 3.1 перепишем её
правую часть в виде t2~к • ^ ап . Для достижения
2
t
a=
0
n
0
0
a
n+1
a
n
cos
(t 0)0
a=
0
2
t
a
n+1
a
n
n=0
желаемой точности достаточно, чтобы было I £
|ап| < 2_к (= е1). Для ускорения вычислений члены ряда вычисляются рекуррентно:
1 an+, t2 2n + 2 - k
2 - k an (2n + 2)(2n + 3) 2n + 4 - k
.(11)
3.4. Формула (6)
Аналогично пункту 3.2 перепишем её правую
cos t sin t часть в виде Г(1 -k)cos(kn/2)--S1--S2 .
tk tk
Здесь суммы S1; S2 такие же, как в пункте 3.2. Поэтому верны те же результаты.
Формулы (7), (8), (9а) комментариев не требуют.
3.5. Формула (9б)
Здесь всё очень похоже на пункт 3.1. Небольшое изменение есть только в рекуррентном вычислении членов ряда:
t_ 4
t
(2n + 1)(2n + 2) n +1'
(12)
кте 3.3 е1 = -
i 2-k
. Больше изменений будет в пун-
получим, что u = u 1 , то есть мнимая часть параметра не влияет на сходимость интеграла. Поэтому все рассуждения предыдущих параграфов сохраняют силу, если в неравенствах, определяющих диапазон значений параметра и в формулах из параграфа 3 для е1 заменить k на k1 . Единственная сложность будет только в аналитическом представлении степенного ряда из пункта 2.3, когда k1 - целое число. Если требуется только одна из двух функций, то в половине случаев проблемы не будет, но благоприятные случаи для разных функций не перекрываются, поэтому если нужны обе функции, то проблема появится при любом целом k1 - см. (7б) и (8б). Причина в том, что уже было упомянуто в конце параграфа 2, мы не можем вычислить cos I(t,1 + ik2), так как нечем заменить формулу (9б). Большую часть слагаемых заменяющей формулы мы выписать можем:
cosI (t,1 + ito) = J
cosu - ,1+iro
4. ОСОБЕННОСТИ ВЫЧИСЛЕНИИ ДЛЯ СЛУЧАЯ ВЕЩЕСТВЕННОГО к И КОМПЛЕКСНОГО t
В основном верны рассуждения из параграфа 3. В пунктах 3.1 (то же в 3.5) и 3.3 изменение
незначительное - ставится знак модуля комплек-
£
сного числа. В пункте 3.1 будет £1 = , а в пун-
ктах 3.2 (то же в 3.4). Одно незначительное: 10 = ]abs(t)[, другое намного существеннее и оно связано с тем, что при комплексном аргументе нельзя считать, что синус и косинус по модулю не превосходят единицы. Вместо этого можно доказать, что выполняются приближённые неравенства: |sint| < et2/-J2, |cost| < et2 /-J2 , где 12 = |lmt| . Поэтому изменится выражение для
k
е-1 г-
е 1 2 .
e2
5. ИЗМЕНЕНИЯ ДЛЯ СЛУЧАЯ КОМПЛЕКСНОГО k
Исходя из равенства uk = uk1+ik2 = uk1 • eik2lnu,
du = sin(colnt) + i cos(alnt) + A + V-(13)
„=1(2и)!-(2и >
Но в этой формуле неизвестна величина A , не зависящая от t. Пока в программной реализации эта проблема никак не решается. В принципе можно использовать приём, который иногда используется для вычисления функции Бесселя Yv при целом индексе. Пользуясь свойством непрерывности, пишут приближённое равенство Yv = (Zv_5 + Yv+s)/2 , а при нецелом индексе вычисления намного проще. Правда, в отличие от функций Бесселя, есть исключение - такой приём подходит для всех целых k1 , кроме k1 =1 .
СПИСОК ЛИТЕРАТУРЫ
1. Борн М, Вольф Э. Основы оптики. М.: Наука, 1973. 720 с.
2. Хонина С.Н. Волотовский С.Г. Фраксикон - дифракционный оптический элемент с конической фокальной областью // Компьютерная оптика. 2009. Т.33. №4. С. 401-411.
3. Анализ осевого распределения, формируемого фрак-сиконом в параксиальном и непараксиальном случаях / С.Н. Хонина, А.В. Устинов, А.В. Карсаков // Известия Самарского научного центра РАН. 2013. № 15(4). С. 18-25.
4. Дагуров П.Н., Дмитриев А.В. Дифракция Френеля-Кирхгофа на щели в проводящем экране при скользящем падении волн // III Всероссийская конференция "Радиолокация и радиосвязь" - ИРЭ РАН, 26-30 октября 2009 г. С.693-696.
5. Generalized Fresnel integrals and fractal properties of related spirals / L. Korkut, D. Vlah, D. Zubrinic, V. Zupanovic // Applied Mathematics and Computation. 2008. №206. P. 236-244.
6. Устинов А.В., Хонина С.Н. Обобщённая линза: анализ осевого и поперечного распределения // Компьютерная оптика. 2013. Т.37, №3. С. 307-315.
7. Диткин В.А., Прудников А.П. Интегральные преоб-
a=
0
a
n
n+1
a
n
е
разования и операционное исчисление. М.: Гос. Из- 8. Специальные функции / Е. Янке, Ф. Эмде, Ф. Лёш. во физико-математической литературы, 1961. 524 с. М.: Наука, 1977. 344 с.
DESCRIPTION OF THE CALCULATION WAY FOR ONE CLASS OF SPIRAL FUNCTIONS
© 2013 A.V. Ustinov
Images Processing Systems Institute of the RAS, Samara
This paper is devoted to issue of nonelementary functions value calculating. These functions are the generalization of integral sine, integral cosine and Fresnel functions. We note a possibility to preset the complex value of function parameter and/or argument. As well we consider use of these functions in the practical computations and aspects of their calculation. Key words: oscillating function; helix; power series; asymptotic series.
Andrey Ustinov, Leading Programmer of the Laser Measurements Laboratory. E-mail: andr@smr.ru
