I'm trying to use the np.gradient for a multidimensional irregularly spaced array. I've got a naive solution working, but it requires a for loop which will be a bottleneck for real computations.
import numpy as np
x = np.geomspace(2*np.pi, 8*np.pi, 500)[None,:]
y = np.geomspace(2*np.pi, 8*np.pi, 500)[:,None]
XX, YY = np.meshgrid(x,y)
f = np.sin(x) * np.cos(y)
dfdx = np.cos(x) * np.cos(y)
grad = np.zeros_like(f)
for row in range(f.shape[0]):
grad[row] = np.gradient(f[row], XX[row], edge_order=2)
print(np.linalg.norm(grad-dfdx, ord=2))
Is there a way to do this computation using np.apply_along_axis? The reason for needing this because np.gradient does not accept multidimensional arrays for the distance argument, they must either be scalars or a 1-D array.