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

Finding Zeros of Bessel Functions of the First and Second Kinds

Synopsis

Включить<Усиление/Математика/Специальные функции/Бессель.hpp>

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

Сигнатурой единичных функций значений являются:

template <class T>
T cyl_bessel_j_zero(
         T v,            // Floating-point value for Jv.
         int m);         // 1-based index of zero.
template <class T>
T cyl_neumann_zero(
         T v,            // Floating-point value for Jv.
         int m);         // 1-based index of zero.

Для нескольких нулей:

template <class T, class OutputIterator>
OutputIterator cyl_bessel_j_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of first zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it);    // Destination for zeros.
template <class T, class OutputIterator>
OutputIterator cyl_neumann_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of zero.
                     unsigned number_of_zeros,  // How many zeros to generate
                     OutputIterator out_it);    // Destination for zeros.

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

 template <class T>
 T cyl_bessel_j_zero(
          T v,            // Floating-point value for Jv.
          int m,          // 1-based index of zero.
          const Policy&); // Policy to use.
 template <class T>
 T cyl_neumann_zero(
          T v,            // Floating-point value for Jv.
          int m,          // 1-based index of zero.
          const Policy&); // Policy to use.
template <class T, class OutputIterator>
OutputIterator cyl_bessel_j_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of first zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it,     // Destination for zeros.
                     const Policy& pol);        // Policy to use.
template <class T, class OutputIterator>
OutputIterator cyl_neumann_zero(
                     T v,                       // Floating-point value for Jv.
                     int start_index,           // 1-based index of zero.
                     unsigned number_of_zeros,  // How many zeros to generate.
                     OutputIterator out_it,     // Destination for zeros.
                     const Policy& pol);        // Policy to use.

Каждый реальный порядок ν цилиндрические функции Бесселя и Неймана имеют бесконечное число нулей на положительной реальной оси. Реальные нули на положительной реальной оси можно найти, решая для корней

