Why my np.gradient calculation in R^2 doesn't fit with the analytical gradient calculation?

Question

I'm trying to compute a gradient on a map using np.gradient, but I'm encountering issues. To simplify my problem I am trying on an analytical function

z = f(x,y) = -(x - 2)**2 - (y - 2)**2

np.gradient is not providing the expected results; the vectors should point towards the center. What am I doing wrong?

Here is the code that I am running:

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)
Z = -(X - 2)**2 - (Y - 2)**2

# Compute gradients numerically
dz_dx, dz_dy = np.gradient(Z, x, y)


# Downsampling to reduce the density of arrows
step = 10

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[::step, ::step], Y[::step, ::step], dz_dx[::step, ::step], dz_dy[::step, ::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()

Answer 1

the problem is not in the gradient function, it is in the different indexing order of np.meshgrid and np.gradient.

by default np.gradient assumes the indexing is the same order as the arguments, ie:

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

whereas np.meshgrid default indexing results in the opposite indexing,

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

you didn't notice this bug because both X and Y are identical, if you made X and Y have different number of points you would get an error.

I like to use Z[y,x], so just swap the order of arguments and return of np.gradient and you will get the correct result.

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

Answer 2

I'm not sure why np.gradient is not working as your expect it, but you could instead manually calculate the gradients, e.g.,

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()