我支持jrennie的第一句话--这要看情况而定。numpy.gradient函数要求数据均匀分布（尽管多维时允许每个方向上的距离不同）。如果您的数据不符合这一点，那么numpy.gradient就没有多大用处。实验数据可能（好的，会有）噪声，此外，不一定是所有均匀的间隔。在这种情况下，它可能是更好的使用scipy之一。插值样条函数（或对象）。这些函数可以采用不均匀间隔的数据，允许平滑，并且可以返回最大为k-1的导数，其中k是所请求的样条拟合的阶数。k的默认值为3，因此二阶导数就足够了。例如：

spl = scipy.interpolate.splrep(x,y,k=3) # no smoothing, 3rd order spline
ddy = scipy.interpolate.splev(x,spl,der=2) # use those knots to get second derivative

面向对象的样条函数，如scipy.interpolate.UnivariateSpline，有导数的方法。注意导数方法在Scipy 0.13中实现，在0.12中没有。
请注意，正如@JosephCottham在2018年的评论中所指出的，这个答案（至少对Numpy 1.08有效），从（至少）Numpy 1.14开始就不再适用了。

数值梯度计算没有一个通用的正确答案，在计算样本数据的梯度之前，您必须对生成该数据的底层函数做一些假设。从技术上讲，您可以使用np.diff进行梯度计算。使用np.gradient是一种合理的方法。我看不出您所做的有什么根本性的错误---它'这是一维函数的二阶导数的一种特殊近似。

由于梯度函数考虑了向左和向右的一个数据点，因此在多次应用时，这种情况会继续/扩散。
另一方面，二阶导数可通过以下公式计算：

d^2 f(x[i]) / dx^2 = (f(x[i-1]) - 2*f(x[i]) + f(x[i+1])) / h^2

比较here。这具有仅考虑两个相邻像素的优点。

在图中，双np.gradient方法（左）和上述公式（右）（由np.diff实现）进行了比较。由于f（x）在零点只有一个扭结，二阶导数（绿色）应该只有一个峰值。由于双梯度解决方案在每个方向上考虑了2个相邻点，这导致在+/- 1处有有限的二阶导数值。

• 但在某些情况下，您可能希望使用双梯度解决方案，因为它对噪声更稳健。*

我不知道为什么会有np.gradient和np.diff，但原因可能是，np.gradient的第二个参数定义了像素距离（对于每个维度），对于图像，它可以同时应用于两个维度gy, gx = np.gradient(a)。

代码

import numpy as np
import matplotlib.pyplot as plt

xs = np.arange(-5,6,1)
f = np.abs(xs)
f_x = np.gradient(f)
f_xx_bad = np.gradient(f_x)
f_xx_good = np.diff(f, 2)
test = f[:-2] - 2* f[1:-1] + f[2:]

# lets plot all this

fig, axs = plt.subplots(1, 2, figsize=(9, 3), sharey=True)

ax = axs[0]
ax.set_title('bad: double gradient')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs, f_xx_bad, marker='o', label='d^2 f(x) / dx^2')
ax.legend()

ax = axs[1]
ax.set_title('good: diff with n=2')
ax.plot(xs, f, marker='o', label='f(x)')
ax.plot(xs, f_x, marker='o', label='d f(x) / dx')
ax.plot(xs[1:-1], f_xx_good, marker='o', label='d^2 f(x) / dx^2')
ax.plot(xs[1:-1], test, marker='o', label='test', markersize=1)
ax.legend()

当我不断地以一种形式或另一种形式一次又一次地解决这个问题时，我决定写一个函数gradient_n，它给np.gradient增加了一个微分或功能。不是所有的np.gradient功能都被支持，比如多轴微分。
与np.gradient类似，gradient_n返回与输入相同形状的微分结果。此外，它还支持像素距离参数（d）。

import numpy as np

