Bessel Functions of the First and Second Kinds . Boost , Math Toolkit 2.5.0 , Bessel Functions

Bessel Functions of the First and Second Kinds

Synopsis

<#include<boost/math/special_functions/bessel.hpp>>

template <class T1, class T2>
calculated-result-type cyl_bessel_j(T1 v, T2 x);
template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_j(T1 v, T2 x, const Policy&);
template <class T1, class T2>
calculated-result-type cyl_neumann(T1 v, T2 x);
template <class T1, class T2, class Policy>
calculated-result-type cyl_neumann(T1 v, T2 x, const Policy&);

Description

Функцииcyl_bessel_jиcyl_neumannвозвращают результат Бесселя функции первого и второго рода соответственно:

cyl_bessel_j(v, x) = J_v(x)

cyl_neumann(v, x) = Y_v(x) = N_v(x)

где:

Тип возврата этих функций вычисляется с помощью.правила расчета типа результата, когда Т1 и Т2 являются различными типами. Функции также оптимизированы для относительно распространенного случая, когда T1 является целым числом.

Конечный аргументПолитикаявляется необязательным и может использоваться для контроля поведения функции: как она обрабатывает ошибки, какой уровень точности использовать и т. д. См. документациюдля более подробной информации.

Функции возвращают результатdomain_errorвсякий раз, когда результат не определен или сложен. Дляcyl_bessel_jэто происходит, когда<x< 0>и v не является целым числом, или когда<x== 0>и<v !=0>. Дляcyl_neumannэто происходит, когда<x<= 0>.

Следующий график иллюстрирует циклическую природу J_v:

Следующий график показывает поведение Y_v: Это также циклично для большихx, но имеет тенденцию к -∞ для маленькихx:

Testing

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

Accuracy

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

Table 6.40. Error rates for cyl_bessel_j (integer orders)

	Microsoft Visual C++ версия 12.0 Win32 двойная	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux Double	Солнечный компилятор версии 0x5130 Солнечный солярис
