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

Modified 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_i(T1 v, T2 x);
template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_i(T1 v, T2 x, const Policy&);
template <class T1, class T2>
calculated-result-type cyl_bessel_k(T1 v, T2 x);
template <class T1, class T2, class Policy>
calculated-result-type cyl_bessel_k(T1 v, T2 x, const Policy&);

Description

Функцииcyl_bessel_iиcyl_bessel_kвозвращают результат модифицированных функций Бесселя первого и второго рода соответственно:

cyl_bessel_i(v, x) = I_v(x)

cyl_bessel_k(v, x) = K_v(x)

где:

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

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

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

Следующий график иллюстрирует экспоненциальное поведение I_v.

Следующий график иллюстрирует экспоненциальный распад K_v.

Testing

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

Accuracy

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

Table 6.44. Error rates for cyl_bessel_i (integer orders)

	Microsoft Visual C++ версия 12.0 Win32	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux double	Солнечный компилятор версии 0x5130 Солнечный солярис
Бессел I0: Mathworld Data (целая версия)	Макс = 0,877ε (Средний = 0,549ε)	Макс = 4.57ε (Средний = 2.1ε) :Макс = 8.49ε (Средний = 3.46ε)И другие неудачи.	Max = 0ε (Mean = 0ε) GSL 1.16:Max = 0.79ε (Mean = 0.482ε)] Rmath 3.0.2:Макс = +INFε (Средний = +INFε)И другие сбои.) Цефес:Макс = 2,55e+43ε (Средний = 8,06e+42ε)И другие неудачи.	Макс = 4.54ε (Средний = 2.11ε)
Бессель I1: Mathworld Data (целая версия)	Макс = 0.885ε (Средний = 0.55ε)	Макс = 7.83ε (Средний = 2.79ε) :Макс = 5ε (Средний = 2.15ε)И другие неудачи.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 0.82ε (Mean = 0.456ε)] Rmath 3.0.2:Макс = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 1.28e+43ε (Средний = 4.05e+42ε)И другие сбои.	Макс = 6,52ε (Средний = 2,25ε)
Бессел In: Mathworld Data (целая версия)	Макс = 3.46ε (Средний = 1.32ε)	Max = 1.8ε (Mean = 1.33ε) :Max = 430ε (Mean = 163ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 5.15ε (Mean = 2.13ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 3,67e+177ε (Средний = +INFε)И другие неудачи.	Макс = 463ε (Средний = 140ε)

Table 6.45. Error rates for cyl_bessel_i

	Microsoft Visual C++ версия 12.0 Win32	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux double	Солнечный компилятор версии 0x5130 Солнечный солярис
Бессел I0: Mathworld Data	Макс = 0,877ε (Средний = 0,549ε)	Max = 4.57ε (Mean = 2.1ε) :Max = 8.49ε (Mean = 3.46ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 270ε (Mean = 91.6ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 2,55e+43ε (Средний = 8,06e+42ε)И другие сбои.	Макс = 4.54ε (Средний = 2.11ε)
Бессель I1: Mathworld Data	Макс = 0.885ε (Средний = 0.55ε)	Max = 7.83ε (Mean = 2.79ε) :Max = 5ε (Mean = 2.15ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 128ε (Mean = 41ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 1.28e+43ε (Средний = 4.05e+42ε)И другие сбои.	Макс = 6,52ε (Средний = 2,25ε)
Бессел In: Mathworld Data	Макс = 3.46ε (Средний = 1.32ε)	Max = 1.8ε (Mean = 1.33ε) :Max = 430ε (Mean = 163ε)И другие сбои.	Max = 0ε (Mean = 0ε) GSL 1.16:Max = 2.31ε (Mean = 0.838ε)И другие сбои.) Rmath 3.0.2:Max = +INFε (Mean = +INFε)И другие неудачи.) Цефес:Макс = 3,67e+177ε (Средний = +INFε)И другие неудачи.	Макс = 463ε (Средний = 140ε)
Бессел Ив: Данные математического мира	Макс = 2.97ε (Средний = 1.33ε)	Max = 4.12ε (Mean = 1.85ε) :Max = 616ε (Mean = 221ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 5.95ε (Mean = 2.08ε)И другие сбои.) Rmath 3.0.2:Max = 1e+04ε (Mean = 3.18e+03ε) Cephes:Max = +INFε (Mean = +INFε)И другие неудачи.	Макс = 4.12ε (Средний = 1,95ε)
Бессел В: Случайные данные	Макс = 9,67ε (Средний = 1,89ε)	Макс = 6.79ε (Mean = 1.15ε) :Макс = 645ε (Mean = 132ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 261ε (Mean = 53.2ε)И другие неудачи.) Rmath 3.0.2:Max = 7.37ε (Mean = 2.4ε)] Цефес:Max = 4.22e+06ε (Mean = 2.26e+05ε)	Макс = 176ε (Средний = 39.2ε)
Бессель Ив: Случайные данные	Макс = 7,46ε (Средний = 1,54ε)	Макс = 8.35ε (Средний = 1,49ε) :Макс = 1,05e+03ε (Средний = 224ε)И другие сбои.	Макс = 0.661ε (Mean = 0.0441ε) GSL 1.16:Макс = 6.18e+03ε (Mean = 1.55e+03ε)И другие сбои.) Rmath 3.0.2:Max = 4.28e+08ε (Mean = 2.85e+07ε) Цефес:Max = 6e+30ε (Mean = 4e+29ε)И другие сбои.	Макс = 283ε (Средний = 88.4ε)
Бессель Ив: Данные математического мира (большие значения)	Макс = 3,67ε (Средний = 1,64ε)	Max = 14.7ε (Mean = 6.57ε) :Max = 118ε (Mean = 57.2ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 37ε (Mean = 18ε)И другие сбои.) Rmath 3.0.2:Max = 3,77e+168ε (Mean = 2,39e+168ε)И другие сбои.) Цефес:Макс = 73,7ε (Средний = 58,5ε)	Макс = 14,7ε (Средний = 6,59ε)

Table 6.46. Error rates for cyl_bessel_k (integer orders)

	Microsoft Visual C++ версия 12.0 Win32	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux double	Солнечный компилятор версии 0x5130 Солнечный солярис
Бессел K0: Данные Mathworld (целая версия)	Макс = 1.55ε (Средний = 0.837ε)	Макс = 4.16ε (Средний = 1,46ε) :Макс = 9,33ε (Средний = 3,25ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 1.2ε (Mean = 0.733ε)] Rmath 3.0.2:Макс = 0.833ε (Mean = 0.601ε)] Цефес:Макс = 1.1e+06ε (Mean = 3.68e+05ε)	Макс = 4.16ε (Средний = 1,49ε)
Бессел К1: Данные Mathworld (целая версия)	Макс = 1ε (Средний = 0,573ε)	Max = 1.8ε (Mean = 1.02ε) :Max = 8.94ε (Mean = 3.19ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 0.626ε (Mean = 0.333ε)] Rmath 3.0.2:Макс = 0.894ε (Mean = 0.516ε)] Цефес:Макс = 5.38e+05ε (Mean = 1.79e+05ε)	Макс = 1,8ε (Средний = 1,02ε)
Бессел Кн: Данные Mathworld (целая версия)	Макс = 3,63ε (Средний = 1,46ε)	Max = 4.48ε (Mean = 2.14ε) :Max = 12.9ε (Mean = 4.91ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 168ε (Mean = 59.5ε)] Rmath 3.0.2:Макс = 8.48ε (Mean = 2.98ε)] Цефес:Макс = +INFε (Mean = +INFε)И другие неудачи.	Макс = 4,48ε (Средний = 1,98ε)

