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

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

但是,我收到以下错误:

在MuPAD命令中使用mupadmex时出错:DOUBLE无法将输入表达式转换为double数组.如果输入表达式包含符号变量,请改用VPA函数.

sym/double中的错误(第514行)Xstr = mupadmex('symobj :: double',Ss,0)`;

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

Answer 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)不常见.

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) + log(1.02 ) = 4.7318e+03。因此，对于 nu >= 750 且 z=1，它至少有 5 位有效数字！如果这对你来说足够好，这将比符号数学快得多。