title: Correctly compute the divergence of a vector field in Python
body: <p>I am trying to compute the divergence of a vector field:</p>
  <pre><code>Fx  = np.cos(xx + 2*yy)
  Fy  = np.sin(xx - 2*yy)
  F = np.array([Fx, Fy])
  </code></pre>
  <p><a href="https://i.stack.imgur.com/FgC6T.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/FgC6T.png" alt="enter image description here" /></a></p>
  <h1>Analytic Solution</h1>
  <p>This is how the divergence (div(F) = dF/dx + dF/dy ) should look like, based on the analytic computation of the divergence (see Wolfram Alpha <a href="https://www.wolframalpha.com/input/?i=divergence%20of%20Cos%5Bx%2B2*y%5D%2CSin%5Bx-2*y%5D" rel="nofollow noreferrer">here</a>):</p>
  <ul>
  <li>dFx/dx = d/dx cos(x+2y) = -sin(x+2y)</li>
  <li>dFy/dy = d/dy sin(x-2y) = -2*cos(x-2y)</li>
  </ul>
  <p>The divergence:</p>
  <pre><code>div_analy = -np.sin(xx + 2*yy) - 2*np.cos(xx - 2*yy)
  </code></pre>
  <p><a href="https://i.stack.imgur.com/4TvUY.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4TvUY.png" alt="enter image description here" /></a></p>
  <p>The code:</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  # Number of points (NxN)
  N = 50
  # Boundaries
  ymin = -2.; ymax = 2.
  xmin = -2.; xmax = 2.
  
  
  # Create Meshgrid
  x = np.linspace(xmin,xmax, N)
  y = np.linspace(ymin,ymax, N)
  xx, yy = np.meshgrid(x, y)
  
  # Analytic computation of the divergence (EXACT) 
  div_analy = -np.sin(xx + 2*yy) - 2*np.cos(xx - 2*yy)
  
  # PLOT
  plt.imshow(div_analy , extent=[xmin,xmax,ymin,ymax],  origin=&quot;lower&quot;, cmap=&quot;jet&quot;)
  </code></pre>
  <h1>Numeric Solution</h1>
  <p>Now, I am trying to obtain the same numerically, so I used <a href="https://stackoverflow.com/questions/67970477/compute-divergence-with-python/67971515#67971515">this function</a> to compute the divergence</p>
  <pre><code>def divergence(f,sp):
      &quot;&quot;&quot; Computes divergence of vector field 
      f: array -&gt; vector field components [Fx,Fy,Fz,...]
      sp: array -&gt; spacing between points in respecitve directions [spx, spy,spz,...]
      &quot;&quot;&quot;
      num_dims = len(f)
      return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])
  </code></pre>
  <p>When I plot the divergence using this function:</p>
  <pre><code> # Compute Divergence
   points = [x,y]
   sp = [np.diff(p)[0] for p in points]
   div_num = divergence(F, sp)
      
   # PLOT
   plt.imshow(div_num, extent=[xmin,xmax,ymin,ymax], origin=&quot;lower&quot;,  cmap=&quot;jet&quot;)
  </code></pre>
  <p>... I get:</p>
  <p><a href="https://i.stack.imgur.com/C0c26.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/C0c26.png" alt="enter image description here" /></a></p>
  <h1>The issue</h1>
  <p>The numeric solution is different from the analytic solution! What am I doing wrong?</p>