# Prolongation and restriction operators in multigrid for high order PDEs

If I have the Poisson equation $$\Delta u = f$$ a standard transfer operator (for a regular grid) is the full weighting/bilinear interpolation scheme:

$$K = \frac{1}{4}\begin{bmatrix}\frac{1}{4} & \frac{1}{2} & \frac{1}{4} \\ \frac{1}{2} & 1 & \frac{1}{2} \\ \frac{1}{4} & \frac{1}{2} & \frac{1}{4}\end{bmatrix}.$$

Generally for a PDE of order $$k$$ one needs an interpolation scheme of order $$k$$ ("On the order of prolongations and restrictions in multigrid procedures", Hemker). For example the $$m$$-harmonic equation $$\Delta^m u = f$$ would require interpolation of order $$2m$$ (degree $$2m-1$$). How would I go about deriving this? As I understand it the bilinear kernel may be derived as a convolution of the 1D linear interpolation kernel $$\begin{bmatrix} \frac{1}{2} & 1 & \frac{1}{2} \end{bmatrix}$$ with itself. Does the same hold for the bi-cubic, bi-quintic etc. interpolation kernels?

For reference I have gone through Hemker's paper, Briggs' tutorial, Hackbush's book, and Trottenberg's book, but they do not go into much details there. Worse yet the provided schemes are only valid for $$H=2h$$ with the canonical fine and grid configuration. The best that I found is a remark from Trottenberg's book stating that the full weighting operator can be derived from a discrete version of the condition:

$$\int_{[x-h,x+h]\times[y-h,y+h]}w(x,y)\,dxdy = \int_{[x-h,x+h]\times[y-h,y+h]}(I^{2h}_hw)(x,y)\,dxdy$$

where the midpoint rule is used to approximate the rhs and the trapezoidal rule is used to approximate the lhs.

I found out how the kernels are typically constructed - the coefficients arise from sampling on the finer grid the interpolation synthesis function $$\phi_k$$ defined on the coarser grid. The synthesis function $$\phi_k$$ is the function used for $$k$$-degree polynomial interpolation on a regular grid with spacing $$h$$: $$f_{int}(x) = \sum_j f_j\phi_k(\frac{x}{h}-j)$$. The synthesis function can be derived by summing up the Lagrange basis polynomials on different intervals on the grid. I haven't found a simpler way of deriving it. Here are examples for linear and cubic interpolation:
$$\phi_1(x) = \begin{cases}1-|x| & |x|<1 \\ 0 & |x| \geq 1 \end{cases} \quad \phi_3(x) = \begin{cases} \frac{3}{2}|x|^3-\frac{5}{2}x^2+1 & |x|<1 \\ -\frac{1}{2}|x|^3+\frac{5}{2}x^2-4|x|+2 & 1 \leq |x| < 2 \\ 0 & |x| \geq 2\end{cases}$$
For simplicity I will consider only a downscaling of $$\approx 2$$ times. If the vector on the fine level has an odd number of elements $$N$$, a vertex-centered scheme can be used: $$p^{2h}_i = p^{h}_{2i-1}, \, i=1\ldots \lceil N/2 \rceil$$. That is, a vertex on the coarser level corresponds to the location of every other vertex on the finer level. If the vector on the fine level has an even number of elements $$N$$, a cell-centered scheme may be used: $$p^{2h}_i = \frac{1}{2}(p^h_{2i-1}+p^h_{2i}), \, i=1,\ldots, N/2$$. That is, coarse vertices are situated between two fine vertices.
Assume a $$k$$-degree polynomial interpolation is required. For a vertex-centered scheme one samples $$\phi_k$$ at $$x_k = \frac{i-1-k}{2}, \, i=1, \ldots, 2k+1$$ to get the weights in the stencil. E.g. for linear using $$\phi_1$$ it can be verified that $$\begin{bmatrix}\frac{1}{2} & 1 & \frac{1}{2}\end{bmatrix}$$ is the kernel, while for cubic using $$\phi_3$$ one gets $$\begin{bmatrix} -\frac{1}{16} & 0 & \frac{9}{16} & 1 & \frac{9}{16} & 0 & -\frac{1}{16}\end{bmatrix}$$. In the cell centered case the only difference are the sampling points which become: $$y_i = \frac{i-0.5-k}{2}, \, i=1,\ldots,2k+2$$. For linear this results in $$\begin{bmatrix} \frac{1}{4} & \frac{3}{4} & \frac{3}{4} & \frac{1}{4} \end{bmatrix}$$.
I have a detailed derivation for the interpolating functions here: Constructing $$C^{k−2} and$$C^{k−1} splines on a regular grid through convolution.