带周期边界条件的二维插值
我正在一个二维空间中进行模拟，这个空间有周期性边界条件。一个连续的函数通过在网格上的值来表示。我需要能够在空间中的任何点上评估这个函数及其梯度。其实，这个问题本身并不难，或者说，几乎已经解决了。我们可以用scipy.interpolate.RectBivariateSpline来通过三次样条插值来处理这个函数。不过，之所以说它是几乎解决的，是因为RectBivariateSpline无法处理周期性边界条件，似乎在scipy.interpolate中没有其他方法可以做到这一点，按照我从文档中了解到的情况。
有没有可以做到这一点的Python库？如果没有，我能否调整scipy.interpolate来处理周期性边界条件？比如说，是否只需要在整个空间周围加上四个网格元素的边框，并在上面明确表示周期性条件就可以了？
[附加说明] 再多说一点，以防有用：我正在模拟动物在化学梯度中的运动。我提到的那个连续函数是它们所吸引的化学物质的浓度。这个浓度会随着时间和空间的变化而变化，遵循一个简单的反应/扩散方程。每只动物都有一个x,y坐标（不能假设在网格点上）。它们沿着吸引物的梯度移动。我使用周期性边界条件是为了简单地模拟一个无限的空间。
4 个回答
我一直在使用一个函数，这个函数可以对输入的数据进行扩展，从而创建出有效的周期性边界条件的数据。扩展数据相比于修改已有的算法有一个明显的好处：扩展后的数据可以用任何算法轻松进行插值处理。下面是一个例子。
def augment_with_periodic_bc(points, values, domain):
"""
Augment the data to create periodic boundary conditions.
Parameters
----------
points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, )
The points defining the regular grid in n dimensions.
values : array_like, shape (m1, ..., mn, ...)
The data on the regular grid in n dimensions.
domain : float or None or array_like of shape (n, )
The size of the domain along each of the n dimenions
or a uniform domain size along all dimensions if a
scalar. Using None specifies aperiodic boundary conditions.
Returns
-------
points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, )
The points defining the regular grid in n dimensions with
periodic boundary conditions.
values : array_like, shape (m1, ..., mn, ...)
The data on the regular grid in n dimensions with periodic
boundary conditions.
"""
# Validate the domain argument
n = len(points)
if np.ndim(domain) == 0:
domain = [domain] * n
if np.shape(domain) != (n,):
raise ValueError("`domain` must be a scalar or have the same "
"length as `points`")
# Pre- and append repeated points
points = [x if d is None else np.concatenate([x - d, x, x + d])
for x, d in zip(points, domain)]
# Tile the values as necessary
reps = [1 if d is None else 3 for d in domain]
values = np.tile(values, reps)
return points, values
例子
下面的例子展示了在一维情况下如何进行带有周期性边界条件的插值，但上面的函数可以应用于任意维度。
rcParams['figure.dpi'] = 144
fig, axes = plt.subplots(2, 2, True, True)
np.random.seed(0)
x = np.linspace(0, 1, 10, endpoint=False)
y = np.sin(2 * np.pi * x)
ax = axes[0, 0]
ax.plot(x, y, marker='.')
ax.set_title('Points to interpolate')
sampled = np.random.uniform(0, 1, 100)
y_sampled = interpolate.interpn([x], y, sampled, bounds_error=False)
valid = ~np.isnan(y_sampled)
ax = axes[0, 1]
ax.scatter(sampled, np.where(valid, y_sampled, 0), marker='.', c=np.where(valid, 'C0', 'C1'))
ax.set_title('interpn w/o periodic bc')
[x], y = augment_with_periodic_bc([x], y, domain=1.0)
y_sampled_bc = interpolate.interpn([x], y, sampled)
ax = axes[1, 0]
ax.scatter(sampled, y_sampled_bc, marker='.')
ax.set_title('interpn w/ periodic bc')
y_sampled_bc_cubic = interpolate.interp1d(x, y, 'cubic')(sampled)
ax = axes[1, 1]
ax.scatter(sampled, y_sampled_bc_cubic, marker='.')
ax.set_title('cubic interp1d w/ periodic bc')
fig.tight_layout()
这些功能可以在我的GitHub上找到，链接是 master/hmc/lattice.py：
- 周期性边界条件 这里的
Periodic_Lattice()类的详细信息可以在 这里 找到。 - 格点导数 在这个仓库里，你会找到一个拉普拉斯函数，一个平方梯度（如果你只想要梯度，可以开平方）以及一个重载的
np.ndarray版本。 - 单元测试 测试案例可以在同一个仓库的
tests/test_lattice.py找到。
另一个可以使用的函数是 scipy.ndimage.interpolation.map_coordinates。这个函数可以进行样条插值，并且支持周期性边界条件。虽然它不直接提供导数，但你可以通过数值方法来计算导数。
看起来，最接近我需求的Python函数是scipy.signal.cspline2d。这个函数正是我想要的，但是它假设边界条件是镜像对称的。因此，我似乎有三个选择：
自己写一个立方样条插值函数，能够处理周期性边界条件，或许可以用cspline2d的源代码作为起点（这些源代码是基于C语言写的）。
一种变通的方法：数据在位置i对位置j的样条系数的影响是按r^|i-j|来计算的，其中r = -2 + sqrt(3) ~ -0.26。所以如果我在周围加一个宽度为20的边框来复制周期性值，边缘的影响就会降到r^20 ~ 10^-5，像这样：
bzs1 = np.array( [zs1[i%n,j%n] for i in range(-20, n+20) for j in range(-20, n+20)] )
bzs1 = bzs1.reshape((n + 40, n + 40))然后我在整个数组上调用cspline2d，但只使用中间部分。这样应该可以，但看起来不太美观。
改用Hermite插值。在一个二维规则网格中，这对应于双三次插值。缺点是插值函数的二阶导数不连续。优点是（1）相对容易编码，(2) 对于我的应用来说，计算效率高。目前，我更倾向于这个解决方案。
我按照@mdurant的建议，用三角函数而不是多项式做了插值的数学计算。结果发现这和立方样条非常相似，但需要更多的计算，并且效果更差，所以我不打算这样做。
编辑：一位同事告诉我还有第四种解决方案：