We know a function $f$ on the vertices of a regular hexagon, as follows

$$\left( 1, \ 0, \ f_{0}\right), \ \left( \frac{1}{2}, \ \frac{\sqrt{3}}{2}, \ f_{1}\right), \ \left( - \frac{1}{2}, \ \frac{\sqrt{3}}{2}, \ f_{2}\right), \ \left( -1, \ 0, \ f_{3}\right), \ \left( - \frac{1}{2}, \ - \frac{\sqrt{3}}{2}, \ f_{4}\right), \ \left( \frac{1}{2}, \ - \frac{\sqrt{3}}{2}, \ f_{5}\right)\, ,$$

and we want to interpolate a polynomial.

If we propose a quadratic polynomial of the form

$$p(x, y) = a_{0} x^{2} + a_{1} x y + a_{2} y^{2} + a_{3} x + a_{4} y + a_{5}\, ,$$

we end up with the following system

$$\left[\begin{matrix} 1 & 0 & 0 & 1 & 0 & 1\\ \frac{1}{4} & \frac{\sqrt{3}}{4} & \frac{3}{4} & \frac{1}{2} & \frac{\sqrt{3}}{2} & 1\\ \frac{1}{4} & - \frac{\sqrt{3}}{4} & \frac{3}{4} & - \frac{1}{2} & \frac{\sqrt{3}}{2} & 1\\1 & 0 & 0 & -1 & 0 & 1\\ \frac{1}{4} & \frac{\sqrt{3}}{4} & \frac{3}{4} & - \frac{1}{2} & - \frac{\sqrt{3}}{2} & 1\\ \frac{1}{4} & - \frac{\sqrt{3}}{4} & \frac{3}{4} & \frac{1}{2} & - \frac{\sqrt{3}}{2} & 1\end{matrix} \right] \left[\begin{matrix}a_{0}\\a_{1}\\a_{2}\\a_{3}\\a_{4}\\a_{5}\end{matrix}\right] = \left[\begin{matrix}f_{0}\\f_{1}\\f_{2}\\f_{3}\\f_{4}\\f_{5}\end{matrix}\right]\, ,$$

where the Vandermonde matrix is singular, and the system does not have any solution.

I suppose that this behavior has something to do with the symmetries present in the sampling points, but I am not sure, though.


  1. What is the reason for this behavior?

  2. Is there any known polynomial interpolator for the vertices of a regular hexagon?

  • 1
    $\begingroup$ For that set of points, the polynomial basis you are using is not unisolvent. Particularly, if you want to fit a quadratic polynomial using the Lagrange basis, I believe one of the nodes have to be inside the hexagon and the rest of them should be 5 out of 6 vertices of the hexagon. I don't have further intuition why that must be. IIRC, the proper basis would contain polynomials similar to the hat functions, i.e. $\phi_i(\vec{x}_j)$ should be zero if $j\neq i$ and should be one if $j=i$. I quickly searched and found the paper "On the hexagonal Shepard method" on the topic. Hope it helps. $\endgroup$ Mar 22, 2021 at 20:23
  • 1
    $\begingroup$ Why not split the hexagon into 6 triangles (with the additional vertex at the center of the hexagon), assigning the average of the 6 values to the center, and do linear interpolation on each triangle? $\endgroup$ Mar 22, 2021 at 21:04
  • $\begingroup$ @WolfgangBangerth, the idea came up trying to find the gradient of the function in the vertices of a triangular mesh. So, we are considering an interpolation over the dual mesh defining the value on each centroid. There are ways of solving the problem, but I am still curious about this case. $\endgroup$
    – nicoguaro
    Mar 22, 2021 at 21:14
  • 1
    $\begingroup$ @nicoguaro Take a look at how the "virtual element method" deals with this. They solve equations on arbitrary polygonal meshes, so a regular hexagon is a typical cell shape for them. $\endgroup$ Mar 23, 2021 at 9:01
  • $\begingroup$ Might be worthwhile instead minimizing some objective function (say a measure of the surface curvature, maybe the Gaussian curvature at the midpoint) subject to the interpolation constraints above. Lagrange multipliers, etc. $\endgroup$ Mar 23, 2021 at 22:55

5 Answers 5


As noticed by @AbdullahAliSivas the set of basis functions is not unisolvent.

In 1d case to obtain an polynomial interpolation of degree $N$ you just need of $N+1$ distinct points, $x_i \neq x_j$ if $i \neq j$. Changing the points disposition change the quality of the interpolation, because change the Lebesgue's costant. Think to Runge's phenomenon.

Going to dimension $\ge 2$ there is a drastic change. An arbitrary set of points does not guarantee the possibility of interpolation. To remark not the quality, but the possibility to have got an interpolation.

This come from a Haar's theorem. I report here:

Def Consider the linear finite dimensional subset $\mathbf{B} \subseteq C(\Omega)$ with the basis $\{B_1, \dots , B_N\}$. $\mathbf{B}$ is an Haar space on $\Omega$ if $$ \det (B_i(x_j)) \neq 0 $$ for any set of distinct $x_1, \dots, x_N$ in $\Omega$.

Note that your interpolation matrix $A$ is $A_{ij} = B_j(x_i)$.

Theorem If $\Omega \subset \mathbb{R}^d$ with $d \geq 2$ contains an interior point, the there exist no Haar spaces of continuous functions except for 1d case.