Bessel J0: Mathworld Data (целая версия)	Макс = 2,52ε (Средний = 1,2ε) :Макс = 1,89ε (Средний = 0,988ε)]	Max = 6.55ε (Mean = 2.89ε) :Max = 5.04ε (Mean = 1,78ε)И другие сбои.	Max = 0ε (Mean = 0ε) GSL 1.16:Max = 1.12ε (Mean = 0.488ε)] Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 1,12ε (Средний = 0,568ε)	Макс = 6.55ε (Средний = 2.86ε)
Bessel J0: Mathworld Data (Tricky cases) (целая версия)	Max = 1e+007ε (Mean = 4.09e+006ε) :Max = 2.54e+008ε (Mean = 1.04e+008ε)]	Макс = 1,63e+08ε (Средний = 6,67e+07ε) :Макс = 4,79e+08ε (Средний = 1,96e+08ε)]	Макс = 7.98e+04ε (Средний = 3.26e+04ε) GSL 1.16:Макс = 1e+07ε (Средний = 4.11e+06ε)] [154e+07ε] Макс = 2,54e+08ε (Средний = 1.04e+08ε]]	Макс = 1,64e+08ε (Средний = 6,69e+07ε)
Бессел J1: Mathworld Data (целая версия)	Макс = 1,73ε (Средний = 0,976ε) :Макс = 11,4ε (Средний = 4,15ε)]	Max = 2.66ε (Mean = 1.38ε) :Max = 6.1ε (Mean = 2.95ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 1,89ε (Mean = 0,721ε)] Rmath 3.0.2:Макс = +INFε (Mean = +INFε)И другие неудачи.) Цефе:Макс = 2,88ε (Средний = 1,12ε)	Макс = 1,44ε (Средний = 0,637ε)
Bessel J1: Mathworld Data (неудобные случаи) (целая версия)	Макс = 3,23e+004ε (Средний = 1,45e+004ε) :Макс = 1,44e+007ε (Средний = 6,5e+006ε)]	Макс = 2.18e+05ε (Средний = 9.76e+04ε) :Макс = 2.15e+06ε (Средний = 1,58e+06ε])	Макс = 106ε (Mean = 47.5ε) GSL 1.16:Макс = 1.26e+06ε (Mean = 6.28e+05ε) Макс = 2.93e+06ε (Mean = 1.7e+06ε]] Цефес:Макс = 9.56e+05ε (Mean = 4.99e+05ε]]	Макс = 2,18e+05ε (Средний = 9,76e+04ε)
Бессель JN: Mathworld Data (целая версия)	Макс = 14,7ε (Mean = 5,4ε) :Макс = +INFε (Mean = +INFε)И другие неудачи.	Max = 6.85ε (Mean = 3.41ε) :Max = 2.13e+19ε (Mean = 5.16e+18ε)И другие сбои.	Max = 0ε (Mean = 0ε) GSL 1.16:Max = 6.9e+05ε (Mean = 2.53e+05ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = +INFε (Средний = +INFε)И другие неудачи.	Макс = 463ε (Средний = 112ε)

Table 6.41. Error rates for cyl_bessel_j

	Microsoft Visual C++ версия 12.0 Win32 двойная	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux Double	Солнечный компилятор версии 0x5130 Солнечный солярис
Бессел J0: Mathworld Data	Макс = 2,52ε (Средний = 1,2ε)	Макс = 6.55ε (Средний = 2.89ε) :Макс = 5.04ε (Средний = 1,78ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 0.629ε (Mean = 0.223ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 1.12ε (Средний = 0,568ε)	Макс = 6.55ε (Средний = 2.86ε)
Бессел J0: Mathworld Data	Макс = 1e+007ε (Средний = 4.09e+006ε)	Макс = 1,63e+08ε (Средний = 6,67e+07ε) :Макс = 4,79e+08ε (Средний = 1,96e+08ε)]	Макс = 7.98e+04ε (Mean = 3.26e+04ε) GSL 1.16:Макс = 6.5e+07ε (Mean = 2.66e+07ε)] [404e+07ε] Макс = 1.04e+08ε) Цефес:Макс = 2,54e+08ε (Mean = 1.04e+08ε]	Макс = 1,64e+08ε (Средний = 6,69e+07ε)
Бессел J1: Mathworld Data	Макс = 1,73ε (Средний = 0,976ε)	Max = 2.66ε (Mean = 1.38ε) :Max = 6.1ε (Mean = 2.95ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 6.62ε (Mean = 2.35ε)И другие сбои. Rmath 3.0.2:Макс = +INFε (Средний = +INFε)И другие неудачи.) Цефес:Макс = 2,88ε (Средний = 1,12ε)	Макс = 1,44ε (Средний = 0,637ε)
Bessel J1: Mathworld Data (неудобные случаи) (целая версия)	Макс = 3,23e+004ε (Средний = 1,45e+004ε)	Макс = 2.18e+05ε (Средний = 9.76e+04ε) :Макс = 2.15e+06ε (Средний = 1,58e+06ε])	Макс = 106ε (Mean = 47.5ε) GSL 1.16:Макс = 8.75e+05ε (Mean = 5.32e+05ε)] [06ε] [Cephes:Макс = 9.56e+05ε (Mean = 4.99e+05ε]]	Макс = 2,18e+05ε (Средний = 9,76e+04ε)
Бессел JN: Mathworld Data	Макс = 14,7ε (Средний = 5,4ε)	Max = 6.85ε (Mean = 3.41ε) :Max = 2.13e+19ε (Mean = 5.16e+18ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 6.9e+05ε (Mean = 2.15e+05ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 5,53e+05ε (Средний = 1,9e+05ε)	Макс = 463ε (Средний = 112ε)
Бессел Дж: Данные математического мира	Макс = 14,9ε (Средний = 3,82ε)	Макс = 14,7ε (Mean = 4.05ε) :Макс = 3,49e+05ε (Mean = 7,89e+04ε)И другие сбои.	Max = 10ε (Mean = 2.19ε) GSL 1.16:Max = 2.39e+05ε (Mean = 5.24e+04ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 5,47e+05ε (Средний = 1,3e+05ε)	Макс = 14,7ε (Средний = 4,12ε)
Бессель Дж: Данные математического мира (большие значения)	Макс = 9.31ε (Средний = 5.52ε)	Max = 607ε (Mean = 305ε) :Max = 34.9ε (Mean = 17.4ε)И другие сбои.	Max = 0,536ε (Mean = 0,268ε) GSL 1.16:Max = 4,91e+03ε (Mean = 2,46e+03ε)И другие сбои.) [Rmath 3.0.2:Max = 35.9ε (Mean = 18.1ε)] Цефес:Max = +INFε (Mean = +INFε)И другие сбои.	Макс = 607ε (Средний = 305ε)
Бессел JN: Случайные данные	Макс = 17.5ε (Средний = 1,46ε)	Макс = 50.8ε (Средний = 4.15ε) :Макс = 1.12e+03ε (Средний = 88.7ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 75.7ε (Mean = 5.36ε)] Макс = 3.93ε (Mean = 1.22ε) Цефес:Макс = 91.4ε (Mean = 6.47ε)	Макс = 99,6ε (Средний = 22ε)
Бессель Дж: Случайные данные	Макс = 9.24ε (Средний = 1.36ε)	Max = 9.81ε (Mean = 1.59ε) :Max = 501ε (Mean = 52.3ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 15.5ε (Mean = 3.33ε)И другие сбои.) Rmath 3.0.2:Max = 6.74ε (Mean = 1.3ε)] Цефес:Max = 16.7ε (Mean = 2.5ε)]	Макс = 260ε (Средний = 34ε)
Бессель Дж. Данные (Tricky big values)	Макс = 59,2ε (Средний = 8,67ε)	Макс = 785ε (Mean = 94.2ε) :Макс = 5.01e+17ε (Mean = 6.23e+16ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 2.48e+05ε (Mean = 5.11e+04ε)] Rmath 3.0.2:Макс = 71.6ε (Mean = 11.7ε)] Цефес:Макс = 2.48e+05ε (Mean = 3.02e+04ε]]	Макс = 785ε (Средний = 97,4ε)

