Mass conservation for heat equation with Neumann conditions

Question

I'm trying to implement Neumann boundary conditions to solve the heat equation with an explicit scheme. I checked mass conservation, but it doesn't seem to hold. Could someone check if my boundary conditions are correctly implemented?

The equation to be computed on the interval $[0,1]$ is:

$$m_t=\nu\Delta m$$ with $m(\cdot,0)=m(\cdot,1)=0$ and the initial condition is Gaussian.

In discretised form:

$$m^{i+1}_j=m^i_j+p(m_{j+1}^i-2m_j^i+m_{j-1}^i)$$ where $p=\nu\frac{dt}{h^2}$. For the BC I used $$\frac{m_1^i-m_{-1}^i}{h^2}=0\quad \text{so}\quad m_1^i=m_{-1}^i$$ and hence $$m_0^{i+1}=m_0^i+2p(m_1^i-m_0^i)$$ So I want to solve: $$m^{i+1}=Am^i$$ where A is the matrix $$ \begin{pmatrix} 1-2p & 2p & 0 & & \cdots& 0 \\ p & 1-2p& p && \\ 0 & p & 1-2p &p &\cdots&0 \\ \vdots & \vdots & \ddots& \ddots & \ddots &\vdots \\ & && p & 1-2p&p\\ 0 & 0 & \cdots& 0& 2p&1-2p \end{pmatrix}$$

Here's the Matlab code

N=15; M=5; 
nu=0.07;

dt=1/(N+1); h=1/(M+1);
p=nu*dt/h^2;

A=eye(M+1) + p * ( diag(ones(M,1),-1) - 2*diag(ones(M+1,1),0) + diag(ones(M,1),1) );

A(1,2)=2*p; A(end,end-1)=2*p; 
mk=zeros(N+1,M+1);
%initial condition
sigma=0.07;
x2=linspace(0,1,M+1);
mk(1,:)=1/sqrt(2*pi*sigma^2)*exp(-(x2-0.5).^2/(2*sigma^2));
for i=1:N
    mk(i+1,:)=A*mk(i,:)';
end

sum(mk,2)

The output is

    4.1097
    4.1099
    4.2142
    4.3373
    4.4505
    4.5464
    4.6257
    4.6895
    4.7428
    4.7836
    4.8155
    4.8423
    4.8624
    4.8793
    4.8884
    4.8993

so mass doesn't seem to be conserved.

I would appreciate if someone could help.

Answer 1

Your scheme will preserve the mass when it is computed by the trapezoidal rule: $$ \text{mass}(i) = \tfrac12(m^i_0 + m^i_M) + \sum_{j=1}^{M-1} m^i_j .$$ One way to see this is that if $w = [0.5,1,\dots,1,0.5]$, then $w A = w$. Hence $w\cdot m^{i+1} = w A m^{i} = w \cdot m^{i}$.