Calcule la función de Bessel en MATLAB usando JM+1 = 2MJ (M) -J (M -1) Fórmula

Question

Traté de implementar la función Bessel usando esa fórmula, este es el código:

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;

Pero si uso la función Bessel de Matlab para compararla con este, obtengo valores demasiado altos. Por ejemplo, si escribo Bessel (20) me da 3.1689e+005 como resultado, si en su lugar escribo Bessel (20,1) me da 3.8735e-025, un resultado totalmente diferente.

Answer 1

Tales relaciones de recurrencia son buenas en matemáticas pero numéricamente inestables al implementar algoritmos utilizando representaciones de precisión limitadas de números de punto flotante.

Considere la siguiente comparación:

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

Entonces puede ver cómo los valores calculados comienzan a diferir significativamente después de 9.

Según Matlab:

Besselj utiliza una interfaz MEX para una biblioteca Fortran de De Amos.

y da lo siguiente como referencias para su implementación:

De Amos, "Un paquete de subrutina para funciones de Bessel de un argumento complejo y orden no negativo", Sandia National Laboratory Report, Sand85-1018, mayo de 1985.

De Amos, "Un paquete portátil para funciones de Bessel de un argumento complejo y orden no negativo", trans. Matemáticas. Software, 1986.

Answer 2

La relación de recurrencia hacia adelante que está utilizando no es estable. Para ver por qué, considere que los valores de Besselj (N, X) se vuelven cada vez más pequeños por aproximadamente un factor 1/2n. Puedes ver esto mirando el primer término de la serie Taylor para J.

Entonces, lo que estás haciendo es restar un número grande de un múltiplo de un número algo más pequeño para obtener un número aún menor. Numéricamente, eso no va a funcionar bien.

Míralo de esta manera. Sabemos que el resultado es del orden de 10^-25. Comienza con números que son del orden de 1. Por lo tanto, para obtener incluso un dígito preciso de esto, tenemos que conocer los dos primeros números con al menos 25 dígitos de precisión. Claramente no lo hacemos, y la recurrencia en realidad diverge.

Usando la misma relación de recurrencia para ir hacia atrás, desde órdenes altas hasta órdenes bajas, es estable. Cuando comienza con los valores correctos para J (20,1) y J (19,1), también puede calcular todos los pedidos a 0 con una precisión total. ¿Por qué funciona esto? Porque ahora los números se están haciendo más grandes en cada paso. Estás restando un número muy pequeño de un múltiplo exacto de un número mayor para obtener un número aún mayor.

Answer 3

Simplemente puede modificar el código a continuación que es para la función esférica de Bessel. Está bien probado y funciona para todos los argumentos y rango de pedidos. Lo siento, es en 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;
        }
    }