Table 6.42. Error rates for cyl_neumann (integer orders)

Microsoft Visual C++ версия 12.0
Win32
двойная

GNU C++ версия 5.1.0
Linux
Double

GNU C++ версия 5.1.0
Linux
длинный двойной

Солнечный компилятор версии 0x5130
Солнечный солярис

Y0: Mathworld Data (целая версия)

Max = 4.61ε (Mean = 2.29ε)

:Max = 5.37e+003ε (Mean = 1,81e+003ε)]

Макс = 0ε (Mean = 0ε)

GSL 1.16:Макс = 6.46ε (Mean = 2.38ε)]
Rmath 3.0.2:Макс = 167ε (Mean = 56.5ε)]
Цефес:Макс = 5.37e+03ε (Mean = 1.81e+03ε)

Макс = 5.59ε (Mean = 2.54ε)

:Макс = 2.05e+05ε (Mean = 6.87e+04ε)]

Макс = 5.53ε (Средний = 2.4ε)

Y1: Mathworld Data (целая версия)

Макс = 4,75ε (Средний = 1,72ε)

:Макс = 1,86e+004ε (Средний = 6,2e+003ε)]

Макс = 0ε (Mean = 0ε)

GSL 1.16:Макс = 1.51ε (Mean = 0.839ε)]
Макс = 193ε (Mean = 64.4ε)]
Цефес:Макс = 1.86e+04ε (Mean = 6.2e+03ε)

Макс = 12,7ε (Mean = 4,34ε)

:Макс = 9,71e+03ε (Mean = 4.08e+03ε)

Макс = 6.33ε (Средний = 2,29ε)

Yn: Mathworld Data (целая версия)

Max = 35ε (Mean = 11.8ε)

:Max = 2.49e+005ε (Mean = 8.14e+004ε)]

Макс = 0.993ε (Средний = 0.314ε)
GSL 1.16:Макс = 2.41e+05ε (Средний = 7.62e+04ε)]
Цефес:Макс = 2.49e+05ε (Средний = 8.14e+04ε)

Max = 55.2ε (Mean = 17.7ε)

:Max = 2.2e+20ε (Mean = 6.97e+19ε)И другие сбои.

Макс = 55.2ε (Средний = 17.8ε)

Table 6.43. Error rates for cyl_neumann

	Microsoft Visual C++ версия 12.0 Win32 двойная	GNU C++ версия 5.1.0 Linux Double	GNU C++ версия 5.1.0 Linux длинный двойной	Солнечный компилятор версии 0x5130 Солнечный солярис
