2D Interpolation with periodic boundary conditions

Question

I'm running a simulation on a 2D space with periodic boundary conditions. A continuous function is represented by its values on a grid. I need to be able to evaluate the function and its gradient at any point in the space. Fundamentally, this isn't a hard problem -- or to be precise, it's an almost already solved problem. The function can be interpolated using a cubic spline with scipy.interpolate.RectBivariateSpline. The reason it's almost solved is that RectBivariateSpline cannot handle periodic boundary conditions, nor can anything else in scipy.interpolate, as far as I can figure out from the documentation.

Is there a python package that can do this? If not, can I adapt scipy.interpolate to handle periodic boundary conditions? For instance, would it be enough to put a border of, say, four grid elements around the entire space and explicitly represent the periodic condition on it?

[ADDENDUM] A little more detail, in case it matters: I am simulating the motion of animals in a chemical gradient. The continuous function I mentioned above is the concentration of a chemical that they are attracted to. It changes with time and space according to a straightforward reaction/diffusion equation. Each animal has an x,y position (which cannot be assumed to be at a grid point). They move up the gradient of attractant. I'm using periodic boundary conditions as a simple way of imitating an unbounded space.

Answer 1

It appears that the python function that comes closest is scipy.signal.cspline2d. This is exactly what I want, except that it assumes mirror-symmetric boundary conditions. Thus, it appears that I have three options:

1. Write my own cubic spline interpolation function that works with periodic boundary conditions, perhaps using the cspline2d sources (which are based on functions written in C) as a starting point.

2. The kludge: the effect of data at i on the spline coefficient at j goes as r^|i-j|, with r = -2 + sqrt(3) ~ -0.26. So the effect of the edge is down to r^20 ~ 10^-5 if I nest the grid within a border of width 20 all the way around that replicates the periodic values, something like this:

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

   Then I call cspline2d on the whole array, but use only the middle. This should work, but it's ugly.

3. Use Hermite interpolation instead. In a 2D regular grid, this corresponds to bicubic interpolation. The disadvantage is that the interpolated function has a discontinuous second derivative. The advantages are it is (1) relatively easy to code, and (2) for my application, computationally efficient. At the moment, this is the solution I'm favoring.

I did the math for interpolation with trig functions rather than polynomials, as @mdurant suggested. It turns out to be very similar to the cubic spline, but requires more computation and produces worse results, so I won't be doing that.

EDIT: A colleague told me of a fourth solution:

1. The GNU Scientific Library (GSL) has interpolation functions that can handle periodic boundary conditions. There are two (at least) python interfaces to GSL: PyGSL and CythonGSL. Unfortunately, GSL interpolation seems to be restricted to one dimension, so it's not a lot of use to me, but there's lots of good stuff in GSL.

Answer 2

Another function that could work is scipy.ndimage.interpolation.map_coordinates. It does spline interpolation with periodic boundary conditions. It does not not directly provide derivatives, but you could calculate them numerically.