8

MATLAB で 2 番目のタイプの修正ベッセル関数の対数を計算しようとしています。つまり、次のようなものです。

log(besselk(nu, Z))

例えば

nu = 750;
Z = 1;

の値がlog(besselk(nu, Z))無限大になるので問題besselk(nu, Z)があります。しかし、log(besselk(nu, Z)) 確かに小さいはずです。

私は次のようなものを書こうとしています

f = double(sym('ln(besselk(double(nu), double(Z)))'));

ただし、次のエラーが表示されます。

mupadmex の使用エラー MuPAD コマンドのエラー: DOUBLE は入力式を double 配列に変換できません。入力式にシンボリック変数が含まれている場合は、代わりに VPA 関数を使用してください。

sym/double のエラー (514 行目) Xstr = mupadmex('symobj::double', Ss, 0)`;

このエラーを回避するにはどうすればよいですか?

4

3 に答える 3

5

あなたはいくつかのことを間違っています。double2 つの引数を to に使用しbesselk、出力をシンボリックに変換するのは意味がありません。への古い文字列ベースの入力も避ける必要がありますsym。代わりに、besselkシンボリックに評価し (約 1.02×10 2055を返し、 よりもはるかに大きいrealmax)、結果のログをシンボリックに取得してから、倍精度に変換し直します。

以下で十分です。1 つ以上の入力引数がシンボリックである場合、シンボリック バージョンのbesselkが使用されます。

f = double(log(besselk(sym(750), sym(1))))

または古い文字列形式:

f = double(sym('log(besselk(750, 1))'))

パラメータをシンボリックに保ち、後で評価する場合:

syms nu Z;
f = log(besselk(nu, Z))
double(subs(f, {nu, Z}, {750, 1}))

大量の注文 ( ) はあまり一般的ではないため、数学でnuとの値を反転していないことを確認してください。Znu

于 2015-09-09T19:35:45.663 に答える
4

njuffa が指摘したように、DLMF は大きな nu に対して K_nu(z) の漸近展開を与えます。10.41.2 から、真の正の引数 z について次のことがわかります。

besselk(nu,z) ~ sqrt(pi/(2nu)) (e z/(2nu))^-nu

これは、いくつかの単純化の後に与えられます

log( besselk(nu,z) ) ~ 1/2*log(pi) + (nu-1/2)*log(2nu) - nu(1 + log(z))

つまり、O(nu log(nu)) です。nu > 750 の場合、直接計算が失敗することは驚くことではありません。

この近似がどれほど正確かはわかりません。おそらく、besselk が無限大よりも小さい値を比較して、目的に合っているかどうかを確認できますか?

編集: nu=750 と z=1 で試してみました: 上記の概算では 4.7318e+03 が得られますが、horchler の結果では log(1.02*10^2055) = 2055*log(10) + log(1.02 ) = 4.7318e+03. したがって、nu >= 750 および z=1 の場合、少なくとも有効数字 5 桁までは正しいです! これで十分な場合は、記号演算よりもはるかに高速になります。

于 2015-09-09T18:54:58.863 に答える