Y0: Mathworld Data	Макс = 4.61ε (Средний = 2,29ε)	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 60.9ε (Mean = 20.4ε)] Rmath 3.0.2:Макс = 167ε (Mean = 56.5ε)]	Макс = 5.59ε (Mean = 2.54ε) :Макс = 2.05e+05ε (Mean = 6.87e+04ε)]	Макс = 5.53ε (Средний = 2.4ε)
Y1: Mathworld Data	Макс = 4,75ε (Средний = 1,72ε)	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 23.4ε (Mean = 8.1ε)] Rmath 3.0.2:Макс = 193ε (Mean = 64.4ε)	Макс = 12,7ε (Mean = 4,34ε) :Макс = 9,71e+03ε (Mean = 4.08e+03ε)	Макс = 6.33ε (Средний = 2,29ε)
Yn: Mathworld Data	Макс = 35ε (Средний = 11.8ε)	Max = 0.993ε (Mean = 0.314ε) GSL 1.16:Max = 2.41e+05ε (Mean = 7.62e+04ε)И другие сбои.) Rmath 3.0.2:Max = 1.24e+04ε (Mean = 4e+03ε)	Max = 55.2ε (Mean = 17.7ε) :Max = 2.2e+20ε (Mean = 6.97e+19ε)И другие сбои.	Макс = 55.2ε (Средний = 17.8ε)
Yv: Mathworld Data	Макс = 7.89ε (Средний = 3.27ε)	Max = 10ε (Mean = 3.02ε) GSL 1.16:Max = 1,07e+05ε (Mean = 3,22e+04ε)И другие сбои.) Rmath 3.0.2:Max = 1,05e+03ε (Mean = 326ε)	Max = 10,7ε (Mean = 4,92ε) :Max = 3,49e+15ε (Mean = 1,05e+15ε)И другие сбои.	Макс = 10,7ε (Средний = 5,1ε)
Yv: Mathworld Data (большие значения)	Макс = 0,682ε (Средний = 0,35ε)	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 60.8ε (Mean = 23ε)И другие неудачи. Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.	Макс = 1.57ε (Mean = 1.17ε) :Макс = 43.2ε (Mean = 16.3ε)И другие неудачи.	Макс = 1.57ε (Средний = 1.24ε)
Y0 и Y1: Случайные данные	Макс = 4.17ε (Средний = 1.24ε)	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 34.4ε (Mean = 8.9ε)] Rmath 3.0.2:Макс = 83ε (Mean = 14.2ε)]	Макс = 11.8ε (Средний = 3.28ε) :Макс = 2,59e+03ε (Средний = 500ε)]	Макс = 10,8ε (Средний = 3,04ε)
Yn: Случайные данные	Макс = 117ε (Средний = 10.2ε)	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 500ε (Mean = 47.8ε)] Rmath 3.0.2:Макс = 691ε (Mean = 67.9ε)]	Max = 338ε (Mean = 28.2ε) :Max = 4.01e+03ε (Mean = 348ε)	Макс = 338ε (Средний = 27,5ε)
Yv: Случайные данные	Макс = 1.23e+003ε (Средний = 69.9ε)	Макс = 1,53ε (Средний = 0,102ε) GSL 1.16:Макс = 1,41e+06ε (Средний = 7,67e+04ε)] Rmath 3.0.2:Макс = 1,79e+05ε (Средний = 9,64e+03ε)	Макс = 2.08e+03ε (Mean = 149ε) :Макс = +INFε (Mean = +INFε)И другие неудачи.	Макс = 2.08e+03ε (Средний = 149ε)

Обратите внимание, что для большиххэти функции во многом зависят от точности функций<std::sin>и<std::cos>.

Сравнение с GSL иCephesинтересно: какCephes, так и эта библиотека оптимизируют случай целого порядка — приводя к идентичным результатам — просто используя общий случай, по большей части немного более точный, хотя, как отмечено лучшей точностью GSL в случаях целых аргументов. Эта реализация, как правило, работает намного лучше, когда аргументы становятся большими,Цефес, в частности, дает некоторые удивительно неточные результаты с некоторыми тестовыми данными (без существенных правильных цифр), и даже GSL работает плохо с некоторыми входами в J_v. Обратите внимание, что при двойной проверке этих результатов были пересчитаны наихудшие показателиЦефеси случаи GSL с использованием функций.wolfram.com, и результат был проверен на основе наших тестовых данных: ошибок в тестовых данных обнаружено не было.

Implementation

Реализация в основном заключается в отфильтровке различных особых случаев:

Еслиxявляется отрицательным, то порядокvдолжен быть целым числом или результатом является ошибка домена. Если порядок является целым числом, то функция является нечетной для нечетных ордеров и даже для четных ордеров, поэтому мы отражаем вx >0.

Если порядокvотрицательный, то формулы отражения можно использовать для перехода кv >0:

Обратите внимание, что если порядок является целым числом, то эти формулы сводятся к:

J_-n= (-1)ⁿJ_n

Y_-n= (-1)ⁿY_n

Однако в целом отрицательный порядок подразумевает, что нам нужно будет вычислить и J, и Y.

Когдаxявляется большим по сравнению с порядкомv, то используются асимптотические расширения для большихxв М. Абрамовице и И. А. Стегуне,9.2.19 (они оказались более надежными, чем те, что в A&S 9.2.5).

