5

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

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 というまったく異なる結果が得られます。

4

3 に答える 3

3

recurrence_bessel

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

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

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、「複雑な引数と非負の順序のベッセル関数のサブルーチン パッケージ」、サンディア国立研究所レポート、SAND85-1018、1985 年 5 月。

DE Amos、「複雑な引数と非負の順序の Bessel 関数のポータブル パッケージ」、Trans. 算数。ソフトウェア、1986年。

于 2011-10-27T05:46:47.147 に答える
2

使用している前方再帰関係は安定していません。その理由を理解するには、BesselJ(n,x) の値が約 1 倍ずつ小さくなることを考慮してください1/2n。これは、J のテイラー級数の第 1 項を見ればわかります。

したがって、あなたがしているのは、やや小さい数の倍数から大きな数を引いて、さらに小さな数を取得することです。数値的には、うまくいきません。

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

同じ再帰関係を使用して、上位から下位へと逆方向に移動すると、安定します。J(20,1) と J(19,1) の正しい値から始めると、0 までのすべての次数を完全な精度で計算することもできます。なぜこれが機能するのですか？今では、各ステップで数値が大きくなっているからです。大きな数の正確な倍数から非常に小さな数を引いて、さらに大きな数を取得しています。

于 2011-10-27T07:38:31.217 に答える
0

球面ベッセル関数用の以下のコードを変更するだけです。これは十分にテストされており、すべての引数と次数範囲で機能します。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;
        }
    }
于 2013-03-12T13:20:49.590 に答える