Calcolare la funzione di Bessel in MATLAB utilizzando Jm + 1 = 2mj (m) -j formula (m-1)

Question

Ho cercato di realizzare bessel funzione utilizzando tale formula, questo è il codice:

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;

Ma se io uso la funzione di Bessel di MATLAB per confrontarlo con questo, ho troppo alti valori diversi. Per esempio, se si digita Bessel (20) Mi dà 3.1689e + 005, come conseguenza, se invece di tipo I di Bessel (20,1) che mi dà 3.8735e-025, un risultato completamente diverso.

Answer 1

tali relazioni di ricorrenza sono belle in matematica ma numericamente instabile nell'attuazione algoritmi mediante rappresentazioni precisione limitato di numeri in virgola mobile.

Si consideri il seguente paragone:

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

Quindi, si può vedere come i valori calcolati iniziano a differire in modo significativo dopo il 9.

Secondo MATLAB:

BESSEL.J utilizza un'interfaccia MEX a una libreria Fortran di D. E. Amos.

e dà il seguente come riferimento per la loro attuazione:

D. E. Amos, "pacchetto Una subroutine per le funzioni di Bessel di un complesso argomento e l'ordine non negativo", relazione Sandia National Laboratory, SAND85-1018, maggio, 1985.

D. E. Amos, "Un pacchetto portatile per funzioni di Bessel di un complesso argomento e l'ordine non negativo", Trans. Math. Software, 1986.

Answer 2

La relazione di ricorrenza in avanti che si sta utilizzando non è stabile. Per capire perché, ritengono che i valori di BESSEL.J (n, x) diventano sempre più piccole di circa un fattore 1/2n. Si può vedere questo guardando il primo termine della serie di Taylor per J.

Quindi, quello che stai facendo è sottrarre un gran numero da un multiplo di un numero leggermente inferiore per ottenere un numero ancora più piccolo. Numericamente, che non sta andando a lavorare bene.

guardare in questo modo. Sappiamo che il risultato è dell'ordine di 10^-25. Si inizia con i numeri che sono dell'ordine di 1. Quindi, al fine di ottenere anche una cifra accurata fuori da questo, dobbiamo conoscere i primi due numeri con almeno 25 cifre di precisione. Noi chiaramente non, e la ricorrenza in realtà diverge.

Utilizzando la stessa relazione di ricorrenza tornare indietro, da alti ordini a basse ordini, è stabile. Quando si avvia con i valori corretti per J (20,1) e J (19,1), è possibile calcolare tutti gli ordini fino a 0, con tutta esattezza pure. Perché questo lavoro? Perché ora i numeri sono sempre più grandi in ogni passaggio. Si sta sottraendo un numero molto piccolo di un multiplo esatto di un numero maggiore di ottenere un numero ancora maggiore.

Answer 3

Si può semplicemente modificare il codice qui sotto, che è per la funzione sferica di Bessel. E 'ben testato e funziona per tutti gli argomenti e la gamma ordine. Mi dispiace che è in 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;
        }
    }