如何用R计算第二类修正贝塞尔函数对数的一阶与二阶导数？

一阶导数计算 #

先明确一阶导数的核心表达式：
$$\frac{\partial}{\partial x} \log(K_\nu(x)) = \frac{K'\nu(x)}{K\nu(x)}$$

替代现有方案的优化选项 #

• 单包完成计算（无需fOptions）
  Bessel包的besselK支持直接返回导数，设置deriv=1即可一次性拿到函数值和一阶导数，不用分开调用两个函数：

  library(Bessel)
  # 输入x和ν，同时获取K_v(x)和它的一阶导数
  calc_result <- besselK(x = 2, nu = 1.5, deriv = 1)
  # 直接计算对数的一阶导数
  log_k_1st_deriv <- calc_result$deriv / calc_result$value
  

• 手动递推计算
  第二类修正贝塞尔函数的导数有现成递推公式：
  $$K'\nu(x) = -\frac{\nu}{x}K\nu(x) - K_{\nu-1}(x)$$
  用基础的besselK就能实现，不需要依赖导数专用函数：

  x <- 2
  nu <- 1.5
  k_nu <- besselK(x, nu)
  k_nu_minus1 <- besselK(x, nu - 1)
  log_k_1st_deriv <- (-nu/x * k_nu - k_nu_minus1) / k_nu
  

二阶导数计算 #

对数的二阶导数可以通过一阶导数进一步推导得到：
$$\frac{\partial^2}{\partial x^2} \log(K_\nu(x)) = \frac{K''\nu(x)K\nu(x) - [K'_\nu(x)]^{2}{[K_\nu(x)]}2}$$

解析递推法（精度更高） #

利用第二类修正贝塞尔函数的二阶导数递推公式：
$$K''\nu(x) = \frac{\nu^2 - x^2}{x2}K\nu(x) + \frac{1}{x}\left(K_{\nu-1}(x) + K_{\nu+1}(x)\right)$$
对应的R代码实现：

x <- 2
nu <- 1.5
# 计算所需的贝塞尔函数值
k_nu <- besselK(x, nu)
k_nu_minus1 <- besselK(x, nu - 1)
k_nu_plus1 <- besselK(x, nu + 1)

# 计算K_v(x)的二阶导数
k_nu_2nd_deriv <- ((nu^2 - x^2)/x^2)*k_nu + (k_nu_minus1 + k_nu_plus1)/x
# 计算K_v(x)的一阶导数（用递推式）
k_nu_1st_deriv <- (-nu/x * k_nu - k_nu_minus1)

# 最终计算对数的二阶导数
log_k_2nd_deriv <- (k_nu_2nd_deriv * k_nu - k_nu_1st_deriv^2) / (k_nu^2)

数值微分法（快速验证） #

如果对精度要求不高，用numDeriv包的数值微分工具可以快速算出结果，不用推导公式：

library(numDeriv)
x <- 2
nu <- 1.5

# 定义对数贝塞尔函数
log_k <- function(x) log(besselK(x, nu))

# 一阶导数（数值法）
log_k_1st_deriv_num <- grad(log_k, x)
# 二阶导数（数值法）
log_k_2nd_deriv_num <- hessian(log_k, x)

内容的提问来源于stack exchange，提问作者S S