Consider the equation: $$ \partial_tu=D\partial_{xx}^2u $$ with reflecting boundary condition at $x=0$ and with $u(x,0)=\delta(x)$ as an initial distribution.
First question: How should I understand a $\delta(x)$ in a positive domain with reflecting boundary conditions at $x=0$?
If we use a forward Euler integration scheme: $$ u^{t+1}_x=u^t_x+\Delta t(\Delta u)_x^t $$ where: $$ (\Delta u)_x^t=D\frac{u^t_{x+1}+u^t_{x-1}-2u^t_{x}}{(\Delta x)^2} $$ A very important point (in my opinion) is how we introduce to this scheme the reflecting boundary conditions:
- $u_1=u_{-1}$ (central difference $=0$)
- $u_1=u_0$ (forward difference $=0$)
- $u_0=u_{-1}$ (backward difference $=0$)
For each option I would like to compute the conservation of mass/probability of this scheme: $$ U(t)=\sum_x u^t_x\Delta x $$ From now on, I consider for simplicity $\mathbf{\Delta x = \Delta t = 1}$ and $\mathbf{D=1}$. Thus, we have: $$ u^0_0 = 1;\quad u^0_x=0 \quad x\geq1 \implies U(0)=1 $$ If we expect from reflecting boundary conditions conservation of mass/probability, after one time-step, we would like: $$ U(0)= u^0_0 + u^0_1 = u^1_0 + u^1_1 = u^0_0 + (\Delta u)^0_0+ u^0_1 +(\Delta u)^0_1 = U(1)\implies (\Delta u)^0_0+(\Delta u)^0_1 = 0 $$ Let's check this for the three possible ways to implement the reflecting boundary conditions.
- $u_1=u_{-1}$ $$ (\Delta u)^0_0= u_{1}+u_{-1}-2u_{0} = 0 + 0 -2\cdot 1=-2; \qquad (\Delta u)^1_0=u_{2}+u_{0}-2u_{1}= 0 + 1 -2\cdot 0 = 1 $$ $$ \implies (\Delta u)^0_0+(\Delta u)^0_1 = -1 \neq 0 $$
- $u_1=u_0$ (I assume $u_{-1}=0$) $$ (\Delta u)^0_0= 1 + 0 -2\cdot 1=-1; \qquad (\Delta u)^1_0= 0 + 1 -2\cdot 1 = -1 \implies (\Delta u)^0_0+(\Delta u)^0_1 = -2 \neq 0 $$
- $u_0=u_{-1}$ $$ (\Delta u)^0_0= 0 + 1 -2\cdot 1=-1; \qquad (\Delta u)^1_0= 0 + 1 -2\cdot 0 = 1 \implies (\Delta u)^0_0+(\Delta u)^0_1 = 0 $$
Therefore, the third way is the only one that conserve the total mass/probability after one timestep, but I would have suggested that option 1. is the correct way to do it, as it more close to the mathematical definition of a derivative. Also, I had simulated many times with 1. with no problem (actually, after few time-steps you match the expected analytical solution).
Second question: Could anyone elucidate why this could seems very unintuitive?
Third question: What do you believe is the appropriate way to implement reflecting boundary conditions?
