How would you explain this failure of numerical twice-differentiation with the numDeriv package ?
Below is a function Y and its first derivative 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)
}
And below is a successful illustration of the grad function of the numDeriv package:
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
Twice-differentiation below works fine too:
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
But when replacing Y with a sharp approximate, twice-differientation fails:
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
However when replacing dYnum with an approximate, it works again:
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
You forgot the inner derivative in the last term. The a in the denominator cancels with the inner derivative. So it should be
-a*A*sin(a*x) + a*B*cos(a*x) - cos(a*x)/2/a + x/2*sin(a*x)
As for the error, I would suspect that the culprit is the approximation. Depending on the type of the approximation, the error between the sample points in the function values can be arbitrarily large, which only magnifies when computing derivatives. You might get better results by using far less sample points, 6 instead of 99.
RTFM: Per R-manual approxfun only provides piecewise constant or linear approximations. Both are not differentiable. This means that the derivative of the piecewise linear approximation is piecewise constant with jumps and the second derivative is piecewise zero with delta peaks at the sample points. Re-interpolating the first derivative via a piecewise linear function will, while not the most mathematical procedure, lead back to reasonable results.
Use R-manual splinefun to get an interpolating function that actually is twice differentiable.