Gradient function on R and Python

Question

I was testing a code in R and then I wanted to use a similar code in Python. I have this potential:Z = 1/sqrt(Y^2+X^2+2*d*X+d^2). I tried to plot it with R using the pracma package like this:

require(pracma)

range = seq(-15,15,by=1)
mesh = meshgrid(range,range)

X = mesh$X
Y = mesh$Y    
d = 5

# Now I define the potential:
Z = 1/sqrt(Y^2+X^2+2*d*X+d^2)
contour(range,range,t(Z),col="black",xlab="x",ylab="y")

# Now the gradient:
grX = gradient(Z,range,range)$X
grY = gradient(Z,range,range)$Y

# I wanted the vectors to be normalized so I did this:
grXN = grX/sqrt(grX^2+grY^2)
grYN = grY/sqrt(grX^2+grY^2)

# Finally I draw the vector field:
quiver(X,Y,grXN,grYN,scale=0.5,col="black")

When I run this code I get this: Quiver in R Which is more or less what I wanted, except it is a bit ugly.

Then I made this code in Python like this:

import numpy as np
import matplotlib.pyplot as plt

rng = np.arange(-15,15,1)
X,Y = np.meshgrid(rng,rng)
d = 5

# Now I define the potential:
Z = 1/np.sqrt(Y**2+X**2+2*d*X+d**2)

# Now the gradient:
grX,grY = np.gradient(Z)

# Since I want the vectors normalized I did this:
grXN = grX/np.sqrt(grX**2+grY**2)
grYN = grY/np.sqrt(grX**2+grY**2)

# In this case I made a different contour:
contor = plt.contourf(X,Y,Z,50)
cbar = plt.colorbar(contor)

# Finally the arrows:
Q = plt.quiver(X,Y,grXN,grYN,color="w",scale=30)

When I run this code I get this: Quiver in Python Which is cute, but totally different from the result obtained from R. Why?

Answer 1

Your Python arrows are simply twisted through 90° relative to your R arrows. The reason is the difference in convention between numpy and matplotlib when it comes to the order of the dimensions "X" and "Y". For reasons of (debatable) intuitiveness, matplotlib functions—like those of its mad uncle, Matlab—typically ask for X, then Y, in that order. Likewise, taking its inspiration from Matlab's well-worn meshgrid function, numpy.meshgrid also adopts that convention, in that its indexing argument defaults to 'xy'. But I believe numpy.gradient must be following the more usual numpy-esque assumption of 'ij' indexing (and when visualizing arrays, Y is usually the row-to-row 'i' dimension, which in numpy arrays is dimension 0, so then X is usually the column-to-column 'j' dimension, dimension 1). Therefore, you should actually be saying:

grY,grX = np.gradient(Z)

instead of

grX,grY = np.gradient(Z)