numpy.gradient可以对函数进行求导,在设置edge_order=1时其逆操作可以由以下代码实现。改写自stackoverflow的回答。该程序仅适用于一维导数且dx为常数的情况,但是可以指定axis:
import numpy as np
def integrate(dydx, y0, dx=1, axis=0):
'''
Inverse of numpy.gradient, only for 1d, evenly spaced x and edge_order=1.
dydx: array_like, (N1, N2, ..., N_axis, ...), output of gradient
y0: array_like, (N1, N2, ...), the initial value, one dimension less than dydx
dx: float, the dx used in np.gradient
axis: int, the axis used in np.gradient
return: ndarray, the same as the data before np.gradient
'''
dydx = np.asarray(dydx)
dydx = np.moveaxis(dydx, axis, 0)
y1 = np.append(dydx[:1]*0., dydx[1:-1:2].cumsum(axis=0), axis=0)
y2 = dydx[::2].cumsum(axis=0) - dydx[:1] / 2
y1 = np.expand_dims(y1, 1)
y2 = np.expand_dims(y2, 1)
y3 = np.hstack([y1, y2])
y3 = y3.reshape(-1, *y1.shape[2:])
out = y0 + 2 * y3[:dydx.shape[0]] * dx
out = np.moveaxis(out, 0, axis)
return out
if __name__ == '__main__':
y = np.random.rand(2, 10, 3)
dydx = np.gradient(y, 0.143, edge_order=1, axis=1)
y1 = integrate(dydx, y[:,0,:], 0.143, axis=1)
print(y)
print(y1)
print('=============')
print(np.abs(y1-y).max())



