It is known that the 2D Poisson equation defined on a domain $\Omega$ (let's say $\Omega := (0,1)^2$) with Dirichlet boundary conditions $u(x,y)_{|\partial \Omega}=g(x,y)$,
$$u_{xx} + u_{yy}=f$$
can be discretized, using finite differences, to obtain a system of linear equations of the form
$$A\vec{u} = \vec{f}+\vec{b}$$
where $A$ is a coefficient matrix, $\vec{u}$ the discrete solution to solve for, and $\vec{b}$ contains terms from the boundary conditions.
The matrix $A$ for the 2D Poisson equation has the block form
\begin{bmatrix}
D & I & \dots &0 \\
I & \ddots & \ddots & \vdots\\
\vdots & \ddots & \ddots & I\\
0 & \dots & I & D
\end{bmatrix}
where
$$D=\begin{bmatrix}
-4 & 1 & 1 &0 & \dots& 0 \\
1 & -4 & 1 & 1 & \dots & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\
\vdots & \ddots & \vdots & \dots & -4 & 1\\
0 & \dots & \dots & 1 & 1 & -4
\end{bmatrix}$$
I want to discretize the 3D Poisson equation
$$u_{xx} + u_{yy}+u_{zz}=f$$
on $\Omega := (0,1)^3$.
So, I suppose, that the discretized set $\Omega_{\Delta x}$ must be a 3D mesh. The resulting discretization equations will be
$$\Delta_h u_{i,j,k} = \frac{u_{i+1,j,k}+u_{i,j+1,k}+u_{i,j,k+1}-6u_{i,j,k}+u_{i-1,j,k}+u_{i,j-1,k}+u_{i,j,k-1}}{(\Delta x)^2}=f_{i,j,k}+b_{i,j,k}$$
So that $A$ should now look like this:
$A=\begin{bmatrix}
-6&1 & 1 & 1 &0 & 0 & 0&\dots& 0 \\
1 & 1 & -6 & 1 & 1 & 1 &0&\dots & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\vdots&\ddots&0\\
\vdots & \vdots & \vdots & \vdots & \dots & \dots & \dots&-6 & 1\\
0 & \dots & 0 & 0 & 0 & 1 & 1 & 1 & -6
\end{bmatrix}$
I'm not quite sure, however, what the form of $D$ exactly is. I.e. is $A$ still block-tridiagonal, with two block diagonals comprised of $I$'s, and the main diagonal comprised of $D$'s, but with higher dimensions?
Would appreciate some insight.
Best Answer
I like a lot the Kronecker product representation. If $A_1$ is the one-dimensional second-derivative approximation, i.e. $$ A_1 = \begin{pmatrix} 2 & -1 & \\ -1 & 2 & \ddots \\ & \ddots & \ddots \end{pmatrix}, $$ then you can write the $3d$ Laplacian as $$ A_3 = I \otimes I \otimes A_1 + I \otimes A_1 \otimes I + A_1 \otimes I \otimes I. $$ For the $2d$ case it is the same, i.e., $A_2 = I \otimes A_1 + A_1 \otimes I$, and this generalizes to higher dimensions. You can have a look into a small instance of this matrix, and its properties (e.g., sparsity pattern) will be clear.