Градиент от градиента функции в Python

Question

Имеетя поле U, представленное в виде трёхмерного массива mnk. Подскажите, как в Python посчитать следующее выражение:

где k - массив коэффициентов такой же размерности, как и U.

Очевидно, что на выходе должен получиться массив размера mnk, т.к. двойной градиент преобразуется в лапласиан:

Применение дважды метода np.gradient(k*np.gradient(U)) не даёт нужного результата, т.к. не понимает этого и продолжает работать с вектором градиента как с массивом поля, просто большей размерности.

Answer 1

Если у вас в первой формуле вашего вопроса ошибка, и вместо двукратного оператора Лапласа там стоит двукратная набла (оператор градиента), то вашу формулу можно упростить вот так:

Для градиента есть встроенная функция, для скалярного произведения тоже. Вам остаётся только вычислить оператор Лапласа. Для этого нужно вручную вычислить градиент от каждого компонента градиента и взять из результата соответствующую производную.

def laplacian(f, dx):
    """
    Вычисляет значение оператора Лапласа на сетке.
    
    f - значения функции в n-мерной сетке.
    dx - кортеж, элементы кортежа есть шаги сетки.
    """
    result = np.zeros_like(f)
    g = np.gradient(f, *dx, edge_order=2)
    for i, g_i in enumerate(g):
        gg_i = np.gradient(g_i, *dx, edge_order=2)
        result += gg_i[i]
    return result

Для примера вычислим оператор Лапласа вручную и функцией laplacian. В качестве подопытной свинки будет функция

import numpy as np
# Сетка
dx = 0.01
X,Y,Z = np.mgrid[-1:1:dx,-1:1:dx,-1:1:dx]
# подопытная функция
F = np.exp(-(X*X+2*Y*Y+3*Z*Z))
# Производные
dFdx = -2*X*F
d2Fdx2 = F*(4*X*X-2)

dFdy = -4*Y*F
d2Fdy2 = F*(16*Y*Y - 4)

dFdz = -6*Z*F
d2Fdz2 = F*(36*Z*Z - 6)
# Лапласиан
L_manual = d2Fdx2 + d2Fdy2 + d2Fdz2
L_computed = laplacian(F, (dx,dx,dx))

# среднеквадратичное отклонение вычисленного от точного
print(np.std(L_manual - L_computed))

Результат 0.0026064524156072335

В основном ошибка набирается по краям. Если отбросить полосу вдоль краёв сетки, то ошибка уменьшается в два раза: np.std(L_manual[3:-3,3:-3,3:-3] - L_computed[3:-3,3:-3,3:-3]) -> 0.0011257917625283698

БЛОГ НА HUSL

Градиент от градиента функции в Python

Ответы (1 шт):