Discretization of a heat equation using finite-difference method

Let $\overline{\Omega}=\left([0,2]\times [0,1]\right)\cup \left([1,2]\times [1,2]\right)$ and $\Gamma=\Gamma_1\cup\ldots\cup\Gamma_6$ be defined as shown in the following picture:

I want to discretize the following heat equation: $$\begin{split} -D\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) &= 0, \text{ for all }(x,y)\in\Omega\\ D\left(\frac{\partial u}{\partial x}n_1+\frac{\partial u}{\partial y}n_2\right) &= 0, \text{ for all }(x,y)\in\Gamma_2\cup\Gamma_3\cup\Gamma_5\cup\Gamma_6\\ u &= u_0, \text{ in }\Gamma_1\\ u &= u_1, \text{ in }\Gamma_4 \end{split}$$ where $h$ is the step width of the used lattice and $D,u_0,u_1$ are given constants. I want to find a matrix $K_h$ and a vector $b_h$ such that $K_hu_h=b_h$, where $u_h$ is the approximate solution vector.

If $(x_i,y_j)$ is a point on the lattice and $u_{ij}=u(x_i,y_j)$, I want to use the approximations $$\Delta u_{ij}=\nabla^2u_{ij}\approx\frac{1}{h^2}\left(u_{i+1,j}+u_{i-1,j}-4u_{ij}+u_{i,j+1}+u_{i,j-1}\right)$$ and $$\frac{\partial}{\partial x}u_{ij}\approx\frac{u_{i+1,j}-u_{i-1,j}}{2h}$$ as well as $$\frac{\partial}{\partial y}u_{ij}\approx\frac{u_{i,j+1}-u_{i,j-1}}{2h}$$

Suppose we first consider $\Omega_1=(0,1)\times (0,1)$. We construct our grid $G^{(h)}$ by the tensor product of two grids of $(0,1)$: $$G_x^{(h)}:=\left\{x_i=(i-1)h:1\le i\le N\right\}\text{ and }G_y^{(h)}:=\left\{y_j=(i-1)h:1\le j\le N\right\}$$ where N:=\frac{h+1}{h}. Using our approximations, we've got the following equations:

$$\begin{split} D\frac{1}{h^2}\left(4u_{ij}-u_{i+1,j}-u_{i-1,j}-u_{i,j+1}-u_{i,j-1}\right)&=0\text{ , for all }1\le i,j\le N\\ -D\frac{u_{i,2}-u_{i,0}}{2h}&=0\text{ , for all } 0\le i\le N+1\\ D\frac{u_{i,N+1}-u_{i,N-1}}{2h}&=0\text{ , for all } 0\le i\le N+1\\ u_{1j}&=u_1\text{ , for all } 1\le j\le N \end{split}$$

Please note that I've split the second equations for the half of $\Gamma_2$ and $\Gamma_6$. There some problems left: First of all, $u_{i,0}$ and $u_{i,N+1}$ are not well-defined. We can use the first equation for $j=1$ and $j=N$ in combination with the second and third equation, respectively. Insertion and scaling by a factor of $1/2$ yields $$D\frac{1}{h^2}\left(2u_{i,1}-\frac12 u_{i+1,1}-\frac12 u_{i-1,1}-u_{i,2}\right)=0$$ and $$D\frac{1}{h^2}\left(2u_{iN}-\frac12 u_{i+1,N}-\frac12 u_{i-1,N}-u_{i,N-1}\right)=0$$ Now there are still problems with the corner nodes $(0,0)$, $(0,1)$ and (maybe?) $(1,1)$.

For $(0,0)$ we consider $(i,j)=(0,0)$ and there are 3 equations involved: From the first equation we've got $$4u_{1,1}-u_{2,1}-u_{0,1}-u_{1,2}-u_{1,0}=0\tag{1}$$ From the second equation we've got: $$u_{1,2}-u_{1,0}=0\tag{2}$$ and the third equations yields $$u_{1,1}=u_1\tag{3}$$ However, we've got now information about $u_{0,1}$.

So there are three questions left:

1. What can we say about $u_{0,1}$?
2. How do we describe all these equations on $\Omega_1$ in a linear equation system "$Ax=b$"
3. Suppose we've got all linear three linear equation systems for $\Omega_1$ as well as $\Omega_2=(1,2)\times (0,1)$ and $\Omega_3=(1,2)\times (1,2)$. How do we combine these systems into one system for $\Omega$?

Make a rectangular lattice on $[0,1]^2$ and on $[1,2] \times [0,2]$ separately. Merge them (in the process making the lattice points on $\{ 1 \} \times [0,1]$ only count once, of course). The awkward thing in this circumstance is that the "ghost points" have to be defined in a rather asymmetric fashion.
• @oxbadfood That depends on the data structure which is convenient for your application. If you're just writing down a matrix equation, then you'll need to reindex. One convenient way to do this might be to actually split the problem into $[0,1] \times [0,1]$, $[1,2] \times [0,1]$, $[1,2] \times [1,2]$. Then you have three square lattices with exactly the same number of mesh points (modulo the ghost point problem), which makes designing a reindexing scheme not very difficult. – Ian Oct 26 '14 at 21:58
• @oxbadfood A lazy approach would be to extend your lattice to the whole square $[0,2]^2$ (or actually, $[-h,2+h]^2$, because of the Neumann conditions on the left, bottom, and right). Then you write your equations for the interior by going across rows of the lattice, with equations for the exterior nodes being simply $u_{i,j} = 0$. Just be careful not to do this for the points which are not "truly" outside the domain, i.e. your Neumann BC ghost points. Also be sure to use a sparse matrix structure if you do this, because otherwise it will be more than twice as slow as it has to be. – Ian Oct 27 '14 at 11:32
• @oxbadfood You can deal with the boundary by using ghost points. You extend the DE to the boundary of your domain proper, so at for example $u_{i,0}$ you have $u_{i,0} = \text{average}(u_{i,-1},u_{i,1},u_{i+1,0},u_{i-1,0})$. Then you formally insert points such as this $u_{i,-1}$ to use for the Neumann boundary condition. That is, the normal derivative at $u_{i,0}$ is proportional to $u_{i,1}-u_{i,-1}$, and similar elsewhere. I forget what exactly you have to do in order to manage the corners (or any other point where analytically speaking the outward normal is not defined). – Ian Oct 27 '14 at 17:28
• As for writing everything in a linear system, you have to perform some sort of interlacing in order to turn your 2D array $\{ u_{i,j} \}_{i=1,j=1}^{i=N_x,j=N_y}$ into a 1D array. Then you just write down each equation after reindexing using this interlacing scheme. – Ian Oct 27 '14 at 17:33