In Octave, I have
x = -0.8:0.2:1;
y = -0.8:0.2:1;
z = -0.8:0.8:1;
[xx,yy,zz] = meshgrid(x, y, z);
u = sin(pi .* xx) .* cos(pi .* yy) .* cos(pi .* zz);
dx = xx(1,:,1)(:);
dy = yy(:,1,1)(:);
dz = zz(1,1,:)(:);
[a, b, c] = gradient (u, dx, dy, dz);
b(1,:,1)
which outputs
1.18882 1.92355 1.92355 1.18882 0.00000 -1.18882 -1.92355 -1.92355 -1.18882 -0.00000
With Python, I tried to replicate it,
import numpy as np
xx, yy, zz = np.meshgrid(np.arange(-0.8, 1.2, 0.2),
np.arange(-0.8, 1.2, 0.2),
np.arange(-0.8, 1.2, 0.8))
u = np.sin(np.pi * xx) * np.cos(np.pi * yy) * np.cos(np.pi * zz)
dx = xx[0,:,0]
dy = yy[:,0,0]
dz = zz[0,0,:]
a,b,c = np.gradient (u, dx, dy, dz)
print (b[0,:,0])
I get
[-1.18882065 -0.59441032 0.59441032 1.55618643 1.92355221 1.55618643
0.59441032 -0.59441032 -1.55618643 -1.92355221]
which looks different. How do I get the Numpy gradient call match the Octave results?
My Numpy version is 1.19.0, Octave is 4.2.2.
I saw this post asking another Matlab / Python gradient question, that call works for me, which I guess is the low dimensional case.
So you need to inverse the axis order of the first two dimensions:
a,b,c = np.gradient (u, dx, dy, dz, axis=[1,0,2])