# Find extension into higher dimensions for a matrix equation representing finite difference formulation of the heat equation

On page 9 of this pdf describing the finite difference formulation for the heat equation, there is a convenient tridiagonal matrix equation to represent equation 17 (which is on page 8). This represents the finite difference solution for 1D space heat conduction over time. However, the extension into 2D space would make this equation not applicable because there would be 2 more terms added for the y-dimension. I visualized this as the $(1-2r)$ term being changed to $(1-4r)$ and then surrounded on all 4 sides by $r$, so not only left and right but also into and out of the page. And then the $\phi$ term is no longer a vector but represents a 2D sheet of temperature values.

From what I can find, it seems that using tensor math would be the way to properly extend this equation into 2D, but I don't have a full grasp of it yet.

My question is, given my limited understanding of tensors:

How would I (1) organize a tensor equation that would properly collect the 5 terms in the 2D heat equation and (2) create a Matlab implementation of the tensor equation ?

To extend the discretisation to $$u_t = \alpha(u_{xx} + u_{yy}),$$ let's first leave time continuous and concentrate on the spatial terms. If $D_{xx}$ and $D_{yy}$ are matrices approximating the second derivative operator in the $x$ and $y$ directions respectively, then an appropriate discretisation is $$\mathbf{u}_t = \alpha \left[(D_{xx} \otimes I_y) + (I_x \otimes D_{yy}) \right]\mathbf{u}.$$ Here, $\otimes$ denotes the Kronecker product, defined for the $m \times n$ matrix $A$ and the $p \times q$ matrix $B$ as the $mp \times nq$ matrix $$A \otimes B = \begin{pmatrix} a_{11} B & \cdots & a_{1n} B \\ \vdots & \ddots & \vdots \\ a_{m1} B & \cdots & a_{mn} B \end{pmatrix},$$ and $I_{x,y}$ are identity matrices of appropriate sizes. If we let $i = 1, \dots, n_x$ and $j = 1, \dots, n_y$ be indices counting the grid points in the $x$ and $y$ directions respectively, then the elements in $\mathbf{u}$ are ordered in the way that would be obtained if the columns of the matrix $(u_{ij})$ were stacked on top of each other, i.e. $$\mathbf{u} = (u_{11}, \dots, u_{n_x1}, u_{12}, \dots, u_{n_x2}, \dots, u_{1n_y}, \dots, u_{n_xn_y})^\top.$$
In Matlab, the Kronecker product is performed through the command kron(). Thus, the matrix on the right-hand side of the discretisation can be obtained by typing
alpha*(kron(Dxx, eye(ny)) + kron(eye(nx), Dyy)).
Rather than having me type out a big matrix here, I encourage you to play around with this formulation on your own and get a feeling for what the resulting matrix looks like for different values of $n_x, n_y, \Delta x, \Delta y, \alpha$ and so on (start with small matrices). Once you feel comfortable with the Kronecker product and the resulting discretisation, you can discretise in time in the same way as in the pdf you linked.