Implicit 2D finite difference linear system

Question

I am familiar with the 1D implicit method which solves the heat equation with homogeneous Dirichlet conditions, $$u_i^{k-1} = \big(1+2\lambda)u_i^k - \lambda \big( u_{i+1}^k + u_{i-1}^k \big) $$ which can be rearranged in to a linear system $$\underline{u} ^{k-1} = A \underline{u} ^{k}$$ where $$A = \begin{pmatrix} 1+2\lambda & -\lambda & 0 & \dots & 0\\ -\lambda & 1+2\lambda & -\lambda & \ddots & \vdots\\ 0 & -\lambda & \ddots & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\lambda\\ 0 & \dots & 0 & -\lambda & 1+2\lambda\\ \end{pmatrix}.$$ However, I am becoming increasingly confused about what the linear system for the 2 Dimensional case would look like. I am aware that the so called algorithm would now be, $$ u_{i,j}^{k-1} = \big(1+4\lambda)u_{i,j}^k - \lambda \big( u_{i+1,j}^k + u_{i-1,j}^k + u_{i,j+1}^k + u_{i,j-1}^k \big). $$ $\textbf{My question - What is the linear system for this algorithm?}$

I know that it can have the same format as the 1D case but I'm unsure what the vector $\underline{u}$ would be since there are now two different spatial dimensions. (i.e we have both $i$ and $j$ varying). I am also unsure of what the coefficient matrix would be. Thanks in advanced for any help.

Answer 1

Warning : this answer may be overkill for your purposes and have an intimidating abstraction level where we would do well to be equipped with (at least) the background knowledge equivalent to a well digested Honors linear-algebra course.

As long as you use a lexical ordering you can be systematic and use Kronecker products with identity matrices in expressing linearly separable equations.

Scalar products can calculate the index of a lexical order. If you excuse the perhaps confusing notation. For two dimensions of your named indices $i$ and $j$ :

$(i,j) \in \{0,\cdots,N_i-1\}\times\{0,\cdots,N_j-1\} \to [i,j] \cdot [1,N_i]^T$

Now in the general case of dimensionality $\geq 3$ the vector to the right would have element 3 = $N_i\cdot N_j$ systematically multiplying up the number of elements of previously passed dimensions. In matlab this would be a "cumprod" operation on the size vector padded with a $1$ at the start.

The nice thing with being systematic like this is that we easily can switch between (i,j) and vectorized index.

Now separable linear operators for such lexical orders are nice to express with Kronecker products.

for operator in the 1 dimensional case represented by matrix $M$:

along i dimension :

$M\otimes I_{N_j}$

along j dimension :

$I_{N_i}\otimes M$

The nice thing here is that with increasing dimensions we can just systematically stack identity matrices and operator matrices as factors where they belong in the Kronecker product.

Answer 2

Basically you have just to make one decision: How do you want to order the vector $\boldsymbol{u}$. That can be done in row-major ordering, i.e., $$\boldsymbol{u} = \begin{pmatrix} u_{1, 1}& u_{2, 1} & \dots & u_{N_x, 1} & u_{1, 2} & u_{2, 2} & \dots & u_{N_x, N_y} \end{pmatrix}^T $$ or in column-major $$\boldsymbol{u} = \begin{pmatrix} u_{1, 1}& u_{1, 2} & \dots & u_{1, N_y} & u_{2, 1} & u_{2, 1} & \dots & u_{N_x, N_y} \end{pmatrix}^T $$ The qualitative structure of the matrix is then the same. You have two more diagonals (one below and one above the main diagonal). The "distance" from the diagonal element to the "far" off diagonals is then determined by your ordering. For row-major, the distance in $\boldsymbol{u}$ from $u_{i,j} $ to $u_{i, j \pm 1}$ is $\pm N_x$, while for column-major, the distance from $u_{i,j} $ to $u_{i \pm 1, j }$ is $\pm N_y $. Thus, the matrix A reads: $$A = \begin{pmatrix} 1+2\lambda & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 0 & \dots & 0 \\ -\lambda & 1+2\lambda & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} -\lambda & 0 & \dots & 0\\ 0 & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots &\\ -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 1+2\lambda & -\lambda \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 0 \\ 0 & & \dots & -\lambda & \underbrace{0}_{\{N_x, N_y\} - 2} & -\lambda & 1+2\lambda \end{pmatrix}.$$