积分二维矢量场阵列(反转np.gradient)
Integrate a 2D vectorfield-array (reversing np.gradient)
我有以下问题:
我想整合一个二维数组,所以基本上是反转一个梯度算子。
假设我有一个非常简单的数组如下:
shape = (60, 60)
sampling = 1
k_mesh = np.meshgrid(np.fft.fftfreq(shape[0], sampling), np.fft.fftfreq(shape[1], sampling))
然后我将我的向量场构造为一个复值数组(x 向量 = 实部,y 向量 = 虚部):
k = k_mesh[0] + 1j * k_mesh[1]
所以真实的部分看起来像这样
现在我取渐变:
k_grad = np.gradient(k, sampling)
然后我使用傅立叶变换来反转它,使用以下函数:
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)
不幸的是,结果在 k
突然变化的位置不是很准确,如下图所示,它显示了 k
和 k_int
.
有什么提高准确性的想法吗?有没有办法让它完全一样?
我居然找到了解决办法。积分本身会产生非常准确的结果。
但是numpy的梯度函数计算的是二阶准确的中心差,也就是说梯度本身已经是一个近似值了。
当您将上述问题替换为二维高斯等解析公式时,可以解析地计算导数。当对这个解析导出的函数进行积分时,误差在 10^-10 的数量级(取决于高斯的宽度,这会导致混叠效应)。
长话短说:上面提出的集成功能按预期工作!
我有以下问题: 我想整合一个二维数组,所以基本上是反转一个梯度算子。
假设我有一个非常简单的数组如下:
shape = (60, 60)
sampling = 1
k_mesh = np.meshgrid(np.fft.fftfreq(shape[0], sampling), np.fft.fftfreq(shape[1], sampling))
然后我将我的向量场构造为一个复值数组(x 向量 = 实部,y 向量 = 虚部):
k = k_mesh[0] + 1j * k_mesh[1]
所以真实的部分看起来像这样
现在我取渐变:
k_grad = np.gradient(k, sampling)
然后我使用傅立叶变换来反转它,使用以下函数:
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)
不幸的是,结果在 k
突然变化的位置不是很准确,如下图所示,它显示了 k
和 k_int
.
有什么提高准确性的想法吗?有没有办法让它完全一样?
我居然找到了解决办法。积分本身会产生非常准确的结果。 但是numpy的梯度函数计算的是二阶准确的中心差,也就是说梯度本身已经是一个近似值了。
当您将上述问题替换为二维高斯等解析公式时,可以解析地计算导数。当对这个解析导出的函数进行积分时,误差在 10^-10 的数量级(取决于高斯的宽度,这会导致混叠效应)。
长话短说:上面提出的集成功能按预期工作!