Рассчитайте функцию Bessel в MATLAB с использованием формулы JM+1 = 2MJ (M) -J (M -1)

Question

Я попытался реализовать функцию Бесселя, используя эту формулу, это код:

function result=Bessel(num);


if num==0
    result=bessel(0,1);
elseif num==1
    result=bessel(1,1);
else
    result=2*(num-1)*Bessel(num-1)-Bessel(num-2);
end;

Но если я использую функцию Bessel от Matlab, чтобы сравнить ее с этим, я получаю слишком высокие различные значения. Например, если я набираю Бессель (20), он дает мне 3.1689E+005 в результате, если вместо этого я введите Бессель (20,1), он дает мне 3.8735E-025, совершенно другой результат.

Answer 1

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

Рассмотрим следующее сравнение:

x = 0:20;
y1 = arrayfun(@(n)besselj(n,1), x);   %# builtin function
y2 = arrayfun(@Bessel, x);            %# your function
semilogy(x,y1, x,y2), grid on
legend('besselj','Bessel')
title('J_\nu(z)'), xlabel('\nu'), ylabel('log scale')

Таким образом, вы можете увидеть, как вычисленные значения начинают значительно различаться после 9.

Согласно Matlab:

Бессельдж использует интерфейс MEX в библиотеку Фортрана от De Amos.

и дает следующее в качестве ссылок для их реализации:

Де Амос, «Подпрограмма для функций Бесселя сложного аргумента и неотрицательного порядка», Национальный лабораторный отчет Сандии, Sand85-1018, май 1985 г.

Де Амос, «портативный пакет для функций Бесселя сложного аргумента и неотрицательного порядка», транс. Математика Программное обеспечение, 1986.

Answer 2

Переднее отношение рецидивов, которое вы используете, не является стабильным. Чтобы понять почему, учитывайте, что значения Besselj (n, x) становятся меньше и меньше примерно в факторе 1/2n. Анкет Вы можете увидеть это, посмотрев на первый срок серии Taylor для J.

Таким образом, то, что вы делаете, это вычитание большого количества из кратного несколько меньшего числа, чтобы получить еще меньшее число. Численно, это не будет работать хорошо.

Рассмотрим этот вариант. Мы знаем, что результат - порядок 10^-25. Анкет Вы начинаете с номеров, которые имеют порядок 1. Поэтому, чтобы получить из этого хотя бы одну точную цифру, мы должны знать первые два числа с точностью не менее 25 цифр. Мы явно этого не делаем, и рецидив на самом деле расходится.

Используя то же рецидивовое отношение, чтобы вернуться назад, от высоких заказов до низких заказов, является стабильный. Когда вы начинаете с правильных значений для J (20,1) и J (19,1), вы также можете рассчитать все заказы до 0 с полной точностью. Почему это работает? Потому что теперь числа становятся больше на каждом шаге. Вы вычитаете очень небольшое число из точного кратного большего числа, чтобы получить еще большее число.

Answer 3

Вы можете просто изменить код ниже, который предназначен для функции сферической Бесселя. Он хорошо протестирован и работает для всех аргументов и диапазона заказов. Мне жаль, что это в c#

     public static Complex bessel(int n, Complex z)
    {
        if (n == 0) return sin(z) / z;
        if (n == 1) return sin(z) / (z * z) - cos(z) / z;

        if (n <= System.Math.Abs(z.real))
        {
            Complex h0 = bessel(0, z);
            Complex h1 = bessel(1, z);
            Complex ret = 0;
            for (int i = 2; i <= n; i++)
            {
                ret = (2 * i - 1) / z * h1 - h0;
                h0 = h1;
                h1 = ret;
                if (double.IsInfinity(ret.real) || double.IsInfinity(ret.imag)) return double.PositiveInfinity;
            }
            return ret;
        }
        else
        {
            double u = 2.0 * abs(z.real) / (2 * n + 1);

            double a = 0.1;
            double b = 0.175;
            int v = n - (int)System.Math.Ceiling((System.Math.Log(0.5e-16 * (a + b * u * (2 - System.Math.Pow(u, 2)) / (1 - System.Math.Pow(u, 2))), 2)));

            Complex ret = 0;
            while (v > n - 1)
            {
                ret = z / (2 * v + 1.0 - z * ret);
                v = v - 1;
            }


            Complex jnM1 = ret;
            while (v > 0)
            {
                ret = z / (2 * v + 1.0 - z * ret);
                jnM1 = jnM1 * ret;
                v = v - 1;
            }
            return jnM1 * sin(z) / z;
        }
    }