If you look the proof you can see that precisely the geometry of a space with dimension $\geq 2$ permit to swap in continuous way two points, without crash, of an arbitrary set. In in a reasoning by contradiction the determinant of matrix $A$ change sign and by continuity must go to zero at same time (here the contradiction).

This does not happens with 1d because the points are in a straight line and is not possible to swap two points without crash. To choose function $B_i$ with a dependence by the data $x_j$ allows to avoid the implication of Haar theorem and this is the door for the branch of scattered data approximation.


As requested I left a reference for the theorem:

Mairhuber, John C., On Haar’s theorem concerning Chebychev approximation problems having unique solutions, Proc. Am. Math. Soc. 7, 609-615 (1956). ZBL0070.29101.

  • $\begingroup$ Could you add a reference to this theorem? $\endgroup$ Apr 2, 2021 at 12:45
  • $\begingroup$ @AmitHochman Hi, I update the post. $\endgroup$ Apr 2, 2021 at 13:29

Since you are trying to find the gradient of the function in the vertices of a triangular mesh. So, we are considering an interpolation over the dual mesh defining the value on each centroid, you can try those area based interpolation function e.g, Wachspress Interpolation function.

The benefit of this interpolation function is that it is general enough for a polygon of any side and for n=4 it reduces to standard linear interpolation of quad element. You also do not need to place any node inside. You can construct the interpolation function just by using the coordinate of the vertices, no need of any matrix inversion to get the coefficient values. The only limitation is that the interpolation function is not in polynomial form rather a rational function.

enter image description here



Let's use the polar coordinates ($\rho,\theta$), then for the six vertices of the regular hexagon defined by $\rho$=1 and polar angles $\theta_i$ there are the following equations

$ f(\theta_i) = a_0 \cos^2(\theta_i) + a_1 \sin^2(\theta_i) + a_2 \cos(\theta_i) \sin(\theta_i) + \\ + a_3 \cos(\theta_i) + a_4 \sin(\theta_i) + a_5 = f_i $

where $f_i$ is the value at the i$^{th}$ vertex $\theta_i = 2\pi i/6,$ i $\in$ {0...5}; and $a_k$, k$\in${0...5}, are the unknown coefficients to be found. These equations are compactly written in terms of the Vandermonde matrix that is shown in the problem statement.

It looks like there are six degrees of freedom and six constraints here, so it should be a well posed problem. However, there is an identity linking the three coefficients,

$ a_0 \cos^2(\theta) + a_1 \sin^2(\theta) + a_5 = (a_0-a_1) \cos^2(\theta) + (a_1+a_5) = \\ (a_1-a_0) \sin^2(\theta) + (a_5+a_0) $

So there is some arbitrariness in the coefficients $a_i$, and that's the crux of the problem, that's what makes the original linear system of six equations ill-conditioned. The number of independent degrees of freedom is one less, whether it is a symmetric arrangement on a unit circle (the regular hexagon) or not.

For clarity of this argument, one can rewrite the original equations like this:

$ f(\theta_i) = b_0 \cos(2\theta_i) + b_1 \sin(2\theta_i) + b_2 \cos(\theta_i) + b_3 \sin(\theta_i) + b_4 = f_i, $

where $b_0=(a_0-a_1)/2$, $b_1=a_2/2$, $b_2=a_3$, $b_3=a_4$, $b_4=a_5+(a_0+a_1)/2$.

Clearly there are five independent degrees of freedom here $b_k$, k$\in${0...4}, and one should have a well posed problem setting the values of the interpolating polynomial at five arbitrary points on a unit circle.


You can treat each of the interpolation nodes as points in the complex plane, i.e., $z_k = \exp(2\pi i k/ 6)$, then form the Vandermonde matrix $V_{kl} = z_k^l$, for $k, l = 1, 2, ..., 6$, and solve $Vc = f$ for the vector of coefficients $f$. Note that in this case, the Vandermonde matrix obeys $V^HV = 6I$ and the interpolation is equivalent to interpolation with trigonometric polynomials or the discrete Fourier Transform. There is no ill-conditioning of $V$ to worry about. The interpolant is simply $\sum_k c_k z^k$. For a real-valued interpolant, just take its real part, which would be a polynomial of degree 6 in $x = \Re(z)$ and $y = \Im(z)$. Below is matlab code that implements this.

z0 = exp(1i*linspace(0, 2*pi*(1-1/6), 6));
V = vander(z0);
f = randn(6, 1);
c = V'*f/6;
x = linspace(-1.1, 1.1, 100);
[X, Y] = meshgrid(x);
Z = X + 1j*Y;
fInterp = polyval(c, Z);
mesh(x, x, real(fInterp));
hold on
plot3(real(z0), imag(z0), f, 'o')
axis equal

And here is the plot

In general, there is no guarantee you can interpolate using an arbitrary set of polynomials, as your example shows. But if you choose your basis of polynomials as the real and imaginary parts of the powers $(x+iy)^k$, for $k = 1, 2, ..., 6,$ then an interpolant is guaranteed to exist.


You need either a cubic term or 3D coordinate system to successfully interpolate in a hexagon domain. I do not know how this can be proven mathematically.

The following two papers might offer some insights.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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