MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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: 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?

share|cite|improve this question
up vote 2 down vote accepted

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


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

share|cite|improve this answer
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.