Tensor product representation for the 9-point finite difference approximations for the Poisson equation

Question

If we use 5-point finite difference approximations in a uniform rectangular grid to solve the Poisson PDE

\begin{align} -\Delta u &= f \ \ \text{en} \ \ (0,1)\times (0,1) \label{P1} \\ u &= 0 \ \ \text{en} \ \partial ((0,1)\times (0,1)) \label{P2} \end{align}

The Laplacian operator $\nabla^2$ can be written using the Kronecker product as follows: $$ \nabla^2 = T\otimes I_{n_{x}} + I_{n_{y}} \otimes T $$

where $I_{n}$ is the identity matrix of size $n$ and the matrix $T$ is given by \begin{equation} T = \frac{1}{h^2}\left[\begin{matrix} -2 & 1 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & -2 \end{matrix} \right] \end{equation}

where $h>0$ is the mesh size.

EDIT: The 9-point finite difference approximation has the following stencil:

If we consider the natural rowwise order

see the book by Randall J. LeVeque - Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems.

It can shown that the 9-point approximations leads to

\begin{equation} \nabla^2 = \frac{1}{6h^2}\left[\begin{matrix} D_{1} & D_{2} & & 0\\ D_{2} & \ddots & \ddots & \\ & \ddots & \ddots & D_{2} \\ 0 & & D_{2} & D_{1} \end{matrix} \right] \end{equation}

where

\begin{equation} D_{1} = \left[\begin{matrix} -20 & 4 & & 0\\ 4 & \ddots & \ddots & \\ & \ddots & \ddots & 4 \\ 0 & & 4 & -20 \end{matrix} \right] \end{equation}

and

\begin{equation} D_{2} = \left[\begin{matrix} 4 & 1 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & 4 \end{matrix} \right] \end{equation}

Now consider the matrices

\begin{equation} K = \frac{1}{6h^2}\left[\begin{matrix} 10 & 1 & & 0\\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ 0 & & 1 & 10 \end{matrix} \right] \end{equation}

\begin{equation} I_{R} = \left[\begin{matrix} -2 & 2/5 & & 0\\ 2/5 & \ddots & \ddots & \\ & \ddots & \ddots & 2/5 \\ 0 & & 2/5 & -2 \end{matrix} \right] \end{equation}

\begin{equation} I_{L} = \left[\begin{matrix} 0 & 3/5 & & 0\\ 3/5 & \ddots & \ddots & \\ & \ddots & \ddots & 3/5 \\ 0 & & 3/5 & 0 \end{matrix} \right] \end{equation}

A quick computation shows that

\begin{equation} K \otimes I_{R} + I_{L} \otimes K = \frac{1}{6h^2}\left[\begin{matrix} D_{1} & D_{2} & & 0\\ D_{2} & \ddots & \ddots & \\ & \ddots & \ddots & D_{2} \\ 0 & & D_{2} & D_{1} \end{matrix} \right] \end{equation}

Questions:

1. Do you know tensor product representations for the 9-point finite difference approximations for the $\nabla^2$ operator?
2. Do you know books or articles where I can study tensor product representations for the 9-point or higher order finite difference approximations for the $\nabla^2$ operator?
3. Do you know books or articles where I can study implementation details for the Poisson equation in 2D?

New questions:

1. Is this a valid tensor representation for the 9-point Laplacian operator or do you know other representation for this operator?
2. Is there a general way to deduce a tensor representation for high-order finite difference approximations for the Laplacian operator?

Thanks in advance.

Answer 1

This is too long for a comment, so I'll post an answer.

If you start from the analytical 2D-Laplace operator, it naturally is already in a (sum of) tensor product form:

$$ \Delta = \partial^2_x \otimes I + I \otimes \partial^2_y\,. $$

By looking at the operator, it seems obvious how to discretize it, namely by using some one-dimensional finite difference approximation for each of the derivatives. This will lead exactly to the desired form -- you seem to know that as well, as you're refering to the 5-point stencil (which is the special case when you pick the 1,-2,1 discretization).

Now, for some reason, you want to go beyond the cross-like stencil and incorporate further points. For this, you can start again from the partial derivatives $\partial^2_x$ and $\partial^2_y$ discretize those separately on, say, a nine-point stencil and sum the result. As a result, you will get a formula like the last one posted in the OP quite naturally -- i.e. without having to apply linear algebra to disentangle the dimensions.

Two further pointers:

1. As said in the comments, it's appropriate to work with polynomials, and on a product grid, that is particularly easy as one can use products of one-dimensional polynomials (e.g. Lagrange polynomials). One can then apply some algorithm to find the 1D-derivatives (e.g. Fornberg's) and simply multiply them together to obtain the 2D-derivative.

2. If you're interested in general decompositions of the form you stated in the last formula of the original question, have a look at the "higher-order singular value decomposition", which (approximately) brings a tensor to the sum of product form you mentioned.