Why is it that when we move from 1D to 2D or 3D, we can represent the whole system in terms of Kronecker products between identity matrices and the coefficient matrix generated for the 1D case? Is it just mathematical coincidence or is there a physical meaning to this?

In the question above, this representation was carried out for Dirichlet conditions.

2 Answers
2

Maybe this isn't a helpful response but the reason this happens for the matrix form of Laplacians is because this actually happens for the true infinite-dimensional Laplacians in some settings. In recangular coordinates, $N$-D Laplacians act exactly like their discretized counterparts in this way. For example, let $\Delta_3$ by the Laplacian on $[0,1]^3$ acting on the space $H^1([0,1]^3)$. One can view this function space as a tensor product of 3 1D space (iirc there is some subtlety here but I believe it works):
$$
H^1([0,1]^3) = \bigotimes_{i=1}^3H^1([0,1]).
$$
You can view each of these spaces as consisting of functions in a different coordinate, i.e., one space is functions of $x$, one is functions of $y$, and one of $z$. With this in mind, we can represent $\Delta_3$ as a combination of different 1D Laplacians, $\Delta_1$. Let $id$ be the identity operator and recall that $\Delta_1$ is just the second derivative, then we have
$$
\Delta_3 = \Delta_1\otimes id\otimes id+id\otimes\Delta_1\otimes id+id\otimes id\otimes \Delta_1 =\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}.
$$
Notice that this is basically identical to Kronecker product representations of discretized Laplacians. Recalling our tensor representation of the function space, we can view the tensor product form of $\Delta_3$ as "Apply $\Delta_1$ to the function of $x$ and leave the others alone, then add $\Delta_1$ applied to the function of $y$ and leave the others alone, etc."

I believe that the key fact allowing us to do this is the fact that both the domain $[0,1]^n$ and the related function spaces are easy to decompose into a tensor product of 1D spaces and that the Laplacian is isotropic, i.e., it acts the same in every direction. It is possible to decompose operators that are not isotropic, i.e., $x\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}$ but then the decomposition depends on the coordinates, whereas the Laplacian is the same for any set of orthonormal coordinates.

This is my attempt at providing some intuition. Everything I state might be obvious, moreover it doesn't have much to do with physics, so this could be a non-answer. I will ignore boundary conditions.

Imagine we have a 2D grid with function values $u_{ij}$ for $i=1\ldots m, j=1\ldots n$. Let the $x$ axis point downward, the $y$ axis to the right. N.b. $z$ would point towards us.

y
---->
|
|
v
x

We will pack these $u_{ij}$ into an $m\times n$ matrix, or an $nm$ sized (column) vector. Both of these representations live in finite dimensional vector spaces. A basis for the first one is
$
e_i \otimes e_j^T,
$
one for the second one is $e_i \otimes e_j$. Here $e_i, e_j$ are the usual basis (column) vectors.

For example, assuming $n=m=2$, $e_2 \otimes e_2^T$ would be

0 0
0 1,

while $e_2 \otimes e_2$ is the transpose of

0 0 0 1.

These spaces are isomorphic, and we have an appropriate pair of bases. While the second representation is more commonly used in applications, the first one is useful in understanding how the discretized Laplacian acts.

We should take a moment to note that given appropriate matrices $A$ and $B$, we can define $A \otimes B$ by
$$
(A \otimes B)(e_i \otimes e_j) = Ae_i \otimes Be_j.
$$

A word about the 3D case. The idea is the same. Just as the definition of $\partial_{11}$ leaves the last $2$ variables fixed, the discrete $D^3_{11}$ will not mix together elements along the $y$ and $z$ axes, which makes sense of the results
$$
D^3_{11} = D^1_{11} \otimes I \otimes I,
$$
and
$$
D^3_{11} = D^1_{11} \otimes I \otimes I + I \otimes D^1_{11} \otimes I + I \otimes I \otimes D^1_{11}.
$$

As a final remark, arranging the function values in accordance with the usual $(x,y)$ plane, the first two terms have to be flipped, e.g.
$$
D^{2}_{11} = I \otimes D^1_{11}, \quad
D^{2}_{22} = D^1_{11} \otimes I.
$$