Когда порядокvявляется целым числом, способ сначала связывает результат с J₀, J₁, Y₀и #160; и Y₁и #160; используя либо прямое, либо обратное повторение (алгоритм Миллера) в зависимости от того, какой из них стабилен. Значения для J₀, J₁, Y₀& #160; и Y₁& #160; вычисляются с использованием рациональных минимаксных приближений на интервалах корневых скобок для малых|x |и асимптотического расширения Ханкеля для больших|x |. Коэффициенты исходят из:

У. Джей КодиАлгоритм 715: SPECFUN - Портативный Фортран Package of Special Function Routines and Test Drivers, ACM Transactions on Mathematical Software, vol 19, 22 (1993).

J.F. Hart et al,Компьютерные приближения, Джон Wiley & Sons, Нью-Йорк, 1968.

Эти приближения точны примерно к 19 десятичным цифрам: поэтому эти методы не используются, когда тип T имеет более 64 двоичных цифр.

Когдаxменьше, чем машинный эпсилон, можно использовать следующие приближения для Y₀(x), Y₁(x), Y₂(x) и Y_n(x) (см.:http://functions.wolfram.com/03.03.06.0037.01,http://functions.wolfram.com/03.03.06.0038.01,http://functions.wolfram.com/03.03.06.0039.01иhttp://functions.wolfram.com/03.03.06.0040.01:

Когдаxневелик по сравнению сvиvне является целым числом, то для Y_v(x) может быть использовано следующее последовательное приближение, это также область, где другие приближения часто слишком медленные, чтобы сходиться для использования (см.http://functions.wolfram.com/03.03.06.0034.01):

Когдаxмало по сравнению сv, J_vx & #160; лучше всего вычисляется непосредственно из серии:

В общем случае мы вычисляем J_v& #160; и Y_v& #160; одновременно.

Чтобы получить начальные значения, пусть μ = ν - пол(ν + 1/2), затем μ является дробной частью ν так что |μ |<= 1/2 (это нам нужно для конвергенции позже). Идея состоит в том, чтобы вычислить J_{& #956;}(x), J_{& #956; +1}(x), Y_{& #956;}(x), Y_{& #956; +1}(x) и использовать их для получения J_{& #957;}(x), Y_{& #957;}(x).

Алгоритм называется методом Стида, который требует двух непрерывных фракций, а также Вронского:

См.: F.S. Acton,Numerical Methods that Work, The Mathematical Association of America, Washington, 1997.

Продолжающиеся фракции вычисляются с использованием модифицированного метода Ленца (W.J. Lentz,).Генерация функций Бесселя в расчетах рассеяния Ми с использованием непрерывных фракций, Applied Optics, vol 15, 668 (1976) . Их коэффициенты конвергенции зависят отx, поэтому нам нужны различные стратегии для большихxи малыхx.

x >v, CF1 нуждается в Ox) итерациях для конвергенции, CF2 быстро сходится

x<= v, CF1 быстро сходится, CF2 не сходится, когдаx<->>0

Когдаxявляется большимx>2), обе непрерывные фракции сходятся (CF1 может быть медленным для действительно большихx. J_{& #956;}, J_{& #956; +1}, Y_{& #956; Y_{& #956; +1}можно вычислить по}

где

J_{& #957;}и Y_{& #956;}затем вычисляются с использованием обратного (алгоритма Миллера) и обратного рецидива соответственно.

Когдаxневеликx<2, конвергенция CF2 может потерпеть неудачу (но CF1 работает очень хорошо). Решением здесь является серия Temme:

где

g_k& #160; и h_k& #160; также вычисляются рекурсиями (с участием гамма-функций), но формулы немного сложны, читатели ссылаются на Н.М. Темме,О численной оценке обычной функции Бесселя второго рода, Journal of Computational Physics, vol 21, 343 (1976). Примечание: Серия Темме сходится только для |μ |<= 1/2.

Как и в предыдущем случае, Y_{& #957;}& #160; вычисляется из форвардного повторения, так и Y_{& #957; +1}. С этими двумя значениями и f_{& #957;}, Вронский дает J_{& #957;}(x) непосредственно без обратного повторения.

Статья Bessel Functions of the First and Second Kinds раздела Math Toolkit 2.5.0 Bessel Functions может быть полезна для разработчиков на c++ и boost.

Bessel Functions of the First and Second Kinds

Boost , Math Toolkit 2.5.0 , Bessel Functions

Boost C++ Libraries

Bessel Functions of the First and Second Kinds

Synopsis

Description

Testing

Accuracy

Implementation