J_{& #957;}(j_{& #957;, m}) = 0

Y_{& #957;}(y_{& #957;, m}) = 0

Здесьj_{и #957;, m}представляетm^thкорень цилиндрической функции Бесселя порядкаи #957;, иу_{и #957;, m}представляетm^thкорень цилиндрической функции Неймана порядкаи #957;.

Нули или корни (значенияx, где функция пересекает горизонтальнуюу=0ось) функций Бесселя и Неймана вычисляются двумя функциями,cyl_bessel_j_zeroиcyl_neumann_zero.

В каждом случае индекс или ранг нуля возвращается на основе 1, что означает:

cyl_bessel_j_zero(v, 1);

Возвращает первый ноль Бесселя Дж.

Прохождениеstart_index<=0приводит к тому, чтоstd::домен_errorподнимается.

Для некоторых параметров, однако, корень нуля определяется и имеет значение нуля. Например, нулевой кореньJподv]xопределен и имеет значение нуля для всех значенийvи0и для отрицательных целых значенийv=-n. Аналогичные случаи описаны в деталях осуществления ниже.

ПорядокvJможет быть положительным, отрицательным и нулевым для функцийcyl_bessel_jиcyl_neumann, но не бесконечным и не NaN.

Examples of finding Bessel and Neumann zeros

Этот пример демонстрирует вычисление нулей функций Бесселя и Неймана. Это также показывает, как увеличить. Математика и рост. Многоточность может быть объединена, чтобы обеспечить точность многих десятичных цифр. Для точности 50 десятичных цифр мы должны включить

#include <boost/multiprecision/cpp_dec_float.hpp>

иtypedefдляfloat_typeможет быть удобным (что позволяет быстро переключаться на повторные вычисления при встроенномdoubleили другой точности)

typedef boost::multiprecision::cpp_dec_float_50 float_type;

Чтобы использовать функции для поиска нулей функций, которые нам нужны

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

Этот файл включает в себя форвардные подписи декларации для функций нулевого поиска:

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

Более подробная информация содержится в полной документации, например, вBoost. Математические функции Бесселя.

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

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

Во-первых, оцените один нулевой Бессель.

Точность контролируется параметром шаблона типа точки поплавкаTv, поэтому этот пример имеетдвойнуюточность, по меньшей мере 15, но до 17 десятичных цифр (для обычного 64-битного двойника).

//    double root = boost::math::cyl_bessel_j_zero(0.0, 1);
//    // Displaying with default precision of 6 decimal digits:
//    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40483
//    // And with all the guaranteed (15) digits:
//    std::cout.precision(std::numeric_limits<double>::digits10);
//    std::cout << "boost::math::cyl_bessel_j_zero(0.0, 1) " << root << std::endl; // 2.40482555769577

Но обратите внимание, что поскольку параметрvконтролирует точность результата,vдолжен быть типом с плавающей точкой. Таким образом, если вы предоставляете целочисленный тип, скажем 0, а не 0,0, то он не будет компилироваться таким образом:

root = boost::math::cyl_bessel_j_zero(0, 1);

с сообщением об ошибке

error C2338: Order must be a floating-point type.

Мы можем использовать политику, чтобы игнорировать ошибки, C-стиль, возвращая некоторую ценность, возможно, бесконечность или NaN, или лучшее, что можно сделать. (См.Обработка ошибок пользователя).

Создать (возможно, неразумную!) политикуигнорирование_все_политики, которая игнорирует все ошибки:

typedef boost::math::policies::policy<
  boost::math::policies::domain_error<boost::math::policies::ignore_error>,
  boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
  boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
  boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
  boost::math::policies::pole_error<boost::math::policies::ignore_error>,
  boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
            > ignore_all_policy;

Примеры использования этогоигнорируют_all_policy

double inf = std::numeric_limits<double>::infinity();
double nan = std::numeric_limits<double>::quiet_NaN();
double dodgy_root = boost::math::cyl_bessel_j_zero(-1.0, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(-1.0, 1) " << dodgy_root << std::endl; // 1.#QNAN
double inf_root = boost::math::cyl_bessel_j_zero(inf, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl; // 1.#QNAN
double nan_root = boost::math::cyl_bessel_j_zero(nan, 1, ignore_all_policy());
std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl; // 1.#QNAN

Другая версияcyl_bessel_j_zeroпозволяет вычислить несколько нулей одним вызовом, помещая результаты в контейнер, частоstd::вектор. Например, сгенерировать и отобразить первые пятьдвойныхкорней J_vдля интегрального порядка 2, как столбецJ₂(x)в таблице 1Wolfram Bessel Function Zeros.

unsigned int n_roots = 5U;
std::vector<double> roots;
boost::math::cyl_bessel_j_zero(2.0, 1, n_roots, std::back_inserter(roots));
std::copy(roots.begin(),
          roots.end(),
          std::ostream_iterator<double>(std::cout, "\n"));

Или мы можем использовать Boost. Многоточность для генерации 50 десятичных цифровых корнейJ_vдля неинтегрального порядкаv=71/19==3,736842, выраженная в виде точной целочисленной фракции для получения максимально точного значения, возможного для всех типов с плавающей точкой.

Мы установили точность выходного потока и показали нулевые значения для отображения фиксированных 50 десятичных цифр.

std::cout.precision(std::numeric_limits<float_type>::digits10); // 50 decimal digits.
std::cout << std::showpoint << std::endl; // Show trailing zeros.
float_type x = float_type(71) / 19;
float_type r = boost::math::cyl_bessel_j_zero(x, 1); // 1st root.
std::cout << "x = " << x << ", r = " << r << std::endl;
r = boost::math::cyl_bessel_j_zero(x, 20U); // 20th root.
std::cout << "x = " << x << ", r = " << r << std::endl;
std::vector<float_type> zeros;
boost::math::cyl_bessel_j_zero(x, 1, 3, std::back_inserter(zeros));
std::cout << "cyl_bessel_j_zeros" << std::endl;
// Print the roots to the output stream.
std::copy(zeros.begin(), zeros.end(),
          std::ostream_iterator<float_type>(std::cout, "\n"));

Using Output Iterator to sum zeros of Bessel Functions

Этот пример показывает суммирование нулей функций Бесселя. Чтобы использовать функции для поиска нулей функций, которые нам нужны

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

Мы используемcyl_bessel_j_zeroпараметр итератора выводаout_itдля создания суммы1/zeros², определяя пользовательский итератор вывода:

template <class T>
struct output_summation_iterator
{
   output_summation_iterator(T* p) : p_sum(p)
   {}
   output_summation_iterator& operator*()
   { return *this; }
    output_summation_iterator& operator++()
   { return *this; }
   output_summation_iterator& operator++(int)
   { return *this; }
   output_summation_iterator& operator = (T const& val)
   {
     *p_sum += 1./ (val * val); // Summing 1/zero^2.
     return *this;
   }
private:
   T* p_sum;
};

Сумма рассчитывается для многих значений, сходящихся на аналитическое точное значение1/8.

using boost::math::cyl_bessel_j_zero;
double nu = 1.;
double sum = 0;
output_summation_iterator<double> it(&sum);  // sum of 1/zeros^2
cyl_bessel_j_zero(nu, 1, 10000, it);
double s = 1/(4 * (nu + 1)); // 0.125 = 1/8 is exact analytical solution.
std::cout << std::setprecision(6) << "nu = " << nu << ", sum = " << sum
  << ", exact = " << s << std::endl;
// nu = 1.00000, sum = 0.124990, exact = 0.125000

Calculating zeros of the Neumann function.

Этот пример также показывает, как увеличить. Математика и рост. Многоточность может быть объединена, чтобы обеспечить точность многих десятичных цифр. Для точности 50 десятичных цифр мы должны включить

#include <boost/multiprecision/cpp_dec_float.hpp>

typedef boost::multiprecision::cpp_dec_float_50 float_type;

Чтобы использовать функции для нахождения нулей функцииcyl_neumann, нам нужно:

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

Нули функции Неймана (Бесселя Y) оцениваются примерно одинаково:

using boost::math::cyl_neumann_zero;
double zn = cyl_neumann_zero(2., 1);
std::cout << "cyl_neumann_zero(2., 1) = " << zn << std::endl;
std::vector<float> nzeros(3); // Space for 3 zeros.
cyl_neumann_zero<float>(2.F, 1, nzeros.size(), nzeros.begin());
std::cout << "cyl_neumann_zero<float>(2.F, 1, ";
// Print the zeros to the output stream.
std::copy(nzeros.begin(), nzeros.end(),
          std::ostream_iterator<float>(std::cout, ", "));
std::cout << "\n""cyl_neumann_zero(static_cast<float_type>(220)/100, 1) = "
  << cyl_neumann_zero(static_cast<float_type>(220)/100, 1) << std::endl;
// 3.6154383428745996706772556069431792744372398748422

Error messages from 'bad' input

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

Чтобы использовать функции для нахождения нулей нужных нам функций:

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

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

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

try
{ // Try a zero order v.
  float dodgy_root = boost::math::cyl_bessel_j_zero(0.F, 0);
  std::cout << "boost::math::cyl_bessel_j_zero(0.F, 0) " << dodgy_root << std::endl;
  // Thrown exception Error in function boost::math::cyl_bessel_j_zero<double>(double, int):
  // Requested the 0'th zero of J0, but the rank must be > 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}

	Note
	Тип, показанный в сообщении об ошибке, - это типпосле продвижения, использующийполитику точностиивнутреннюю политику продвижения, от`плавающего`до`двойного`в этом случае.

В этом примере продвижение идет:

Аргументыплаваютиint.
int«как если бы» это былдвойной, поэтому аргументыплаваютидвойной.
Обычный типдвойной- так что это точность, которую мы хотим (и тип, который будет возвращен).
Оцените внутренне какдвойнойдля полнойплавающейточности.

См. полный код для других примеров, которые способствуют отдвойнойдодлиннойдвойной.

Другие примеры «плохих» входов, таких как бесконечность и NaN, приведены ниже. Некоторые предупреждения компилятора указывают на то, что во время компиляции обнаруживаются «плохие» значения.

try
{ // order v = inf
   std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << std::endl;
   double inf = std::numeric_limits<double>::infinity();
   double inf_root = boost::math::cyl_bessel_j_zero(inf, 1);
   std::cout << "boost::math::cyl_bessel_j_zero(inf, 1) " << inf_root << std::endl;
   // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
   // Order argument is 1.#INF, but must be finite >= 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}
try
{ // order v = NaN, rank m = 1
   std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << std::endl;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double nan_root = boost::math::cyl_bessel_j_zero(nan, 1);
   std::cout << "boost::math::cyl_bessel_j_zero(nan, 1) " << nan_root << std::endl;
   // Throw exception Error in function boost::math::cyl_bessel_j_zero<long double>(long double, unsigned):
   // Order argument is 1.#QNAN, but must be finite >= 0 !
}
catch (std::exception& ex)
{
  std::cout << "Thrown exception " << ex.what() << std::endl;
}

Выводы из других примеров показаны приложенными к полному списку кода.

Полный код (и выход) для этих примеров находится на. Бесселевые нули,Бесселевые нули итератор,Неймановы нули,Бесселевские сообщения об ошибках.

Implementation

Различные методы используются для вычисления начальных оценок дляj_{& #957;, m}иy_{& #957;, m}; подробно описаны ниже.

После нахождения первоначальной оценки данного корня его точность впоследствии уточняется до желаемого уровня с помощью итерации Ньютона-Рафсона от Boost. Математическийкорневой поиск с производнымиутилиты в сочетании с функциямиcyl_bessel_jиcyl_neumann.

Итерация Ньютона требует какJ_{& #957;}(x)илиY_{& #957;}(x), так и его производного. ПроизводныеJ_{& #957;}(x)иY_{& #957;}(x)в отношенииxданы М. Абрамовичем и И. А. Стегуном, Справочник математических функций, NBS (1964). В частности,

d/_dxJ_{& #957;}(x)=J_{& #957;}/ x

& #8193;d/_dxY_{& #957;}(x)=Y_{& #957; & #957; Y_{& #957;}}/ x

Перечисление ранга корня (иными словами, индекс корня) начинается с одного и считается, другими словами.m,=1,2,3,…Значение первого корня всегда больше нуля.

Для некоторых особых параметров цилиндрические функции Бесселя и цилиндрические функции Неймана имеют корень в происхождении. Например,J_{& #957;}[x]имеет корень в начале для каждого положительного порядка& #957; >0, и для каждого отрицательного целого порядка& #957; = -nсn & #8712; & #8469;⁺иn & #8800; 0.

Кроме того,Y_{& #957;}[x]имеет корень в начале для каждого отрицательного полуцелого порядка& #957; = -n/2, сn & #8712; & #8469;⁺иn & #8800; 0.

Для этих специальных значений параметров происхождение со значениемx = 0обеспечивается в виде0^thкорня, генерируемогоcyl_bessel_j_zeroиcyl_neumann_zero.

При вычислении исходных оценок для корней функций Бесселя проводится различие между положительным и отрицательным порядком, и для них используются разные методы. Кроме того, различные алгоритмы используются для первого корням = 1и для последующих корней с более высоким рангомм ≥ 2. Кроме того, оценки корней для функций Бесселя с порядком выше и ниже отсечки прии #957; = 2.2вычисляются различными методами.

Расчеты оценокj_{& #957;,1}иy_{& #957;,1}с0 & #8804; & #957;< 2.2используют эмпирически табулированные значения. Коэффициенты для них генерируются компьютерной системой алгебры.

Расчеты оценокj_{& #957;,1}иy_{& #957;,1}с& #957; & #8805; 2.2используют Eqs.9.5.14 и 9.5.15 в M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS (1964).

В частности,

& #8193;j_{& #957;,1}& #8773; & #957; + 1.85575 & #957;^{& #8531;}+ 1.033150 & #957;^{- & #8531;}- 0,00397 & #957;^{- -1}- 0,0908 & #957;^-5/3+ 0,043 & #957;^-7/3+ & #8230;

и

y_{& #957;,1}& #8773; & #957; + 0.93157 & #957;^{& #8531;}+ 0.26035 & #957;^-⅓+ 0,01198 & #957;^-1- 0,0060 & #957;^-5/3- 0,001 & #957;^-7/3+ & #8230;

Расчеты оценокj_{& #957;, m}иy_{& #957;, m}с рангомm >2и0 & #8804; & #957;< 2.2используют приближение МакМахона, как описано в M. Abramowitz и I. A. Stegan, Section 9.5 and 9.5.12. В частности,

j_{& #957;,m},y_{& #957;,m}& #8773; & #946; - (μ-1) / 8β

- 4(μ-1)(7μ - 31) / 3(8β)³

-32(μ-1)(83μ² - 982μ +3779) / 15(8β)⁵

-64(μ-1)(6949μ³- 153855μ² + 1585743μ- 6277237) / 105(8a)⁷

& #8193; & #8193; & #8193;- & #8230;& #8193; & #8193; (5)

где& #956; = 4ν²и& #946; = (m + & #189; & #957; - & #188; - & #660;дляj_{& #957;,m}и& #946; = (m + & #189; & #957; -¾)π дляy_{& #957;,m}.

Расчеты оценокj_{& #957;, m}иy_{& #957;, m}с& #957; & #8805; 2.2используют один термин в асимптотическом расширении, данном в Eq.9.5.22 и верхней строке Eq.9.5.26 в сочетании с Eq. 9.3.39, все в M. Abramowitz и I. A. Stegun, Handbook of Mathematical Functions, NBS (1964) явное и простое для понимания лечение асимптотического расширения нулей. Последние два уравнения выражены для аргумента(x)больше одного. (Олвер также дает серию уравнений в& #167; 10.21(vi) Макмэхон асимптотические расширения для больших нулей- с использованием немного разных переменных имен.)

Вкратце,

j_{& #957;, m}& #8764; & #957;x(-ζ) + f₁(-ζ/ν)

где-ζ = ν^-2/3a_mиa_mявляется абсолютным значениемm^thкорняAi(x)на отрицательной реальной оси.

Здесьx = x(-ζ)обратная функция

⅔(-ζ)^3/2= √(x² - 1) - cos⁻¹(1/x) (7)

Кроме того,

f₁(-ζ) = ½x(-ζ) {h(-ζ)}² ⋅ b₀(-ζ)

где

h(-ζ) = {4(-ζ) / (x² - 1)}⁴

и

b₀(-ζ) = -5/(48ζ²) + 1/(-ζ)^½⋅ {5/(24(x²-1)^3/2) + 1/(8(x²-1)^½]}

При решении дляx(-ζ)в экв. 7 выше, правая сторона расширена до порядка 2 в серии Тейлора для большихx. Это приводит к

⅔(-ζ)^3/2≈ x + 1/2x - π/2

Положительный корень полученного квадратичного уравнения используется для нахождения начальной оценкиx(-ζ). Эта первоначальная оценка впоследствии уточняется несколькими этапами итерации Ньютона-Рафсона в экв. 7.

Оценки корней цилиндрических Бесселевых функций отрицательного порядка на положительной реальной оси обнаруживаются с помощью переплетения отношений. Например, кореньm^thцилиндрической функции Бесселяj_-ν,mзаключен в скобки корнемm^thи корнем[m+1]^thфункции Бесселя соответствующего положительного целого порядка. Другими словами,

j_nν,m<j_-ν,m<j_nν,m+1

гдеm >1иn_{& #957;}представляет собой интегральный пол абсолютного значения|-ν |.

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

Соотношения скобки не удерживают первый корень цилиндрических функций Бесселя и цилиндрических функций Неймана с отрицательным порядком. Поэтому итеративные алгоритмы в сочетании с корневым поиском через бисекция используются для локализацииj_-ν,1иy_-ν,1.

Testing

Точность оценки нулей была проверена на 50 десятичных цифрах с использованиемcpp_dec_float_50и найдена идентичной с точечными значениями, вычисленнымиWolfram Alpha.

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

Finding Zeros of Bessel Functions of the First and Second Kinds

Boost , Math Toolkit 2.5.0 , Bessel Functions

Boost C++ Libraries

Finding Zeros of Bessel Functions of the First and Second Kinds

Synopsis

Description

Examples of finding Bessel and Neumann zeros

Using Output Iterator to sum zeros of Bessel Functions

Calculating zeros of the Neumann function.

Error messages from 'bad' input

Implementation

Testing