I'm tried to solve numerically the following PDE (forget that an exact solution is known):
\begin{array}{l l} \partial_t \rho - \varepsilon^2 \Delta \rho &= 0 &\text{ on } (0, 1) \times (0,1) \\ \rho(0,x)&=\rho_0(x),&\;\forall x\in [0,1]\\ \partial_x\rho(0,t)&=\partial_x\rho(1,t)=0,&\;\forall t \in [0,1] \end{array}
If $\rho$ is a solution, then it's mass is preserved, meaning $\int \rho(t,x)dx=\int \rho_0(x)dx = \text{const}$.
To solve it I implemented the Crank-Nicolson method for the temporal derivative and central finite differences for the spatial derivative \begin{equation} \frac{\rho^{n+1}_i - \rho^n_i}{\Delta t} = \frac{\varepsilon^2}{2 (\Delta x)^2} \left[ \left( \rho^{n+1}_{i-1} -2\rho^{n+1}_i +\rho^{n+1}_{i+1} \right) + \left( \rho^{n}_{i-1} -2\rho^{n}_i +\rho^{n}_{i+1} \right) \right] \end{equation}
Here $\rho(t_n,x_i)=\rho_i^n$. Which in matrix form is:
\begin{equation} \left( \begin{array}{c} 1 & -1& 0&\cdots&\cdots&\cdots&0\\ -\alpha & 1+2\alpha & -\alpha & \cdots&\cdots&\cdots &0\\ 0 &-\alpha & 1+2\alpha & -\alpha & \cdots &\cdots&0\\ \ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\ 0 &\cdots&\cdots&\cdots&-\alpha & 1+2\alpha & -\alpha \\ 0 &\cdots&\cdots&\cdots&\cdots& -1& 1\\ \end{array} \right) \left( \begin{array}{c} \rho_{0}^{n+1}\\ \vdots\\ \rho_{i}^{n+1}\\ \rho_{i+1}^{n+1}\\ \vdots\\ \rho_{M}^{n+1}\\ \end{array} \right) = \left( \begin{array}{c} 0 & 0& 0&\cdots&\cdots&\cdots&0\\ \alpha & 1-2\alpha & \alpha & \cdots&\cdots&\cdots &0\\ 0 &\alpha & 1-2\alpha & \alpha & \cdots &\cdots&0\\ \ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\ 0 &\cdots&\cdots&\cdots&\alpha & 1-2\alpha & \alpha \\ 0 &\cdots&\cdots&\cdots&\cdots& 0& 0\\ \end{array} \right) \left( \begin{array}{c} \rho_{0}^{n}\\ \vdots\\ \rho_{i}^{n}\\ \rho_{i+1}^{n}\\ \vdots\\ \rho_{M}^{n}\\ \end{array} \right) \end{equation}
When I use this method, I get a solution that doesn't preserve mass. However, this numerical scheme should preserve it. That is, it should we true that $ m^n=\sum_i\rho_i^n $ doesn't depend on n.
Anyone know what am I doing wrong? Maybe I didn't do the modeling properly?
Edit: For mass to be preserved, we need that when we add theirs rows we get the vector $(1,\ldots,1)$, this is because $$ m(Ax)=m((<R_1,x>,\ldots,<R_n,x>))=\\ <\sum_i R_i,x> $$ Where $R_i$ are the rows of $A$.