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

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?

• For question 1, it boils down to writing your matrix T for 7-point or 9-point FD, right? Jun 25, 2021 at 16:06
• If you dont want to waste your time with different FD approximation stencils start to think in polynomials. Once you got it, you can derive the operator ins secondes. Jun 25, 2021 at 20:39
• I don't think that 7- and 9-point stencil matrices can be written as tensor products. I would run a loop similar to scicomp.stackexchange.com/questions/37502/… . You can also check devitoproject.org Jun 26, 2021 at 3:16
• I am taking my comment back. Chen Long has them all in his notes: math.uci.edu/~chenlong/226/FDMcode.pdf Jun 26, 2021 at 3:22
• It appears that you answered your question in the same post and added new questions. I suggest that you add an answer and create a new post if you have new questions. Jun 30, 2021 at 13:12

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.