Блог

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

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

#math #matlab #recurrence #bessel-functions

#математика #matlab #повторение #функции Бесселя

Вопрос:

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

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

Комментарии:

1. Это может быть проблемой точности с плавающей запятой, поскольку вы пытаетесь сравнить два разных способа вычисления одного и того же значения. Я думаю, что функции Бесселя Matlab используют mexfile Fortran, который может выполнять намного меньше вычислений, чем ваша реализация.

2. Повторение Бесселя нестабильно, когда вы делаете это таким образом.

Ответ №1:

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.

Согласно MATLAB:

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

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

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

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

Ответ №2:

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

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

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

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

Комментарии:

1. Решение состоит в том, чтобы использовать алгоритм Миллера и тот факт, что повторение фактически стабильно при n < x. Когда n > x, вы выполняете рекурсию назад, пока не достигнете J_0 или J_1, и вы находите коэффициент нормализации, вычисляя J_0 или J_1 в этот момент.

Ответ №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;
        }
    }