1

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?

Anthony Tan
  • 13
  • 2
  • *"As you can see"* ... did you forget a picture? If not, consider redoing that with a known random seed so that we can reproduce your findings. – r2evans Jan 05 '23 at 02:01
  • what results are you expecting? why order y and sort x? – megmac Jan 05 '23 at 03:10

2 Answers2

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)

enter image description here

Dave2e
  • 22,192
  • 18
  • 42
  • 50
1

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

Ric
  • 5,362
  • 1
  • 10
  • 23