贝塞尔函数的自然对数，溢出

Question

我正在尝试在 MATLAB 中计算第二类修正贝塞尔函数的对数，即类似的东西：

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

但是，我收到以下错误：

Error using mupadmex Error in MuPAD command: DOUBLE cannot convert the input expression into a double array. If the input expression contains a symbolic variable, use the VPA function instead.

Error in sym/double (line 514) Xstr = mupadmex('symobj::double', S.s, 0)`;

如何避免这个错误？

Answer 1

你试过积分表示吗？

Log[Integrate[Cosh[Nu t]/E^(Z Cosh[t]), {t, 0, Infinity}]]

Answer 2

正如 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 ) + 日志 (1.02) = 4.7318e+03。所以对于 nu >= 750 和 z=1，至少 5 位有效数字是正确的！如果这对你来说足够好，这将比符号数学快得多。

Answer 3

你做错了几件事。将 double 用于 besselk 的两个参数并将输出转换为符号是没有意义的。您还应该避免使用旧的基于字符串的输入 sym。相反，您想以符号方式计算 besselk（return 大约 1.02×10²⁰⁵⁵，远大于 realmax），取对数符号化结果，然后转换回双精度。

以下就足够了——当一个或多个输入参数是符号时，将使用 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 和 Z 值，因为大订单 (nu) 并不常见。

贝塞尔函数的自然对数，溢出

Natural Logarithm of Bessel Function, Overflow

algorithm

math

matlab

bessel-functions