[Numerical Mathemtics] How to solve hexagonal central differences

Question

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?

Answer 1

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.$