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

【问题标题】：Natural Logarithm of Bessel Function, Overflow贝塞尔函数的自然对数，溢出
【发布时间】：2015-12-05 17:25:57
【问题描述】：

我正在尝试在 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)))'));

但是，我收到以下错误：

使用 mupadmex 时出错 MuPAD 命令中的错误：DOUBLE 无法将输入表达式转换为双精度数组。如果输入表达式包含符号变量，请改用 VPA 函数。

sym/double 中的错误（第 514 行）Xstr = mupadmex('symobj::double', S.s, 0)`;

我怎样才能避免这个错误？

【问题讨论】：

是否溢出？或者结果实际上是无限的？你的值大于1.7977e+308，你需要玩这么大的数值吗？
旁注：它不会溢出，但会给出Inf。从技术上讲，数字会再次从负数开始溢出，但这只会发生在整数上，而不是浮点数。
@AnderBiguri 从技术上讲，溢出是数字再次从负数开始，但这只会发生在整数上，而不是浮点数。您所指的内容称为“integer wrap around”。溢出是“达到极限”的总称，applies to floating-point numbers as well.
DLMF 为section 10.40 中的大参数和section 10.41 中的大订单提供渐近扩展，这可能适合您的用例。
this question on Mathoverflow 的答案也可能有用。

标签： algorithm matlab math bessel-functions

【解决方案1】：

你做错了一些事情。将double 用于besselk 的两个参数并将输出转换为符号是没有意义的。您还应该避免对sym 进行旧的基于字符串的输入。相反，您想以符号方式计算 besselk（这将返回大约 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) 并不常见。

【讨论】：

【解决方案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 个有效数字是正确的！如果这对您来说足够好，这将比符号数学快得多。

【讨论】：

谢谢。我会尝试使用这个近似值，但是我不确定精度是否足够......
@JamesGreen：我刚刚尝试过 nu=750 和 z=1。查看我的编辑。

【解决方案3】：

你试过积分表示吗？

对数[积分[Cosh[Nut]/E^(Z Cosh[t]), {t, 0, Infinity}]]

【讨论】：

感谢您的建议。集成很慢，我必须做很多迭代......