Вычислите функцию Бесселя в MATLAB, используя формулу Jm+1=2mj(m) -j(m-1)

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

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, совершенно другой результат.


matlab recurrence math bessel-functions
person Ramy Al Zuhouri    schedule 18.10.2011    source источник
comment
Это может быть проблемой точности с плавающей запятой, когда вы пытаетесь сравнить два разных способа вычисления одного и того же значения. Я думаю, что функции Matlab Bessel используют mexfile Fortran, который может выполнять гораздо меньше вычислений, чем ваша реализация.   -  person Aabaz    schedule 19.10.2011
comment
Повторение Бесселя нестабильно, когда вы делаете это таким образом.   -  person Alexandre C.    schedule 27.10.2011

Ответы (3)


arrow_upward
3
arrow_downward

recurrence_bessel

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

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

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')

comparison_plot

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

Согласно МАТЛАБ:

BESSELJ использует интерфейс MEX для библиотеки Fortran от D. E. Amos.

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

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

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

person Amro    schedule 27.10.2011

arrow_upward
2
arrow_downward

Используемое вами прямое рекурсивное отношение не является стабильным. Чтобы понять почему, рассмотрим, что значения BesselJ(n,x) становятся все меньше и меньше примерно в 1/2n раз. Вы можете убедиться в этом, взглянув на первый член ряда Тейлора для Дж.

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

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

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

person Jeffrey Sax    schedule 27.10.2011
comment
Решение состоит в том, чтобы использовать алгоритм Миллера и тот факт, что повторение действительно стабильно, когда n ‹ x . Когда n > x, вы выполняете рекурсию в обратном направлении, пока не достигнете J_0 или J_1, и вы найдете коэффициент нормализации, вычислив J_0 или J_1 в этой точке. - person Alexandre C.; 27.10.2011

arrow_upward
0
arrow_downward

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

     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;
        }
    }
person falopsy    schedule 12.03.2013