Ваш собственный ответ - шаг в правильном направлении, но есть некоторые проблемы как в ответе, так и в коде, генерирующем декартовы производные.
В этих строках есть проблема:
x = np.linspace(start, end, num=n)
dx = (end - start) / n
Размер шага на самом деле (end-start)/(n-1).
Здесь:
x_grid, y_grid, z_grid = np.meshgrid(x, y, z)
df_dx[x, y, z] = (eval_xyz[x + 1, y, z] - eval_xyz[x, y, z]) / dx
вы попали в ловушку настройки meshgrid по умолчанию: meshgrid(np.arange(n1), np.arange(n2)) вернет массивы в форма (n2, n1), если вы не добавите параметр indexing='ij'. Поскольку у вас размер n во всех измерениях, вы не получите предупреждений об ошибках индексации, но вы, возможно, потратите много времени, пытаясь отладить, почему числа не имеют смысла.
Когда вы манипулируете многомерными массивов, рекомендуется установить размеры в разных направлениях на несколько разные значения, чтобы вы могли легко проверить, соответствуют ли формы массива тем, которые вы хотите.
Кроме того, вы должны обычно оценивать производную как (f[i+1]-f[i-1])/(2*dx), что верно до второго порядка по x.
for x in range(eval_xyz.shape[0] - 1):
for y in range(eval_xyz.shape[1] - 1):
for z in range(eval_xyz.shape[2] - 1):
При работе с numpy вы всегда должны пытаться векторизовать операции, а не записывать для циклов, которые потенциально необходимо повторять более тысячи элементов.
Вот код, который вычисляет декартову производную, а затем радиальную производную.
import numpy as np
def get_cartesian_gradient(f, xyzsteps):
"""For f shape (nx, ny, nz), return gradient as (3, nx, ny, nz) shape.
xyzsteps is a (3,) array.
Note: edge points of the gradient array are set to NaN.
(Exercise for the reader to implement those).
"""
fshape = f.shape
grad = np.full((3,) + fshape, np.nan, dtype=np.float64)
sl, sm, sr = slice(0, -2), slice(1, -1), slice(2, None)
# Note: multiplying is faster than dividing.
grad[0, sm, sm, sm] = (f[sr, sm, sm] - f[sl, sm, sm]) * (0.5/xyzsteps[0])
grad[1, sm, sm, sm] = (f[sm, sr, sm] - f[sm, sl, sm]) * (0.5/xyzsteps[1])
grad[2, sm, sm, sm] = (f[sm, sm, sr] - f[sm, sm, sl]) * (0.5/xyzsteps[2])
return grad
def get_dfdr_from_cartesian(grad, x1s, y1s, z1s):
"""Return df/dr array from gradient(f).
grad.shape must be (3, nx, ny, nz)
return shape (nx, ny, nz).
"""
_, nx, ny, nz = grad.shape
# we need sin(theta), cos(theta), sin(phi), and cos(phi)
# rxy: shape (nx, ny, 1)
rxy = np.sqrt(x1s.reshape(-1, 1, 1)**2 + y1s.reshape(1, -1, 1)**2)
# r: shape (nx, ny, nz)
r = np.sqrt(rxy**2 + z1s.reshape(1, 1, -1)**2)
# change zeros to NaN
r = np.where(r==0, np.nan, r)
rxy = np.where(rxy==0, np.nan, rxy)
cos_theta = z1s.reshape(1, 1, -1) / r
sin_theta = rxy / r
cos_phi = x1s.reshape(-1, 1, 1) / rxy
sin_phi = y1s.reshape(1, -1, 1) / rxy
# and the derivative
dfdr = (grad[0]*cos_phi + grad[1]*sin_phi)*sin_theta + grad[2]*cos_theta
return dfdr
x1s = np.linspace(-1, 1, 19)
y1s = np.linspace(-1, 1, 21)
z1s = np.linspace(-1, 1, 23)
xs, ys, zs = np.meshgrid(x1s, y1s, z1s, indexing='ij')
xyzsteps = [x1s[1]-x1s[0], y1s[1]-y1s[0], z1s[1]-z1s[0]]
def func(x, y, z):
return x**2 + y**2 + z**2
def dfdr_analytical(x, y, z):
r = np.sqrt(x**2 + y**2 + z**2)
return 2*r
# grad has shape (3, nx, ny, nz)
grad = get_cartesian_gradient(func(xs, ys, zs), xyzsteps)
dfdr = get_dfdr_from_cartesian(grad, x1s, y1s, z1s)
# test
diff = dfdr - dfdr_analytical(xs, ys, zs)
assert np.nanmax(np.abs(diff)) < 1e-14
Обратите внимание, что я решил возвращать значения NaN для точек на z- оси, потому что df / dr там не определено, если только f (x, y, z) не является вращательно симметричным c вокруг оси z и имеет df / dr = 0 во всех направлениях. Это то, что не гарантируется для произвольного набора данных.
Причина замены нулей в знаменателях на np.nan с использованием np.where заключается в том, что при делении на ноль будут появляться предупреждающие сообщения, а при делении на nan выиграет. т.