Integrate a 2D vectorfield-array (reversing np.gradient)

Question

i have the following problem: I want to integrate a 2D array, so basically reversing a gradient operator.

Assuming i have a very simple array as follows:

shape = (60, 60)
sampling = 1
k_mesh = np.meshgrid(np.fft.fftfreq(shape[0], sampling), np.fft.fftfreq(shape[1], sampling))

Then i construct my vectorfield as a complex-valued arreay (x-vector = real part, y-vector = imaginary part):

k = k_mesh[0] + 1j * k_mesh[1]

So the real part for example looks like this

Now i take the gradient:

k_grad = np.gradient(k, sampling)

I then use Fourier transforms to reverse it, using the following function:

def freq_array(shape, sampling):

    f_freq_1d_y = np.fft.fftfreq(shape[0], sampling[0])
    f_freq_1d_x = np.fft.fftfreq(shape[1], sampling[1])
    f_freq_mesh = np.meshgrid(f_freq_1d_x, f_freq_1d_y)
    f_freq = np.hypot(f_freq_mesh[0], f_freq_mesh[1])

    return f_freq


def int_2d_fourier(arr, sampling):
    freqs = freq_array(arr.shape, sampling)

    k_sq = np.where(freqs != 0, freqs**2, 0.0001)
    k = np.meshgrid(np.fft.fftfreq(arr.shape[0], sampling), np.fft.fftfreq(arr.shape[1], sampling))

    v_int_x = np.real(np.fft.ifft2((np.fft.fft2(arr[1]) * k[0]) / (2*np.pi * 1j * k_sq)))
    v_int_y = np.real(np.fft.ifft2((np.fft.fft2(arr[0]) * k[0]) / (2*np.pi * 1j * k_sq)))

    v_int_fs = v_int_x + v_int_y
    return v_int_fs


k_int = int_2d_fourier(k, sampling)

Unfortunately, the result is not very accurate at the position where k has an abrupt change, as can be seen in the plot below, which displayes a horizontal line profile of k and k_int.

Any ideas how to improve the accuracy? Is there a way to make it exactly the same?

Answer 1

I actually found a solution. The integration itself yields very accurate results. However, the gradient function from numpy calculates second order accurate central differences, which means that the gradient itself already is an approximation.

When you replace the problem above with an analytical formula such as a 2D Gaussian, one can calculate the derivative analytically. When integrating this analytically derived function, the error is on the order of 10^-10 (depending on the width of the Gaussian, which can lead to aliasing effects).

So long story short: The integration function proposed above works as intended!