Рассчитать функцию Бесселя в 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;

Но если я использую функцию бесселя MATLAB, чтобы сравнить ее с этой, я получу слишком высокие значения. Например, если я набираю Bessel (20), это дает мне 3.1689e + 005 в качестве результата, если вместо этого я набираю bessel (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:

BESSELJ используетMEX-интерфейс с библиотекой Fortran от DE Amos.

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

DE Amos, "Пакет подпрограмм для функций Бесселякомплексный аргумент и неотрицательный порядок ", Отчет Национальной лаборатории Сандии, SAND85-1018, май, 1985.

DE Амос," ​​Переносимый пакет для функций Бесселя комплексного аргумента и неотрицательного порядка ", Пер.МатематикаПрограммное обеспечение, 1986.

Answer 2

Используемое вами прямое рекуррентное отношение нестабильно.Чтобы понять почему, рассмотрим, что значения BesselJ (n, x) становятся все меньше и меньше примерно в 100 раз *.Вы можете увидеть это, посмотрев на первый член ряда Тейлора для 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;
        }
    }