Why does a kernel density estimate of a normally distributed vector have a non-smooth second derivative? I had some normally distributed data:
mu <- 3
sigma <- 5
x <- rnorm(1e5, mu, sigma)

I took a kernel density estimate with a fairly high bandwidth:
kernel_density_of_x <- density(x, bw = 5)

Then I differentiated it:
differentiate <- function(x, y)
{
  diffOfX <- diff(x)
  data.frame(
    x      = x[-length(x)] + (diffOfX / 2), 
    dyByDx = diff(y) / diffOfX
  )
}

first_derivative <- with(kernel_density_of_x, differentiate(x, y))

This looked just as expected:
library(ggplot2)
(p1 <- ggplot(first_derivative, aes(x, dyByDx)) + geom_line())


When I differentiated again, I expected another smooth curve, but I saw an odd cyclical effect.
second_derivative <- with(first_derivative, differentiate(x, dyByDx))
(p2 <- p1 %+% second_derivative + ylab("d2yByDx2"))


I tried a few different options for the kernel argument, but the noisiness persisted.
Dropping the bandwidth down to, for example, 0.5 gave a lower frequency noise that dominated the plot (making it nonsense).
Dropping the number of sampling points down from n = 512 to n = 32 stopped the issue, but that causes other problems.
Why does this effect occur? Is it an artifact of the density function, or have I done something silly?

We can redraw the plot using the probability density function of the normal distribution that x was generated from to see the shape I expected:
xx <- seq.int(-20, 20, 0.1)
pdf_of_xx <- dnorm(xx, mu, sigma)
first_derivative_of_xx <- differentiate(xx, pdf_of_xx)
second_derivative_of_xx<- with(first_derivative_of_xx, differentiate(x, dyByDx))
ggplot(second_derivative_of_xx, aes(x, dyByDx)) + geom_line()


 A: As whuber commented, what you're seeing is due to the approximation (via a fast fourier transform) that density uses.
If you calculate the kernel density estimate by brute force and compare it to the estimate that density gives, you'll see this sort of cyclic pattern in the differences.
The 2nd differences in the result from density result in an exaggeration of the effect. The second differences from the brute-force kernel density estimate don't show that pattern, but the calculation is like 3 orders of magnitude slower.
I put some code in this gist. My brute-force kernel density estimate is the following:
mydensity <-
function(dat, x, bw=5) # dat=the data; x=points to calculate density estimate
{
  y <- vapply(x, function(a) mean(dnorm(a, dat, bw)), 1)
  data.frame(x=x, y=y)
}

Here's a picture of the second differences, blue from density and red from brute-force approach.

A: In your codes in the gist you used a hard-entered bandwidth (bw=5) which is very different from a optimal one returned for example as:
library(ks)
h <- hpi(x)

which is actually about h=0.85.
It is interesting to execute your codes where bw is set to h. Then the final figure is as below. Here the differences are not so drastic.  
