By finite differences, the first order derivative of y for each mean value of x over your array is given by :

dy=np.diff(y,1)
dx=np.diff(x,1)
yfirst=dy/dx

And the corresponding values of x are :

xfirst=0.5*(x[:-1]+x[1:])

For the second order, do the same process again :

dyfirst=np.diff(yfirst,1)
dxfirst=np.diff(xfirst,1)
ysecond=dyfirst/dxfirst

xsecond=0.5*(xfirst[:-1]+xfirst[1:])

I think I have figured out the answer. The results have to do with the difference between the integral of a continuous Gaussian function, versus the sum of a discretized Gaussian.

In the ideal case, continuous and discrete sums of Gaussians already differ

As one illustrative example, consider the integral of a Gaussian versus the discrete summation of a Gaussian with unit steps is $$ \int_{-\infty}^\infty e^{-x^2}\!\text{d}x=\sqrt{\pi}\approx 1.77245,\qquad \sum_{x\in\mathbb{Z}} e^{-x^2} = \theta_3(0,e^{-1})\approx 1.77264.$$ Apparently, the two expressions differ in the details (after 4 decimals). Such things matter when we compare a continuous convolution () against a discrete convolution () using Gaussian kernels.

Symmetric versus asymmetric functions

A Gaussian is a symmetric function, i.e., . The first derivative of a Gaussian is antisymmetric, i.e., . The second derivative of a Gaussian is symmetric again, and so on.

Implication for convolutions

We consider truncated (to terms, although we include in this description) convolutions with Gaussians or their -th derivatives () with the following computation: Knowing that Gaussian functions are either anti-symmetric or symmetric, that means we can 'simplify' the summation to only consider non-negative integer values for ,

We now note that for an odd (i.e., first-order, third-order, ...) derivative of a Gaussian, the antisymmetry also implies that . This is not true for an even derivative of a Gaussian, where .

Specific example for constant functions,

We can thus symbolically write down what happens when convolving a Gaussian derivative with a constant function . For a first derivative, with , we compute that is, it evaluates to 0, which we expect from differentiating a constant function.

Conversely, for a second derivative, , it results in This expression only equals the expected result of 0, if the truncated and discrete kernel sums to 0, which is alternatively expressed as the expression above. This strict requirement is not met for discrete sums of a Gaussian derivative, even without truncation (i.e., when letting the sum is still not 0)!

Explaining the observed behavior

So, I can now answer the question posed above (which was 'why do two successive 1st-order derivatives give the correct result, but a single 2nd-order derivative not?'): after applying the first derivative operation, only zeroes remain (per the result for described above), which leads to field of zeroes for the second convolution operation. Conversely, an immediate application of simply returns the sum of the kernel multiplied with the local value. That second definition should be 0 in the ideal and continuous case -- but because we're working with a discrete application of convolutions it ceases to hold. Even if we use a very large truncation (e.g., specify truncate=1000 when calling the Gaussian filter in Python), it will not converge to 0. Instead, it will converge to (the term comes from the Gaussian which is applied in the horizontal direction), which we can confirm with our code:

import numpy as np
import matplotlib.pyplot as plt
A = np.zeros((30,30))
A[:,15] = 1
from scipy.ndimage import gaussian_filter
B = gaussian_filter(A, sigma=1, order=[2, 0], mode='reflect', truncate=10000)
plt.imshow(B)
plt.colorbar()
plt.show()
print(B[0,15]) # >> -8.426948379372423e-08

So, to conclude, the scipy implementation is correct, and it does what it is supposed to do, from a discrete point of view. However, from a continuous point of view, the results are slightly unexpected and seem wrong. I'll ponder what the most appropriate fix is. Perhaps, we must simply ensure that the summation convention holds, by slightly modifying the kernel , after which the discrete kernel will mimic the continuous kernel, at least for constant functions!

Maxim Umansky’s answer describes the storage convention of the FFT frequency components in detail, but doesn’t necessarily explain why the original code didn’t work. There are three main problems in the code:

1. x = linspace(0,2*pi,N): By constructing your spatial domain like this, your x values will range from $0$ to $2\pi$, inclusive! This is a problem because your function y = sin(2*x)+cos(5*x) is not exactly periodic on this domain ($0$ and $2\pi$ correspond to the same point, but they’re included twice). This causes spectral leakage and thus a small deviation in the result. You can avoid this by using x = linspace(0,2*pi,N, endpoint=False) (or x = 2*pi*arange(0,N,1)/N, as Maxim Umansky did; this is what he is referring to with “open-periodic domain”).
2. k = fftshift(k): As Maxim Umansky explained, your k values need to be in a specific order to match the FFT convention. fftshift sorts the values (from small/negative to large/positive), which can be useful e. g. for plotting, but is incorrect for computations.
3. dydx1 = ifft(-k*1j*fft(y)).real: scipy defines the FFT as y(j) = (x * exp(-2*pi*sqrt(-1)*j*np.arange(n)/n)).sum(), i. e. with a factor of $2\pi$ in the exponential, so you need to include this factor when deriving the formula for the derivative. Also, for scipy’s FFT convention, the k values shouldn’t get a minus sign.

So, with these three changes, the original code can be corrected as follows:

from scipy.fftpack import fft, ifft, dct, idct, dst, idst, fftshift, fftfreq
from numpy import linspace, zeros, array, pi, sin, cos, exp
import matplotlib.pyplot as plt

N = 100
x = linspace(0,2*pi,N, endpoint=False) # (1.)

dx = x[1]-x[0]
y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)

k = fftfreq(N,dx)
# (2.)

dydx1 = ifft(2*pi*k*1j*fft(y)).real # (3.)

plt.plot(x,dydx,'b',label='Exact value')
plt.plot(x,dydx1,'r',label='Derivative by FFT')
plt.legend()
plt.show()