NumPy gradient over multidimensional irregularly spaced array

Viewed 75

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.

1 Answers

You can do something like the following:

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.array(
    [np.gradient(f[row], XX[row], edge_order=2) for row in range(f.shape[0])]
)

print(np.linalg.norm(grad-dfdx, ord=2))

which returns 0.08971829701854393.

Related