r - 使用 numDeriv R 包进行二次微分失败

Question

您如何解释使用 numDeriv 包进行数值二次微分的失败？

下面是函数Y及其一阶导数dY。

Y <- function(x, A=1, B=1, a=1){
  A*cos(a*x) + B*sin(a*x) - x/2/a*cos(a*x)
}
dY <- function(x, A=1, B=1, a=1){
  -a*A*sin(a*x) + a*B*cos(a*x) - cos(a*x)/2/a + x/2/a*sin(a*x)
}

下面是 numDeriv 包的 grad 函数的成功说明:

library(numDeriv)
x <- c(0.2,0.4,0.6)
grad(Y, x)
## [1] 0.31123 0.14900 0.01742
dY(x)
## [1] 0.31123 0.14900 0.01742

下面的两次微分也可以很好地工作:

dYnum <- function(x) grad(Y, x)
d2Ynum <- function(x) grad(dYnum, x)
d2Ynum(x)
## [1] -0.8820 -0.7368 -0.5777
grad(dY, x)
## [1] -0.8821 -0.7368 -0.5777

但是当用锐近似替换 Y 时，二次微分失败:

xx <- seq(0.01, 0.99, by=0.01)
Yapprox <- approxfun(xx, Y(xx))
dYnum <- function(x) grad(Yapprox, x)
dYnum(x) # first-order derivative is OK
## [1] 0.31124 0.14901 0.01743
d2Ynum <- function(x) grad(dYnum, x)
d2Ynum(x) # not the good result
## [1] -2698 -1127  -589
grad(dY, x)  # this is the good result
## [1] -0.8821 -0.7368 -0.5777

但是，当用近似值替换 dYnum 时，它会再次起作用:

xx <- seq(0.1, 0.9, by=0.01)
dYnumapprox <- approxfun(xx, dYnum(xx)) # approxfun(xx, grad(Yapprox, xx))
d2Ynum_ <- function(x) grad(dYnumapprox, x)
d2Ynum_(x) 
## [1] -0.8820 -0.7368 -0.5777

Answer 1

你忘记了上一项的内导数。分母中的 a 与内导数相抵消。所以应该是这样

  -a*A*sin(a*x) + a*B*cos(a*x) - cos(a*x)/2/a + x/2*sin(a*x)

---

至于错误，我怀疑罪魁祸首是近似值。根据近似的类型，函数值中样本点之间的误差可以是任意大，只有在计算导数时才会放大。使用更少的样本点(6 个而不是 99 个)可能会获得更好的结果。

---

RTFM:根据 R-manual approxfun仅提供分段常数或线性近似值。两者不可微分。这意味着分段线性近似的导数是带有跳跃的分段常数，并且二阶导数是分段零，并且在样本点处有 delta 峰值。通过分段线性函数重新插值一阶导数虽然不是最数学的过程，但会得到合理的结果。

使用R-manual splinefun获得实际上两次可微的插值函数。