def gradient_n(arr, n, d=1, axis=0):
    """Differentiate np.ndarray n times.

    Similar to np.diff, but additional support of pixel distance d
    and padding of the result to the same shape as arr.

    If n is even: np.diff is applied and the result is zero-padded
    If n is odd: 
        np.diff is applied n-1 times and zero-padded.
        Then gradient is applied. This ensures the right output shape.
    """
    n2 = int((n // 2) * 2)
    diff = arr

    if n2 > 0:
        a0 = max(0, axis)
        a1 = max(0, arr.ndim-axis-1)
        diff = np.diff(arr, n2, axis=axis) / d**n2
        diff = np.pad(diff, tuple([(0,0)]*a0 + [(1,1)] +[(0,0)]*a1),
                    'constant', constant_values=0)

    if n > n2:
        assert n-n2 == 1, 'n={:f}, n2={:f}'.format(n, n2)
        diff = np.gradient(diff, d, axis=axis)

    return diff

def test_gradient_n():
    import matplotlib.pyplot as plt

    x = np.linspace(-4, 4, 17)
    y = np.linspace(-2, 2, 9)
    X, Y = np.meshgrid(x, y)
    arr = np.abs(X)
    arr_x = np.gradient(arr, .5, axis=1)
    arr_x2 = gradient_n(arr, 1, .5, axis=1)
    arr_xx = np.diff(arr, 2, axis=1) / .5**2
    arr_xx = np.pad(arr_xx, ((0, 0), (1, 1)), 'constant', constant_values=0)
    arr_xx2 = gradient_n(arr, 2, .5, axis=1)

    assert np.sum(arr_x - arr_x2) == 0
    assert np.sum(arr_xx - arr_xx2) == 0

    fig, axs = plt.subplots(2, 2, figsize=(29, 21))
    axs = np.array(axs).flatten()

    ax = axs[0]
    ax.set_title('x-cut')
    ax.plot(x, arr[0, :], marker='o', label='arr')
    ax.plot(x, arr_x[0, :], marker='o', label='arr_x')
    ax.plot(x, arr_x2[0, :], marker='x', label='arr_x2', ls='--')
    ax.plot(x, arr_xx[0, :], marker='o', label='arr_xx')
    ax.plot(x, arr_xx2[0, :], marker='x', label='arr_xx2', ls='--')
    ax.legend()

    ax = axs[1]
    ax.set_title('arr')
    im = ax.imshow(arr, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

    ax = axs[2]
    ax.set_title('arr_x')
    im = ax.imshow(arr_x, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

    ax = axs[3]
    ax.set_title('arr_xx')
    im = ax.imshow(arr_xx, cmap='bwr')
    cbar = ax.figure.colorbar(im, ax=ax, pad=.05)

test_gradient_n()

这是从原始文档中摘录的内容（在撰写本文时，可在http://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html中找到）。它指出，除非采样距离为1，否则需要包含一个包含距离的列表作为参数。

numpy.gradient(f, *varargs,**kwargs)

返回N维数组的渐变。
梯度是使用内部的二阶精确中心差和边界处的一阶差或二阶精确单侧（向前或向后）差计算的。因此，返回的梯度具有与输入数组相同的形状。
参数：
f：array_like包含标量函数样本的N维数组。
varargs：标量列表，可选的N个标量，指定每个维度的采样距离，即dx、dy、dz...默认距离：1.
边顺序：{1，2}，使用边界处的N阶精度差计算可选梯度。默认值：1. 1.9.1版中的新增功能。
退货：
gradient：ndarray与f形状相同的N个数组，给出f在每个维度上的导数。

我的解决方案是创建一个类似于np.gradient的函数，该函数从数组数据中计算二阶导数。

import numpy as np

def gradient2_even(y, h=None, edge_order=1):
    """
    Return the 2nd-order gradient i.e. 
    2nd derivatives of y with n samples and k components.

    The 2nd-order gradient is computed using second-order-accurate central differences
    in the interior points and either first or second order accurate one-sided
    (forward or backwards) differences at the boundaries.
    The returned gradient hence has the same shape as the input array.

    Parameters
    ----------
    y : 1d or 2d array_like
        The array containing the samples. If 2d with shape (n,k),
        n is the number of samples at least 2 while k is the number of
        y series/components. 1d input is equivalent to 2d input with shape (n,1).
    h : constant or 1d, optional
        spacing between the y samples. Default unitary spacing for
        all y components. Spacing can be specified using:

            1. Single scalar spacing value for all y components.
            2. 1d array_like of length k specifying the spacing for each y component

    edge_order : {1, 2}, optional
                 Order 1 means 3-point forward/backward finite differences
                 are used to calculate the 2nd derivatves at the edge points while
                 order 2 uses 4-point forward/backward finite differences.

    Returns
    ----------
    d2y : 1d or 2d array
        Array containing the 2nd derivatives. The output shape is the same as y.
    """
    if edge_order!=1 and edge_order!=2:
        raise ValueError('edge_order must be 1 or 2.')
    else:
        pass

    y = np.asfarray(y)
    origshape = y.shape
    if y.ndim!=1 and y.ndim!=2:
        raise ValueError('y can only be 1d or 2d.')
    elif y.ndim==1:
        y = np.atleast_2d(y).T
    elif y.ndim==2:
        if y.shape[0]<2:
            raise ValueError('The number of y samples must be atleast 2.')
        else:
            pass
    else:
        pass
    n,k = y.shape

    if h is None:
        h = 1.0
    else:
        h = np.asfarray(h)
        if h.ndim!=0 and h.ndim!=1:
            raise ValueError('h can only be 0d or 1d.')
        elif h.ndim==0:
            pass
        elif h.ndim==1 and h.size!=n:
            raise ValueError('If h is 1d, it must have the same number as the components of y.')
        else:
            pass
    d2y = np.zeros_like(y)
    if n==2:
        pass
    elif n==3:
        d2y[:] = ( 1/h**2 * (y[0] - 2*y[1] + y[2]) )
    else:
        d2y = np.zeros_like(y)
        d2y[1:-1]=1/h**2 * ( y[:-2] - 2*y[1:-1] + y[2:] )
        if edge_order==1:
            d2y[0]=1/h**2 * ( y[0] - 2*y[1] + y[2] )
            d2y[-1]=1/h**2 * ( y[-1] - 2*y[-2] + y[-3] )
        else:
            d2y[0]=1/h**2 * ( 2*y[0] - 5*y[1] + 4*y[2] - y[3] )
            d2y[-1]=1/h**2 * ( 2*y[-1] - 5*y[-2] + 4*y[-3] - y[-4] )
    return d2y.reshape(origshape)

以您的示例为例，

# After importing the function from the script file or running it

from numpy import *
from matplotlib.pyplot import *

x, h = linspace(0, 10, 17) # use a fairly coarse grid to see the discrepancies better
y = sin(x)
ypp = -sin(x) # analytical 2nd derivatives

# Compute numerically the 2nd derivatives using 2nd-order finite differences at the edge points

d2y = gradient2_even(y, h, 2)

# Compute numerically the 2nd derivatives using nested gradient function

d2y2 = gradient(gradient(y, h, edge_order=2), h, edge_order=2)

# Compute numerically the 2nd derivatives using 1st-order finite differences at the edge points

d2y3 = gradient2_even(y, h, 1)

fig,ax=subplots(1,1)
ax.plot(x, ypp, x, d2y, 'o', x, d2y2, 'o', x, d2y3, 'o'), ax.grid()
ax.legend(['Analytical', 'edge_order=2', 'nested gradient', 'edge_order=1'])
fig.tight_layout()