Table 6.47. Error rates for cyl_bessel_k

	Microsoft Visual C++ версия 12.0 Win32	GNU C++ версия 5.1.0 Linux длинный двойной	GNU C++ версия 5.1.0 Linux double	Солнечный компилятор версии 0x5130 Солнечный солярис
Бессел К0: Данные математического мира	Макс = 1.55ε (Средний = 0.837ε)	Макс = 4.16ε (Средний = 1,46ε) :Макс = 9,33ε (Средний = 3,25ε)]	Max = 0ε (Mean = 0ε) GSL 1.16:Max = 6.04ε (Mean = 2.16ε)] Rmath 3.0.2:Max = 0.833ε (Mean = 0.601ε)]	Макс = 4.16ε (Средний = 1,49ε)
Бессел К1: Данные математического мира	Макс = 1ε (Средний = 0,573ε)	Max = 1.8ε (Mean = 1.02ε) :Max = 8.94ε (Mean = 3.19ε)]	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 6.26ε (Mean = 2.21ε)] Rmath 3.0.2:Макс = 0.894ε (Mean = 0,516ε)]	Макс = 1,8ε (Средний = 1,02ε)
Бессел Кн: Данные математического мира	Макс = 3,63ε (Средний = 1,46ε)	Max = 4.48ε (Mean = 2.14ε) :Max = 12.9ε (Mean = 4.91ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 3.36ε (Mean = 1,43ε)И другие сбои.) Rmath 3.0.2:Max = 8.48ε (Mean = 2.98ε)	Макс = 4,48ε (Средний = 1,98ε)
Бессел Кв: Данные математического мира	Макс = 4.78ε (Средний = 2.2ε)	Макс = 3.58ε (Средний = 2.44ε) :Макс = 13ε (Средний = 4.81ε)И другие неудачи.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 5.47ε (Mean = 2.04ε)И другие сбои.) Rmath 3.0.2:Max = 3.15ε (Mean = 1.35ε)	Макс = 3,58ε (Средний = 2,29ε)
Бессел Кв: Данные математического мира (большие значения)	Макс = 59,8ε (Средний = 26,9ε)	Max = 42.3ε (Mean = 21ε) :Max = 42.3ε (Mean = 19.8ε)И другие сбои.	Макс = 0ε (Mean = 0ε) GSL 1.16:Макс = 308ε (Mean = 142ε)И другие сбои.) Rmath 3.0.2:Max = 84.6ε (Mean = 37.8ε)	Макс = 43.1ε (Средний = 21.3ε)
Бессел Кн: Случайные данные	Макс = 7,47ε (Средний = 1,4ε)	Max = 4.55ε (Mean = 1.09ε) :Max = 13.9ε (Mean = 2.91ε)]	Макс = 0,764ε (Mean = 0,0348ε) GSL 1.16:Макс = 8,71ε (Mean = 1,76ε)И другие сбои.) Rmath 3.0.2:Max = 7.47ε (Mean = 1.34ε)	Макс = 4.55ε (Средний = 1.21ε)
Бессел Кв: Случайные данные	Макс = 8,33ε (Средний = 1,62ε)	Max = 7.88ε (Mean = 1.48ε) :Max = 13.6ε (Mean = 2.68ε)И другие сбои.	Макс = 0,507ε (Mean = 0,0313ε) GSL 1.16:Макс = 9.71ε (Mean = 1,47ε)И другие сбои.) Rmath 3.0.2:Max = 7.37ε (Mean = 1,49ε)	Макс = 7,88ε (Средний = 1,49ε)

