JM+1 = 2MJ（M）-J（M -1）式を使用してMATLABのベッセル関数を計算する

Question

私はその式を使用してベッセル関数を実装しようとしました。これはコードです。

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を与えます。

Answer 1

このような再発関係は、数学では優れていますが、浮動小数点数の限られた精密表現を使用してアルゴリズムを実装する場合は数値的に不安定です。

次の比較を検討してください。

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

したがって、9以降、計算された値がどのように大きく異なるかを見ることができます。

Matlabによると：

Besseljは、De AmosのFortranライブラリにMEXインターフェイスを使用しています。

以下を実装の参照として示します。

De Amos、「複雑な議論と非陰性秩序のBessel機能のためのサブルーチンパッケージ」、Sandia National Laboratory Report、Sand85-1018、1985年5月。

De Amos、「複雑な議論と非陰性秩序のBessel機能のためのポータブルパッケージ」、Trans。算数。ソフトウェア、1986年。

Answer 2

使用している前方再発関係は安定していません。理由を確認するには、Besselj（n、x）の値がほぼ係数によって小さくて小さくなることを考慮してください 1/2n. 。 Jのテイラーシリーズの最初の期間を見ることで、これを見ることができます。

したがって、あなたがしていることは、やや小さい数の倍数から多数を差し引いて、さらに少ない数を得ることです。数値的には、それはうまく機能しません。

このように見てください。結果は次のことを知っています 10^-25. 。 1のオーダーである数字から始めます。したがって、これから1桁の数字を1つでも取得するには、少なくとも25桁の正確な最初の2つの数値を知る必要があります。私たちは明らかにそうではなく、再発は実際に分岐します。

同じ再発関係を使用して、高次から低いオーダーまで、後方に進む、 は 安定。 J（20,1）とJ（19,1）の正しい値から始めると、すべての注文を完全な精度で0に計算できます。なぜこれが機能するのですか？これで、各ステップで数値が大きくなっているためです。さらに多くの数字の正確な倍数から非常に少数を差し引いて、さらに多くの数を取得します。

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