Компоненты numpy.gradient симметричной функции различны

Обычно градиенты и якобианы являются операторами функций

Если вам нужен градиент f = 1/(X*X + Y*Y +1.0), вы должны вычислить его символически. Или оцените его с помощью численных методов, использующих эту функцию.

Я не знаю, что такое градиент постоянного трехмерного массива. numpy.gradient — это одномерное понятие.

Python имеет пакет sympy, который может автоматически вычислять якобианы символически.

Если под second order derivative of a scalar 3d field вы подразумеваете лапласиан, то вы можете оценить это с помощью стандартного 4-точечного трафарета.

Давайте пройдемся по этому шаг за шагом. Итак, во-первых, как правильно упомянул meowgoesthedog numpy вычисляет производные в направлении.

Способ Numpy для вычисления градиентов

Важно отметить, что np.gradient использует значение центральных различий (для простоты мы смотрим только в одном направлении):

grad_f[i] = (f[i+1] - f[i])/2 + (f[i] - f[i-1])/2 =  (f[i+1] - f[i-1])/2

На границе numpy вычисляет (в качестве примера возьмем min)

grad_f[min] = f[min+1] - f[min]
grad_f[max] = f[max] - f[max-1]

В вашем случае граница 0 и 2.

2D случай

Если вы используете более одного измерения, нам необходимо учитывать направление производной. np.gradient вычисляет производные во всех возможных направлениях. Давайте воспроизведем ваши результаты:

Давайте двигаться вдоль столбцов, так что мы считаем с строки векторов

f[1,:] - f[0,:] 

Вывод

array([0.16666667, 0.5       , 0.16666667])

это точно первая строка первого элемента вашего градиента.

Строка рассчитывается с центрированными производными, поэтому:

(f[2,:]-f[1,:])/2 + (f[1,:]-f[0,:])/2

Вывод

array([0., 0., 0.])

Третий ряд:

f[2,:] - f[1,:] 

Вывод

array([-0.16666667, -0.5       , -0.16666667])

Для другого направления просто поменяйте местами : и числа и помните, что сейчас вы рассчитываете столбцы векторов. Это приводит непосредственно к производной транспонированный в случае симметричной функции, как в вашем случае.

3D чехол

x_ = (-1,0,4)
y_ = (-3,0,1)
z_ = (-1,0,12)

x, y, z = np.meshgrid(x_, y_, z_, indexing='ij')
f = 1/(x**2 + y**2 + z**2 + 1)
np.gradient(f)[1]

Вывод

array([[[ *2.50000000e-01,  4.09090909e-01,  3.97702165e-04*],
        [ 8.33333333e-02,  1.21212121e-01,  1.75554093e-04],
        [-8.33333333e-02, -1.66666667e-01, -4.65939801e-05]],

       [[ **4.09090909e-01,  9.00000000e-01,  4.03045231e-04**],
        [ 1.21212121e-01,  2.00000000e-01,  1.77904287e-04],
        [-1.66666667e-01, -5.00000000e-01, -4.72366556e-05]],

       [[ ***1.85185185e-02,  2.03619910e-02,  3.28827183e-04***],
        [ 7.79727096e-03,  8.54700855e-03,  1.45243282e-04],
        [-2.92397661e-03, -3.26797386e-03, -3.83406181e-05]]])

Приведенный здесь градиент рассчитывается по строкам (0 будет по матрицам, 1 по строкам, 2 по столбцам).

Это можно рассчитать по

(f[:,1,:] - f[:,0,:])

Вывод

array([[*2.50000000e-01, 4.09090909e-01, 3.97702165e-04*],
       [**4.09090909e-01, 9.00000000e-01, 4.03045231e-04**],
       [***1.85185185e-02, 2.03619910e-02, 3.28827183e-04***]])

Я добавил звездочки, чтобы было понятно, где найти соответствующие векторы-строки. Поскольку мы рассчитали градиент в направлении 1, нам нужно искать строки векторов.

Если кто-то хочет воспроизвести весь градиент, это делается с помощью

np.stack(((f[:,1,:] - f[:,0,:]), (f[:,2,:] - f[:,0,:])/2, (f[:,2,:] - f[:,1,:])), axis=1)

n-мерный случай

Мы можем обобщить то, что мы узнали здесь, для вычисления градиентов произвольных функций по направлениям.

def grad_along_axis(f, ax):
    f_grad_ind = []
    for i in range(f.shape[ax]):
        if i == 0:
            f_grad_ind.append(np.take(f, i+1, ax) - np.take(f, i, ax))
        elif i == f.shape[ax] -1:
            f_grad_ind.append(np.take(f, i, ax) - np.take(f, i-1, ax))
        else:
            f_grad_ind.append((np.take(f, i+1, ax) - np.take(f, i-1, ax))/2)
    f_grad = np.stack(f_grad_ind, axis=ax)
    return f_grad

где

np.take(f, i, ax) = f[:,...,i,...,:]

и i находится под индексом ax.