Почему мой расчет np.gradient в R^2 не соответствует аналитическому расчету градиента?

Я не уверен, почему np.gradient не работает так, как вы ожидаете, но вместо этого вы можете вручную рассчитать градиенты, например,

import numpy as np
import matplotlib.pyplot as plt

# Define the grids for x and y
x = np.linspace(0, 4, 100)  # 100 points between 0 and 4
y = np.linspace(0, 4, 100)  # 100 points between 0 and 4
X, Y = np.meshgrid(x, y)  # Create a 2D grid

# Define the function f(x, y)
def func(X, Y):
    return -(X - 2)**2 - (Y - 2)**2

def grad(func, X, Y, delta=0.01):
    """
    Calculate gradients of a function.
    """

    dp = delta / 2
    df_dx = (func(X + dp, Y) - func(X - dp, Y)) / delta
    df_dy = (func(X, Y + dp) - func(X, Y - dp)) / delta

    return df_dx, df_dy

Z = func(X, Y)
dz_dx, dz_dy = grad(func, X, Y)

# Downsampling to reduce the density of arrows
step = 10
start = step // 2

plt.figure(figsize=(10, 8))
contour = plt.contourf(X, Y, Z, cmap='viridis', levels=50, alpha=0.8)
plt.colorbar(contour, label='f(x, y)')
plt.quiver(X[start::step, start::step], Y[start::step, start::step], dz_dx[start::step, start::step], dz_dy[start::step, start::step], 
           color='r', headlength=3, headwidth=4)
plt.title('Function $f(x, y) = -(x - 2)^2 - (y - 2)^2$ and its gradients (numerical)')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)

plt.show()

проблема не в функции градиента, а в разном порядке индексации np.meshgrid и np.gradient.

по умолчанию np.gradient предполагает, что индексация имеет тот же порядок, что и аргументы, т.е.:

Z[x,y] -> np.gradient(Z, x, y)

тогда как np.meshgrid индексация по умолчанию приводит к противоположной индексации,

Z[y,x] -> np.meshgrid(x,y)  # default indexing = 'xy'

вы не заметили эту ошибку, потому что X и Y идентичны. Если бы вы сделали X и Y с разным количеством точек, вы получили бы ошибку.

Мне нравится использовать Z[y,x], поэтому просто поменяйте порядок аргументов и возврат np.gradient, и вы получите правильный результат.

dz_dy, dz_dx = np.gradient(Z, y, x)  # Z[y,x]

Как объяснил Ахмед АЕК в своем ответе, проблема в заказах. Вместо того, чтобы менять способ звонка np.gradient, вам следует изменить способ звонка np.meshgrid.

Если вы посмотрите документацию np.meshgrid, вы увидите это по умолчанию indexing = "xy". Другой вариант indexing = "ij". Это контролирует, как определяются/упорядочиваются точки в сетке.

Для "xy" упорядочения вы можете представить массив как декартову координатную плоскость. Внизу слева — (0,0), значение x увеличивается при движении вправо, а значение y увеличивается при движении вверх.

Для порядка "ij" индексация аналогична координатам декартовой плоскости. Это означает, что (0,0), когда он используется в качестве индекса массива, находится вверху слева, увеличение значения x будет означать увеличение первого значения в индексе массива, поэтому перемещение вниз (построчно) увеличивает значение x. Второй индекс — это «координата» y, поэтому перемещение по столбцу (столбец за столбцом) увеличивает значение y.

Итак, все, что вам нужно сделать, это изменить вызов сетки на:

X, Y = np.meshgrid(x, y, indexing = "ij")