How can I calculate the second-order derivative of a vector using finite differences if the interval is non-constant?

Question

Say I have vectors x and y and want to calculate the second derivative of y with respect to x using finite differences.

I'd do

x <- rnorm(2000)
y <- x^2
y = y[order(x)]
x = sort(x)

dydx = diff(y) / diff(x)
d2ydx2 = c(NA, NA, diff(dydx) / diff(x[-1]))
plot(x, d2ydx2)

As you can see, there are a few points which are wildly inaccurate. I believe the problem arises because values in dydx do not exactly correspond to those of x[-1] leading a second differentiation to have inaccurate results. Since the step in x is non-constant, the second-order differentiation is not straight forward. How can I do this?

Answer 1

Each time you are taking the numerical approximation derivative, you are losing one value in the vector and shifting the output value over one spot. You are correct, the error is due to the uneven spacing in the x values (incorrect divisor in dydx & d2ydx2 calculations).

In order to correct, calculate a new set of x values corresponding to the mid point between the adjacent x values at each derivative. This is the value where the slope is calculated.
Thus y'1 = f'((x1+x2)/2).

This method is not perfect but the resulting error is much smaller.

#create the input
x <- sort(rnorm(2000))
y <- x**2

#calculate the first deriative and the new mean x value
xprime <- x[-1] - diff(x)/2
dydx <- diff(y)/diff(x)

#calculate the 2nd deriative and the new mean x value
xpprime <- xprime[-1] - diff(xprime)/2
d2ydx2 <- diff(dydx)/diff(xprime)

plot(xpprime, d2ydx2)

Answer 2

Another way is using splinefun, which returns a function from which you can calculate cubic spline derivatives. Of course, given your example function y= x^2 the second derivatives will be always 2

x <- rnorm(2000)
y <- x^2
y = y[order(x)]
x = sort(x)

fun = splinefun(x,y)

plot(x,fun(x,deriv=2))