Implementation

В первую очередь рассматриваются следующие особые случаи:

При вычислении I_v дляx< 0, то ν должно быть целое число или возникает ошибка домена. Если ν является целым числом, то функция нечетная, если ν нечетная и даже если ν четная, и мы можем отражать доx >0.

Для I_v& #160; с v, равным 0, 1 или 0,5, рассматриваются как особые случаи.

В 0 и 1 случаях используются минимаксные рациональные приближения на конечных и бесконечных интервалах. Коэффициенты исходят из:

Дж.М. Блэр и К.А. ЭдвардсСтабильные рациональные приближения к модифицированным функциям Бесселя I_0(x) и I_1(x), Atomic Energy of Canada Limited Report 4928, Chalk River, 1974.
С. МоширМетоды и программы для математических функций, Ellis Horwood Ltd, Chichester, 1989.

В то время как корпус 0.5 является простой тригонометрической функцией:

I_0,5(x) = sqrt(2 / & #960;x) * sinh(x)

Для K_v& #160; сvцелым числом результат рассчитывается с использованием отношения повторения:

начиная с K₀& #160; и K₁& #160; которые рассчитываются с использованием рациональных приближений выше. Эти рациональные приближения точны примерно до 19 цифр и поэтому используются только тогда, когда T имеет не более 64 двоичных цифр точности.

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

В общем случае мы сначала нормализуем ν до<0, [inf]>с помощью формул отражения:

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

Алгоритм предложен Темме в Н.М. Темме. О численной оценке модифицированной функции осей третьего рода, Journal of Computational Physics, vol 19, 324 (1975), которая требует двух непрерывных фракций, а также Вронского:

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

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

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

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

где

Sтакже представляет собой серию, которая суммируется вместе с CF2, см. I.J. Thompson and A.R. Barnett,Модифицированные функции Бесселя I_v и K_v реального порядка и сложного аргумента с выбранной точностью, Computer Physics Communications, vol 47, 245 (1987).

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

где

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

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

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

Modified Bessel Functions of the First and Second Kinds

Boost , Math Toolkit 2.5.0 , Bessel Functions

Boost C++ Libraries

Modified Bessel Functions of the First and Second Kinds

Synopsis

Description

Testing

Accuracy

Implementation