Berechnen Sie die Bessel -Funktion in MATLAB unter Verwendung von JM+1 = 2MJ (M) -J (M -1) Formel

Question

Ich habe versucht, die Bessel -Funktion mit dieser Formel zu implementieren. Dies ist der 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;

Aber wenn ich die Bessel -Funktion von Matlab benutze, um sie mit diesem zu vergleichen, bekomme ich zu hoch unterschiedliche Werte. Wenn ich beispielsweise Bessel (20) eingeben kann, gibt es mir 3.1689E+005 als Ergebnis, wenn ich stattdessen Bessel (20,1) eingeben, gibt es mir 3.8735E-025, ein völlig anderes Ergebnis.

Answer 1

Solche Wiederholungsbeziehungen sind in der Mathematik nett, aber numerisch instabil, wenn Algorithmen mit begrenzten Präzisionsdarstellungen von Gleitkomma-Zahlen implementiert werden.

Betrachten Sie den folgenden Vergleich:

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

Sie können also sehen, wie sich die berechneten Werte nach 9 erheblich unterscheiden.

Laut Matlab:

Besselj verwendet eine Mex -Schnittstelle zu einer FORTRAN -Bibliothek von De Amos.

und gibt Folgendes als Referenzen für ihre Implementierung an:

De Amos, "Ein Subroutine-Paket für Bessel-Funktionen eines komplexen Arguments und einer nichtnegativen Ordnung", Sandia National Laboratory Report, Sand85-1018, Mai 1985.

De Amos, "ein tragbares Paket für Bessel -Funktionen eines komplexen Arguments und einer nichtnegativen Ordnung", Trans. Mathematik. Software, 1986.

Answer 2

Die von Ihnen verwendete Forward Recurquenz -Beziehung ist nicht stabil. Betrachten Sie, dass die Werte von besselj (n, x) um einen Faktor immer kleiner werden 1/2n. Sie können dies sehen, indem Sie sich die erste Amtszeit der Taylor -Serie für J. ansehen

Sie subtrahieren also eine große Zahl von einem Vielfachen einer etwas geringerer Anzahl, um eine noch kleinere Zahl zu erhalten. Numerisch wird das nicht gut funktionieren.

Schau es dir so an. Wir wissen, dass das Ergebnis der Reihenfolge von ist 10^-25. Sie beginnen mit Zahlen, die in der Reihenfolge von 1 sind. Um nur eine genaue Ziffer davon herauszuholen, müssen wir die ersten beiden Zahlen mit mindestens 25 Ziffern Präzision kennen. Wir tun es eindeutig, und das Wiederauftreten weicht tatsächlich ab.

Verwenden der gleichen Rezidivbeziehung, um rückwärts zu gehen, von hohen Bestellungen bis hin zu niedrigen Bestellungen, ist stabil. Wenn Sie mit den richtigen Werten für J (20,1) und J (19,1) beginnen, können Sie alle Bestellungen auch mit vollem Genauigkeit auf 0 berechnen. Warum funktioniert das? Denn jetzt werden die Zahlen in jedem Schritt größer. Sie subtrahieren eine sehr kleine Zahl von einem exakten Vielfachen einer größeren Zahl, um eine noch größere Zahl zu erhalten.

Answer 3

Sie können einfach den folgenden Code für die sphärische Bessel -Funktion ändern. Es ist gut getestet und funktioniert für alle Argumente und Auftragsbereiche. Es tut mir leid, dass es in C# ist

     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;
        }
    }