Kronecker Sum of Discrete Laplacians
We consider a two-dimensional grid of \(\,m\times n\,\) points
with a spacing \(\,h.\ \) \(\\\)
Application of the finite difference method with the five-point stencil
\[\begin{split}\begin{array}{ccc}
& (x-h,y) &
\\
& | &
\\
(x,y-h) & -\quad (x,y)\quad - & (x,y+h)
\\
& | &
\\
& (x+h,y) &
\end{array}\end{split}\]
leads to the approximate expression for the Laplacian of a function
\(\,f(x,y)\,\) defined on the grid:
\[\Delta f(x,y)\ \approx\
\frac{1}{h^2}\ \
[\,f(x-h,y)+f(x,y-h)+f(x,y+h)+f(x+h,y)-4\,f(x,y)\,].\]
We shall put henceforth \(\,h=1.\ \)
Thus \(\,x=1,2,\ldots,m,\ y=1,2,\ldots,n,\ \) and
\[\Delta f(x,y)\ =\
\Delta^{(xx)} f(x,y)\ +\ \Delta^{(yy)} f(x,y)\]
where the 2nd order partial derivatives correspond to the three-point
stencil in one dimension:
(1)\[ \begin{align}\begin{aligned}\Delta^{(xx)} f(x,y)\ \equiv\
\frac{\partial^{\,2}}{\partial x^2}\ f(x,y) \approx\
f(x-1,y)\ - \ 2\,f(x,y)\ + \ f(x+1,y),\\\Delta^{(yy)} f(x,y)\ \equiv\
\frac{\partial^{\,2}}{\partial y^2}\ f(x,y) \approx\
f(x,y-1)\ - \ 2\,f(x,y)\ +\ f(x,y+1).\end{aligned}\end{align} \]
To write these relations in a matrix form, we use the matrix notation:
\[\begin{split}\begin{array}{rl}
f(x,y) \rightarrow f_{xy}\,, &
\boldsymbol{f}\ :\,=\ [\,f_{xy}\,]_{\,m \times n}\,;
\\
\Delta f(x,y) \rightarrow (\Delta f)_{xy}\,, &
\boldsymbol{\Delta} f\ :\,=\ [\,(\Delta f)_{xy}\,]_{\,m \times n}\,;
\\ \\
\Delta^{(xx)} f(x,y) \rightarrow (\Delta^{(xx)} f)_{xy}\,, &
\boldsymbol{\Delta}^{(xx)} f\ :\,=\ [\,\Delta^{(xx)} f)_{xy}\,]_{\,m \times n}\,;
\\
\Delta^{(yy)} f(x,y) \rightarrow (\Delta^{(yy)} f)_{xy}\,, &
\boldsymbol{\Delta}^{(yy)} f\ :\,=\ [\,\Delta^{(yy)} f)_{xy}\,]_{\,m \times n}\,;
\\ \\
\boldsymbol{D}^{(xx)}\ =\ [\,d^{(xx)}_{vw}\,]_{m\times m}\,, &
\boldsymbol{D}^{(yy)}\ =\ [\,d^{(yy)}_{vw}\,]_{n\times n}\,.
\end{array}\end{split}\]
The index \(\,y\in\{1,2,\ldots,n\}\,\) being fixed,
the first of Eqs. (1) can be rewritten as
\[ \begin{align}\begin{aligned}(\Delta^{(xx)} f)_{xy}\ \ =\ \
\displaystyle\sum_{w=1}^m\ d^{(xx)}_{xw}\ f_{wy}\,,
\qquad \forall\ x = 1,2,\ldots,m;\\\boldsymbol{\Delta}^{(xx)} f\ \ =\ \
\boldsymbol{D}^{(xx)}\cdot\,\boldsymbol{f}\,,\\\boldsymbol{\Lambda}^{mn}(\boldsymbol{\Delta}^{(xx)}f)\ =\
(\boldsymbol{D}^{(xx)}\otimes\boldsymbol{I}_n)\ \cdot\
\boldsymbol{\Lambda}^{mn}(\boldsymbol{f})\,.\end{aligned}\end{align} \]
Here \(\,\boldsymbol{\Lambda}^{mn}(\boldsymbol{A})\,\) is the column
of coordinates of a \(\,m\times n\,\) matrix \(\,\boldsymbol{A}\,\)
in the standard basis of the vector space \(\,M_{m\times n}(K).\ \)
The last formula above is the vectorized form of the preceding matrix equation
(see the section ‘Vectorization of Matrices’ in this chapter).
On the other hand, for a fixed \(\,x\in\{1,2,\ldots,m\},\ \)
the second equation in (1) yields
\[ \begin{align}\begin{aligned}(\Delta^{(yy)} f)_{xy}\ \ =\ \
\displaystyle\sum_{v=1}^n\ d^{(yy)}_{yv}\ f_{xv}\ \ =\ \
\displaystyle\sum_{v=1}^n\ f_{xv}\ \left[\,d^{(yy)}\,\right]_{vy}^{\,T}\,,
\qquad \forall\ y = 1,2,\ldots,n;\\\boldsymbol{\Delta}^{(yy)} f\ \ =\ \
\boldsymbol{f}\,\cdot\,\left[\,\boldsymbol{D}^{(yy)}\,\right]^{\,T}\,,\\\boldsymbol{\Lambda}^{mn}(\boldsymbol{\Delta}^{(yy)}f)\ =\
(\boldsymbol{I}_m\otimes\boldsymbol{D}^{(yy)})\ \cdot\
\boldsymbol{\Lambda}^{mn}(\boldsymbol{f})\end{aligned}\end{align} \]
(once again we used vectorization formulas for a matrix product
from a preceding section).
Gathering together the partial results for \(\ \,\boldsymbol{\Delta}f\,=\,
\boldsymbol{\Delta}^{(xx)}f\ + \ \boldsymbol{\Delta}^{(yy)}f\,,\ \ \)
we obtain
\[\blacktriangleright\quad
\boldsymbol{\Lambda}^{mn}(\boldsymbol{\Delta}f)\ =\
\left(\boldsymbol{D}^{(xx)}\oplus\,\boldsymbol{D}^{(yy)}\right)\ \cdot\
\boldsymbol{\Lambda}^{mn}(\boldsymbol{f})\,,\]
where the \(\,\) Kronecker sum \(\,\) of matrices
\(\,\boldsymbol{D}^{(xx)}\in M_m(K)\ \) and
\(\,\boldsymbol{D}^{(yy)}\in M_n(K)\ \) is defined as
\[\triangleright\quad
\boldsymbol{D}^{(xx)}\oplus\,\boldsymbol{D}^{(yy)}\ :\,=\ \,
\boldsymbol{D}^{(xx)}\otimes\,\boldsymbol{I}_n\ +\
\boldsymbol{I}_m\otimes\,\boldsymbol{D}^{(yy)}\,.\]
To get some idea how the matrices \(\,\boldsymbol{D}^{(xx)}\ \)
and \(\,\boldsymbol{D}^{(yy)}\ \) look, we shall use the tools of Sage.
The function ML(n) (‘Minus Laplacian’) returns the matrix
\(\,\boldsymbol{D}^{(xx)}\ \) or \(\,\boldsymbol{D}^{(yy)}\ \)
of a given size \(\,n:\)
For example, for a \(\,3\times 5\ \) grid of 15 points \(\,\)
(\(m=3,\ n=5\)) \(\,\) the matrices read:
\[\begin{split}\boldsymbol{D}^{(xx)}\ =\
\left[\begin{array}{rrr}
-1 & 1 & 0 \\
1 & -2 & 1 \\
0 & 1 & -1
\end{array}\right]
\,,\qquad
\boldsymbol{D}^{(yy)}\ =\
\left[\begin{array}{rrrrr}
-1 & 1 & 0 & 0 & 0 \\
1 & -2 & 1 & 0 & 0 \\
0 & 1 & -2 & 1 & 0 \\
0 & 0 & 1 & -2 & 1 \\
0 & 0 & 0 & 1 & -1
\end{array}\right]\,.\end{split}\]
Unlike the internal rows, the first and the last rows of these matrices
do not result from (1).
Instead, they have been chosen so that \(\,\boldsymbol{D}^{(xx)}\ \)
or \(\,\boldsymbol{D}^{(yy)}\ \) are the negative Laplacian matrix
of a one-dimensional grid with \(\,m\,\) or \(\,n\,\) points,
respectively.
The code below shows the matrix
\(\,\boldsymbol{D}^{(xx)}\oplus\,\boldsymbol{D}^{(yy)}\ \)
for a \(\,3\times 5\ \) grid of 15 points \(\,\)
(\(m=3,n=5\)):
