# The system matrix and the right hand side for diffusion equation with staggered grid

In the following staggered grid setting,

I want to solve diffusion equation as a linear system. $$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}$$

For discretization, use central difference for space and backward euler for time, you can obtain $$\frac{u^{n+1}-u^{n}}{\Delta t} = A u^{n+1}\\ (I-A\Delta t)u^{n+1} = u^{n}\\ \tilde{A}u^{n+1} = u^{n}$$ Where $A$ is the differential operator.

Then, solving a linear system give the velocity $u$ at the next time step $n+1$. However, to incorporate the boundary condition, I have to modify the matrix $A$ and the right hand side. the boundary conditions are given as Dirichlet condition at the boundaries, $u_N$, $u_S$, $u_W$ and $u_E$, respectively. Since the points $u_1$~$u_4$ and $u_{17}$~$u_{20}$ are located outside the domain, the velocity at those points are given as $$\frac{u_i+u_{i+4}}{2} = u_N \quad \text{and} \frac{u_i+u_{i+4}}{2} = u_S$$ Then the linear system is the following, $$\begin{pmatrix} 1 & 0 & 0 & 0 & 1 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0\\ 0 & 1 & 0 & 0 & 0 & 1 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots& 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & \cdots & \cdots & \cdots & \cdots & \cdots& 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & \cdots & \cdots & \cdots & \cdots& 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & \cdots & \cdots & \cdots & \cdots& 0\\ 0 & \frac{-1\Delta t}{dy^2} & 0 & 0 & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{2\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 0 & 0 & \frac{-1\Delta t}{dy^2} & \cdots & \cdots & 0\\ & & & & \vdots & & & & &\\ 0 & \cdots & \frac{-1\Delta t}{dy^2} & 0 & 0 & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{2\Delta t}{dy^2} & \frac{-1 \Delta t}{dx^2} & 0 & 0 & \frac{-1\Delta t}{dy^2} & \cdots & 0\\ & & & & \vdots & & & & &\\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 1 & 0 & 0 & 0 & 1\\ \end{pmatrix} \begin{pmatrix} u_1^{n+1}\\ u_2^{n+1}\\ u_3^{n+1}\\ u_4^{n+1}\\ u_5^{n+1}\\ u_6^{n+1}\\ \vdots \\ u^{n+1}\\ \vdots \\ u_{20}^{n+1}\\ \end{pmatrix} = \begin{pmatrix} 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ u_{W}\\ u_6^{n}\\ \vdots \\ u_i^{n}\\ \vdots \\ 2u_{N}\\ \end{pmatrix}$$

I want to write this system as $$\tilde{A}u^{n+1}=u^n+f$$ so that I can modify the right hand side by adding something to the velocity vector at the previous time step. Is it possible?

I understand it's possible if the system is for only the internal points $[u_6,u_7,u_{10},u_{11},u_{14},u_{15}]$, $$\begin{pmatrix} 1+\frac{2\Delta t}{dx^2}+\frac{3\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 0 & 0 & 1 & 0\\ \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{3\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 0 & 1 & 0\\ \frac{-1\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{2\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & \frac{-1\Delta t}{dy^2} & 0\\ 0 & \frac{-1\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{2\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & \frac{-1\Delta t}{dy^2}\\ 0 & 0 & \frac{-1\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{3\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2}\\ 0 & 0 & 0 & \frac{-1\Delta t}{dy^2} & \frac{-1\Delta t}{dx^2} & 1+\frac{2\Delta t}{dx^2}+\frac{3\Delta t}{dy^2}\\ \end{pmatrix} \begin{pmatrix} u_6^{n+1}\\ u_7^{n+1}\\ u_{10}^{n+1}\\ u_{11}^{n+1}\\ u_{14}^{n+1}\\ u_{15}^{n+1}\\ \end{pmatrix} = \begin{pmatrix} u_6^{n}\\ u_7^{n}\\ u_{10}^{n}\\ u_{11}^{n}\\ u_{14}^{n}\\ u_{15}^{n}\\ \end{pmatrix} + \begin{pmatrix} 2\frac{\Delta t}{dy^2}u_S+\frac{\Delta t}{dx^2}u_W\\ 2\frac{\Delta t}{dy^2}u_S+\frac{\Delta t}{dx^2}u_E\\ 0\\ 0\\ 2\frac{\Delta t}{dy^2}u_{N}+\frac{\Delta t}{dx^2}u_W\\ 2\frac{\Delta t}{dy^2}u_{N}+\frac{\Delta t}{dx^2}u_E\\ \end{pmatrix}$$

• In short, the answer is yes. I would suggest asking how this is done or what concerns you may have. Also note, the $u^{n+1}$ on the LHS is modified after moving BCs to the RHS. The way it's modified depends on whether the data is cell centered / cell corner data along the applied BC direction and the BC type (dirichlet, Neumann etc.). – Charles Nov 4 '16 at 3:05
• My motivation is to make a subspace from the velocity vectors at all the time steps so that it is convenient to have that form, the right hand side it the vector at the previous time step and the solution is the one at the next time step. In this case, all the points for u are given in the middle of the cell face as in the figure and BC are all Dirichlet. More specifically, u1-u4 and u17-u20 are given outside the boundary and u5,u8,u9,u12,u13 and u16 are give on the boundary. – user26767 Nov 4 '16 at 9:49

I found a way to do it. I can modify the RHS as following, \begin{pmatrix} 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ u_{W}\\ u_6^{n}\\ \vdots \\ u_i^{n}\\ \vdots \\ 2u_{N}\\ \end{pmatrix} into $$\begin{pmatrix} u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\\ u_{5}\\ u_6^{n}\\ \vdots \\ u_i^{n}\\ \vdots \\ u_{20}\\ \end{pmatrix} - \begin{pmatrix} u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\\ u_{5}\\ 0\\ \vdots \\ 0\\ \vdots \\ u_{20}\\ \end{pmatrix} + \begin{pmatrix} 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ 2u_{S}\\ u_{W}\\ 0\\ \vdots \\ 0\\ \vdots \\ 2u_{N}\\ \end{pmatrix}$$