La fonction de Calcul en utilisant MATLAB Jm + 1 = 2 mJ (m) -j (m-1) de formule

Question

J'ai essayé de mettre en œuvre Bessel fonction à l'aide de cette formule, c'est le code:

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;

Mais si j'utilise la fonction de Bessel de Matlab pour le comparer avec celui-ci, je reçois trop élevé des valeurs différentes. Par exemple, si je tape Bessel (20), il me donne 3.1689e + 005 résultat, si au lieu de type I Bessel (20,1), il me donne 3.8735e-025, un résultat totalement différent.

Answer 1

ces relations de récurrence sont bien en mathématiques, mais numériquement instable lors de l'implémentation des algorithmes utilisant des représentations de précision limitées des nombres à virgule flottante.

Considérons la comparaison suivante:

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

Vous pouvez voir comment les valeurs calculées commencent à différer de manière significative après 9.

Selon Matlab:

BESSELJ utilise une interface MEX à une bibliothèque Fortran par D. E. Amos.

et donne les références suivantes pour leur mise en œuvre:

D. E. Amos, « Un ensemble de sous-programme pour les fonctions de Bessel d'un complexe l'argument et l'ordre non négatif », Sandia National Laboratory Rapport, SAND85-1018, mai 1985.

D. E. Amos, « Emballage portatif pour les fonctions de Bessel d'un complexe l'argument et l'ordre non négatif », Trans. Math. Software, 1986.

Answer 2

La relation de récurrence avant que vous utilisez est pas stable. Pour comprendre pourquoi, considèrent que les valeurs de BESSELJ (n, x) deviennent plus petits et plus petits d'environ un 1/2n facteur. Vous pouvez le voir en regardant le premier terme de la série de Taylor J.

Alors, ce que vous faites est la soustraction d'un grand nombre d'un multiple d'un nombre un peu plus petit pour obtenir un nombre encore plus petit. Numériquement, qui ne va pas bien fonctionner.

Regardez cette façon. Nous savons que le résultat est de l'ordre de 10^-25. Vous commencez avec des chiffres qui sont de l'ordre de 1. Ainsi, afin d'obtenir encore un chiffre précis sur ce, nous devons connaître les deux premiers chiffres avec au moins 25 chiffres de précision. Nous faisons clairement pas, et la récurrence fait diverge.

En utilisant la même relation de récurrence pour revenir en arrière, de commandes à haut bas commandes, stable. Lorsque vous commencez avec des valeurs correctes pour J (20,1) et J () 19,1, vous pouvez calculer toutes les commandes à 0 avec une précision totale aussi bien. Pourquoi ce travail? Parce que maintenant les chiffres deviennent plus grandes à chaque étape. Vous retranchant un très petit nombre d'un multiple exact d'un plus grand nombre pour obtenir un plus grand nombre.

Answer 3

Vous pouvez modifier simplement le code ci-dessous qui est pour la fonction de la rotule. Il est bien testé et fonctionne pour tous les arguments et la gamme de commande. Je suis désolé, il est 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;
        }
    }