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} $$
