Correctly compute the divergence of a vector field in Python

Question

I am trying to compute the divergence of a vector field:

Fx  = np.cos(xx + 2*yy)
Fy  = np.sin(xx - 2*yy)
F = np.array([Fx, Fy])

Analytic Solution

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 here):

• dFx/dx = d/dx cos(x+2y) = -sin(x+2y)
• dFy/dy = d/dy sin(x-2y) = -2*cos(x-2y)

The divergence:

div_analy = -np.sin(xx + 2*yy) - 2*np.cos(xx - 2*yy)

The 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="lower", cmap="jet")

Numeric Solution

Now, I am trying to obtain the same numerically, so I used this function to compute the divergence

def divergence(f,sp):
    """ Computes divergence of vector field 
    f: array -> vector field components [Fx,Fy,Fz,...]
    sp: array -> spacing between points in respecitve directions [spx, spy,spz,...]
    """
    num_dims = len(f)
    return np.ufunc.reduce(np.add, [np.gradient(f[i], sp[i], axis=i) for i in range(num_dims)])

When I plot the divergence using this function:

 # 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="lower",  cmap="jet")

... I get:

The issue

The numeric solution is different from the analytic solution! What am I doing wrong?

Answer 1

Both graphs are wrong, because you use np.meshgrid the wrong way.

The other parts of your code are expecting xx[a, b], yy[a, b] == x[a], y[b], where a, b are integers between 0 and 49 in your case.

On the other hand, you write

xx, yy = np.meshgrid(x, y)

which causes xx[a, b], yy[a, b] == x[b], y[a]. Futhermore, the value of div_analy[a, b] becomes -sin(x[b]+2y[a]) - 2cos(x[b]+2y[a]) and the value of div_num[a, b] also becomes some other wrong value.

You can simply fix it by writing

yy, xx = np.meshgrid(x, y)

Then you will see both solutions get the correct result.

By the way, you set N = 50, so that np.linspace(xmin,xmax, N) takes 50 points from interval [-2, 2] and the distance between contiguous points will be 1/49. Setting N=51 may get a more radable line sapce.