[Numerical Mathemtics] How to solve hexagonal central differences

I want to simulate a 2d linear wave equation on a circle ($\displaystyle\frac{\partial^2 z(x,y,t)}{\partial t^2}=v^2\cdot\left(\displaystyle\frac{\partial^2 z(x,y,t)}{\partial x^2}+\displaystyle\frac{\partial^2 z(x,y,t)}{\partial y^2}\right)$).

To have a more significant result, I decided to use a hexagonal pattern (each point has 6 closest points at equidistant distance) as shown here: http://upload.wikimedia.org/wikipedia/en/8/81/Uniform_polyhedron-63-t0.png where the white dots are discrete checkpoints who have describe the actual value of the wave at a certain time.

To solve the problem, I want to use central differences to calculate a new situation out of the previous 2 (in time). How can I convert the central differences (that use $x,y$ values , however you rotate the situation there is maximum one dimension that fits) intro the checpoints based on the hexagonal structure?

I suppose I have to interpolate the point of a square structure out of the hexagonal points, or are there better/faster ways?

-

Write $X_{B,b} = \{\alpha \in \mathbb{Z}^B : \sum_j \alpha_j = b\}$. Now using the convention $0^0 \equiv 1$, define the matrix $W_{\alpha, \alpha'} := \alpha^{\alpha'}$. For arbitrary $f:X_{B,b} \rightarrow \mathbb{K}$ we can write $f_\alpha \equiv \sum_{\alpha'} c_{\alpha'}$ from which it follows that $c = W^{-1}f$. The facts that this procedure is well-defined, and that $W$ possesses an inverse, follow from a result in multivariate interpolation assuring us that the Lagrange interpolation problem on $X_{B,b}$ is "poised".

Actually, although generic discrete point sets admit a specific multivariate Lagrange interpolation protocol that satisfies many desirable properties, only $X_{B,b}$ does it so beautifully. As a result, we obtain a Lagrange interpolation: $f_{\mathfrak{I}}(x) := \sum_\alpha (W^{-1}f)_\alpha x^\alpha$ which satisfies $f_{\mathfrak{I}}(\alpha) = f(\alpha)$.

You can use this to define differencing schemes on a triangular (or hexagonal by suitable dual hand-waving) grid by considering $B = 3$. An example of the interpolation is shown.

Define $d_{\mathfrak{I}} f := d(f_{\mathfrak{I}})|_{X_{B,b}}$. Note (e.g.) that

$\partial_j f_{\mathfrak{I}} = \sum_\alpha (W^{-1} f)_\alpha \partial_j x^\alpha = \sum_\alpha (W^{-1} f)_\alpha \frac{\alpha_j}{x_j} \partial_j x^\alpha$

(suitably interpreted) is easy to compute in silico. Explicitly, set

$\left(W_{(\partial_j)}\right)_{\alpha, \alpha'} := \frac{\alpha'_j}{\alpha_j} \alpha^{\alpha'}, \quad \left(\mathcal{W}_{(\partial_j)}\right)_{x, \alpha'} := \frac{\alpha'_j}{x_j} x^{\alpha'}.$

Then

$\partial_j f_{\mathfrak{I}} = \mathcal{W}_{(\partial_j)} W^{-1} f, \quad \partial_j f \equiv W_{(\partial_j)} W^{-1} f.$

-
BTW, all this can be done with nice periodic (permutohedral) boundary conditions. I'll also mention that since with six neighbors as you say, you're really on the triangular lattice (though you can nevertheless adapt stuff to the hexagonal graph) it may be simpler just to tilt an axis at 60 degrees. This is often done in 2D lattice (gas or Boltzmann) simulations. –  Steve Huntsman Jul 20 '10 at 19:18