Numerical Solution of the Diffusion Equation with Constant Concentration Boundary Conditions - Setup

Let's consider the diffusion equation with boundary conditions $c(0,t) = c(1,t) = 0$, that is, the concentration at the boundaries is held at zero. Physically, this could correspond to our system being in contact at its boundaries with a very large reservoir containing a very small concentration of the chemical. Indeed, imagine that the reservoir is so large that, even when some of our chemical escapes into it, the concentration within the reservoir stays approximately equal to zero. We enforce these boundary conditions in our finite difference scheme by explicitly requiring that $c_{1,j} = c_{N_x,j} = 0$ for all $j$.

Before giving the Matlab code to numerically solve the diffusion equation, let's rewrite this equation (2) as

$$\frac{\partial c}{\partial t} = - \frac{\partial}{\partial x} \left( -\frac{\partial c}{\partial x} \right).$$

We identify the quantity inside the parentheses as the flux $J$:

$$J(x,t) \equiv -\frac{\partial c}{\partial x}.$$ (7)

Since this is a one-dimensional problem, the flux $J(x,t)$ is the number of particles per unit time crossing $x$ at time $t$. A positive flux corresponds to a net flow in the $+x$-direction, while a negative flux corresponds to a net flow in the $-x$-direction. Because it's instructive, let's keep track of the total number of particles which have crossed the boundaries at $x=0,1$. At $x=1$, this is

$$s_1(t) \equiv \int_0^t J(1,t') dt'.$$ (8)

At $x=0$, we need to take into account the fact that a negative flux means that a positive number of particles are leaving the domain. Thus, the total number of particles which have crossed the boundary at $x=0$ is

$$s_0(t) \equiv -\int_0^t J(0,t') dt'.$$ (9)

The total number of particles within the box is given by

$$s(t) \equiv \int_0^1 c(x,t) dx.$$ (10)

By conservation of mass, it is expected that $s(t) + s_0(t) + s_1(t)$ should be a constant.

For definiteness, we'll consider the initial condition

$$c(x,t) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(x-\mu)^2/(2 \sigma^2)},$$ (11)

i.e., a Gaussian with mean $\mu$ and variance $\sigma^2$. This has been normalized so that

$$\int_{-\infty}^\infty c(x,t) dx = 1.$$ (12)

We'll choose this Gaussian to be narrow enough that the initial condition satisfies the boundary conditions at $x=0,1$ for all practical purposes, so that

$$\int_0^1 c(x,t) dx \approx 1.$$