I don't understand the different behaviour of the advection-diffusion equation when I apply different boundary conditions. My motivation is the simulation of a real physical quantity (particle density) under diffusion and advection. Particle density should be conserved in the interior unless it flows out from the edges. By this logic, if I enforce Neumann boundary conditions the ends of the system such as $\frac{\partial \phi}{\partial x}=0$ (on the left and the right sides) then the system should be *"closed"* i.e. if the **flux** at the boundary is zero then no particles can escape.
For all the simulations below, I have applied the Crank-Nicolson discretization to the advection-diffusion equation and all simulation have $\frac{\partial \phi}{\partial x}=0$ boundary conditions. However, for the first and last rows of the matrix (the boundary condition rows) I allow $\beta$ to be changed independently of the interior value. This allows the end points to be fully implicit.
Below I discuss 4 different configurations, only one of them is what I expected. At the end I discuss my implementation.
Diffusion only limit
--------------------
Here the advection terms are turned off by setting the velocity to zero.
**Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points**
![Diffusion only (Neumann boundaries with beta=0.5)][1]
The quantity is not conserved as can be seen by the pulse area reducing.
**Diffusion only, with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries**
![Diffusion only (Neumann boundaries with beta=0.5 for interior, beta=1 fully implicit) the boundaries][2]
By using fully implicit equation on the boundaries I achieve what I expect: **no particles escape**. You can see this by the area being conserved as the particle diffuse. **Why** should the choice of $\beta$ at the boundary points influence the physics of the situation? Is this a bug or expected?
Diffusion and advection
-----------------------
When the advection term is included, the value of $\beta$ at the boundaries does not seem to influence the solution. However, for all cases when the boundaries seem to be *"open"* i.e. particles can escape the boundaries. Why is this the case?
**Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at all points**
![Advection-Diffusion (Neumann boundaries with beta=0.5)][4]
**Advection and Diffusion with $\boldsymbol{\beta}$=0.5 (Crank-Niscolson) at interior points, and $\boldsymbol{\beta}$=1 (full implicit) at the boundaries**
![Advection and Diffusion (Neumann boundaries with beta=0.5 for interior, beta=1 fully implicit) the boundaries][5]
Implementation of the advection-diffusion equation
--------------------------------------------------
Starting with the advection-diffusion equation,
$
\frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial x^2} + \boldsymbol{v}\frac{\partial \phi}{\partial x}
$
Writing using Crank-Nicolson gives,
$
\frac{\phi_{j}^{n+1} - \phi_{j}^{n}}{\Delta t} = D \left[ \frac{1 - \beta}{(\Delta x)^2} \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + \frac{\beta}{(\Delta x)^2} \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) \right] + \boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right) \right]
$
Note that $\beta$=0.5 for Crank-Nicolson, $\beta$=1 for fully implicit, and, $\beta$=0 for fully explicit.
To simplify the notation let's make the substitution,
$
s = D\frac{\Delta t}{(\Delta x)^2} \\
r = \boldsymbol{v}\frac{\Delta t}{2 \Delta x}
$
and move the known value $\phi_{j}^{n}$ of the time derivative to the right-hand side,
$
\phi_{j}^{n+1} = \phi_{j}^{n} + s \left( 1-\beta \right) \left( \phi_{j-1}^{n} - 2\phi_{j}^{n} + \phi_{j+1}^{n} \right) + s \beta \left( \phi_{j-1}^{n+1} - 2\phi_{j}^{n+1} + \phi_{j+1}^{n+1} \right) +
r \left( 1 - \beta \right) \left( \phi_{j+1}^{n} - \phi_{j-1}^{n} \right) + r \beta \left( \phi_{j+1}^{n+1} - \phi_{j-1}^{n+1} \right)
$
Factoring the $\phi$ terms gives,
$
\underbrace{\beta(r - s)\phi_{j-1}^{n+1} + (1 + 2s\beta)\phi_{j}^{n+1} -\beta(s + r)\phi_{j+1}^{n+1}}_{\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}}} = \underbrace{ (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}}_{\boldsymbol{M\cdot}\boldsymbol{\phi^n}}
$
which we can write in matrix form as $\boldsymbol{A}\cdot\boldsymbol{\phi^{n+1}} = \boldsymbol{M}\cdot\boldsymbol{\phi^{n}}$ where,
$ \boldsymbol{A} =
\left(
\begin{matrix}
1+2s\beta & -\beta(s + r) & & 0 \\
\beta(r-s) & 1+2s\beta & -\beta (s + r) & \\
& \ddots & \ddots & \ddots \\
& \beta(r-s) & 1+2s\beta & -\beta (s + r) \\
0 & & \beta(r-s) & 1+2s\beta \\
\end{matrix}
\right)
$
$
\boldsymbol{M} =
\left(
\begin{matrix}
1-2s(1-\beta) & (1 - \beta)(s + r) & & 0 \\
(1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\
& \ddots & \ddots & \ddots \\
& (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\
0 & & (1 - \beta)(s - r) & 1-2s(1-\beta) \\
\end{matrix}
\right)
$
Applying Neumann boundary conditions
------------------------------------
**NB is working through the derivation again I think I have spotted the error. I assumed a fully implicit scheme ($\beta$=1) when writing the finite difference of the boundary condition. If you assume a Crank-Niscolson scheme here the complexity become too great and I could not solve the resulting equations to eliminate the nodes which are outside the domain. However, it would appear possible, there are two equation with two unknowns, but I couldn't manage it. This probably explains the difference between the first and second plots above. I think we can conclude that only the plots with $\beta$=0.5 at the boundary points are valid.**
Assuming the flux at the left-hand side is known (assuming a fully implicit form),
$
\frac{\partial\phi_1^{n+1}}{\partial x} = \sigma_L
$
Writing this as a centred-difference gives,
$
\frac{\partial\phi_1^{n+1}}{\partial x} \approx \frac{\phi_2^{n+1} - \phi_0^{n+1}}{2\Delta x} = \sigma_L
$
therefore,
$
\phi_0^{n+1} = \phi_{2}^{n+1} - 2 \Delta x\sigma_L
$
Note that this introduces a node $\phi_0^{n+1}$ which is outside the domain of the problem. This node can be eliminated by using a second equation. We can write the $j=1$ node as,
$
\beta(r - s)\phi_0^{n+1} + (1+2s\beta)\phi_1^{n+1} - \beta(s+r)\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n}
$
Substituting in the value of $\phi_0^{n+1}$ found from the boundary condition gives the following result for the $j$=1 row,
$
(1+2s\beta)\phi_1^{n+1} - 2s\beta\phi_2^{n+1} = (1-\beta)(s - r)\phi_{j-1}^{n} + (1-2s[1-\beta])\phi_{j}^{n} + (1-\beta)(s+r)\phi_{j+1}^{n} + 2\beta(r-s)\Delta x\sigma_L
$
Performing the same procedure for the final row (at $j$=$J$) yields,
$
-2s\beta\phi_{J-1}^{n+1} + (1+2s\beta)\phi_J^{n+1} = (1-\beta)(s - r)\phi_{J-1}^{n} + (1 - 2s(1-\beta))\phi_{J}^{n} + 2\beta(s+r)\Delta x\sigma_R
$
Finally making the boundary rows implicit (setting $\beta$=1) gives,
$
(1+2s)\phi_1^{n+1} - 2s\phi_2^{n+1} = \phi_{j-1}^{n} + 1\phi_{j}^{n} + 2(r-s)\Delta x\sigma_L
$
$
-2s\phi_{J-1}^{n+1} + (1+2s)\phi_J^{n+1} = \phi_{J}^{n} + 2(s+r)\Delta x\sigma_R
$
Therefore with Neumann boundary conditions we can write the matrix equation, $\boldsymbol{A}\cdot\phi^{n+1} = \boldsymbol{M}\cdot\phi^{n} + \boldsymbol{b_N}$,
where,
$
\boldsymbol{A} =
\left(
\begin{matrix}
1+2s & -2s & & 0 \\
\beta(r-s) & 1+2s\beta & -\beta (s + r) & \\
& \ddots & \ddots & \ddots \\
& \beta(r-s) & 1+2s\beta & -\beta (s + r) \\
0 & & -2s & 1+2s \\
\end{matrix}
\right)
$
$
\boldsymbol{M} =
\left(
\begin{matrix}
1 & 0 & & 0 \\
(1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) & \\
& \ddots & \ddots & \ddots \\
& (1 - \beta)(s - r) & 1-2s(1-\beta) & (1 - \beta)(s + r) \\
0 & & 0 & 1 \\
\end{matrix}
\right)
$
$
\boldsymbol{b_N} = \left(
\begin{matrix}
2 (r - s) \Delta x \sigma_L & 0 & \ldots & 0 & 2 (s + r) \Delta x \sigma_R
\end{matrix}
\right)^{T}
$
My current understanding
------------------------
* I think the difference between the first and second plots is explained by noting the error outlined above.
* <s>Regarding the conservation of the physical quantity. I believe the cause is that, [as pointed out here][3], the advection equation in the form I have written it doesn't allow propagation in the reverse direction so the wave just passes through **even with zero-flux boundary conditions**. My initial intuition regarding conservation only applied when advection term is zero (this is solution in plot number 2 where the area is conserved).</s>
* Even with **Neumann zero-flux** boundary conditions $\frac{\partial \phi}{\partial x} = 0$ the mass can still leave the system, this is because the correct boundary conditions in this case are **Robin** boundary conditions in which the *total* flux is specified $j = D\frac{\partial \phi}{\partial x} + \boldsymbol{v}\phi = 0$
Would you agree?
[1]: https://i.sstatic.net/XGWYQ.gif
[2]: https://i.sstatic.net/un9MJ.gif
[3]: http://scicomp.stackexchange.com/questions/5425/strange-oscillation-when-solving-the-advection-equation-by-finite-difference-wit
[4]: https://i.sstatic.net/1Ck7s.gif
[5]: https://i.sstatic.net/e1WhI.gif
[6]: http://scicomp.stackexchange.com/questions/5355/writing-the-poisson-equation-finite-difference-matrix-with-neumann-boundary-cond
[7]: http://scicomp.stackexchange.com/questions/5425/strange-oscillation-when-solving-the-advection-equation-by-finite-difference-wit