Kronecker sum of discrete Laplacians

Mathematically, using the Kronecker sum:

${\displaystyle L=\mathbf {D_{xx}} \oplus \mathbf {D_{yy}} =\mathbf {D_{xx}} \otimes \mathbf {I} +\mathbf {I} \otimes \mathbf {D_{yy}} ,\,}$

where ${\displaystyle \mathbf {D_{xx}} }$ and ${\displaystyle \mathbf {D_{yy}} }$ are 1D discrete Laplacians in the x- and y-directions, correspondingly, and ${\displaystyle \mathbf {I} }$ are the identities of appropriate sizes. Both ${\displaystyle \mathbf {D_{xx}} }$ and ${\displaystyle \mathbf {D_{yy}} }$ must correspond to the case of the homogeneous Dirichlet boundary condition at end points of the x- and y-intervals, in order to generate the 2D discrete Laplacian L corresponding to the homogeneous Dirichlet boundary condition everywhere on the boundary of the rectangular domain.

Here is a sample OCTAVE/MATLAB code to compute L on the regular 10×15 2D grid:

nx = 10; % number of grid points in the x-direction;
ny = 15; % number of grid points in the y-direction;
ex = ones(nx,1);
Dxx = spdiags([ex -2*ex ex], [-1 0 1], nx, nx); %1D discrete Laplacian in the x-direction ;
ey = ones(ny,1);
Dyy = spdiags([ey, -2*ey ey], [-1 0 1], ny, ny); %1D discrete Laplacian in the y-direction ;
L = kron(Dyy, speye(nx)) + kron(speye(ny), Dxx) ;

${\displaystyle L=\mathbf {D_{xx}} \otimes \mathbf {I} \otimes \mathbf {I} +\mathbf {I} \otimes \mathbf {D_{yy}} \otimes \mathbf {I} +\mathbf {I} \otimes \mathbf {I} \otimes \mathbf {D_{zz}} ,\,}$

where ${\displaystyle \mathbf {D_{xx}} ,\,\mathbf {D_{yy}} }$ and ${\displaystyle \mathbf {D_{zz}} }$ are 1D discrete Laplacians in every of the 3 directions, and ${\displaystyle \mathbf {I} }$ are the identities of appropriate sizes. Each 1D discrete Laplacian must correspond to the case of the homogeneous Dirichlet boundary condition, in order to generate the 3D discrete Laplacian L corresponding to the homogeneous Dirichlet boundary condition everywhere on the boundary. The eigenvalues are

${\displaystyle \lambda _{jx,jy,jz}=-{\frac {4}{h_{x}^{2}}}\sin \left({\frac {\pi j_{x}}{2(n_{x}+1)}}\right)^{2}-{\frac {4}{h_{y}^{2}}}\sin \left({\frac {\pi j_{y}}{2(n_{y}+1)}}\right)^{2}-{\frac {4}{h_{z}^{2}}}\sin \left({\frac {\pi j_{z}}{2(n_{z}+1)}}\right)^{2}}$

where ${\displaystyle j_{x}=1,\ldots ,n_{x},\,j_{y}=1,\ldots ,n_{y},\,j_{z}=1,\ldots ,n_{z},\,}$, and the corresponding eigenvectors are

${\displaystyle v_{ix,iy,iz,jx,jy,jz}={\sqrt {\frac {2}{n_{x}+1}}}\sin \left({\frac {i_{x}j_{x}\pi }{n_{x}+1}}\right){\sqrt {\frac {2}{n_{y}+1}}}\sin \left({\frac {i_{y}j_{y}\pi }{n_{y}+1}}\right){\sqrt {\frac {2}{n_{z}+1}}}\sin \left({\frac {i_{z}j_{z}\pi }{n_{z}+1}}\right)}$

where the multi-index ${\displaystyle {jx,jy,jz}}$ pairs the eigenvalues and the eigenvectors, while the multi-index ${\displaystyle {ix,iy,iz}}$ determines the location of the value of every eigenvector at the regular grid. The boundary points, where the homogeneous Dirichlet boundary condition are imposed, are just outside the grid.