np.gradient()的edge_order参数作用及Python替代Matlab del2函数的方法

1. numpy.gradient函数中edge_order参数的具体作用 #

• edge_order参数用于控制数组边缘位置的梯度计算精度，仅对边缘点生效，可选值为1或2：
  • 当edge_order=1时，边缘点采用一阶差分计算：前边缘用向前差分（(y[1]-y[0])/h），后边缘用向后差分（(y[-1]-y[-2])/h），精度较低。
  • 当edge_order=2时，边缘点采用二阶差分计算：前边缘使用(-3*y[0]+4*y[1]-y[2])/(2*h)，后边缘使用(3*y[-1]-4*y[-2]+y[-3])/(2*h)，精度更高。
• 数组内部的非边缘点，无论edge_order取何值，均采用中心差分计算，精度一致。

2. Python中是否有与Matlab del2函数等价的实现？ #

有两种可靠的实现方式：

手动基于numpy实现 #

Matlab的del2计算的是离散拉普拉斯算子的1/4，对应二维数组的实现如下（支持自定义步长）：

import numpy as np

def del2(arr, dx=1.0, dy=1.0):
    # 初始化拉普拉斯计算
    laplacian = (np.roll(arr, 1, axis=0) + np.roll(arr, -1, axis=0) - 2*arr) / dx**2 + \
                (np.roll(arr, 1, axis=1) + np.roll(arr, -1, axis=1) - 2*arr) / dy**2
    
    # 修正边缘（解决roll带来的wrap问题，匹配Matlab默认边界处理）
    laplacian[0, :] = (arr[1, :] - arr[0, :])/dx**2 + (arr[0, 1:] + arr[0, :-1] - 2*arr[0, :])/dy**2
    laplacian[-1, :] = (arr[-2, :] - arr[-1, :])/dx**2 + (arr[-1, 1:] + arr[-1, :-1] - 2*arr[-1, :])/dy**2
    laplacian[:, 0] = (arr[:, 1] - arr[:, 0])/dy**2 + (arr[1:, 0] + arr[:-1, 0] - 2*arr[:, 0])/dx**2
    laplacian[:, -1] = (arr[:, -2] - arr[:, -1])/dy**2 + (arr[1:, -1] + arr[:-1, -1] - 2*arr[:, -1])/dx**2
    
    return laplacian / 4.0

使用scipy简化实现 #

scipy.ndimage.laplace计算标准离散拉普拉斯，将结果除以4即可与Matlabdel2等价：

from scipy.ndimage import laplace
import numpy as np

def del2_scipy(arr):
    return laplace(arr) / 4.0

注：如果数组基于非单位步长的网格，需要在手动实现中传入对应的dx和dy参数，确保结果与Matlab一致。

内容的提问来源于stack exchange，提问作者Roberto_1986