Numeric Solution of Advection-Diffusion Equations by a Discrete Time Random Walk Scheme

We consider a discrete time random walk on a one dimensional lattice. A particle begins at an initial position and at each time step the particle will jump from its current lattice site to a neighbouring lattice site. The probability that a particle on the $i^{\mathrm{th}}$ lattice site will jump to the right on the $n^{\mathrm{th}}$ time step is given by an arbitrary time and space dependent probability, $P_{r}(i,n)$, and the probability of jumping to the left is given by the complement, $P_{l}(i,n)=1-P_{r}(i,n)$. We can recursively define the probability of the particle being at lattice site $i$ at time step $n$, denoted $U(i,n)$, given some initial condition, $U(i,0)$,

$$U(i,n)=P_{r}(i-1,n-1)U(i-1,n-1)+P_{l}(i+1,n-1)U(i+1,n-1).$$

(1)

This is the master equation governing the time evolution of the probability mass function for the process. Here we have considered an infinite lattice, later we will discuss a finite lattice, in which boundary conditions $U(L,n)$ are defined for a boundary set at $x=L$ at all times $n$.

We wish to consider the continuum space and time limit of the discrete time random walk. The appropriate limit in this case is a diffusion limit in which the space and time scales are together taken to zero with the requirement that the following limit exists,

$$D=\lim_{\Delta x,\Delta t\to 0}\frac{\Delta x^{2}}{2\Delta t},$$

(2)

where $\Delta x$ is the spatial grid spacing and $\Delta t$ is the temporal grid spacing [22]. In the diffusion limit the master equation, Eq. (1), will become an advection-diffusion PDE. As such, the master equation itself may be considered as an approximation to the advection-diffusion equation.

First we begin by associating each of the discrete functions, $U(i,n)$ with a continuous function,$u_{\Delta}(x,t)$, such that,

$$u_{\Delta}(i\Delta x,n\Delta t)=U(i,n).$$

(3)

This association is dependent on the lattice spacings, $\Delta x$ and $\Delta t$. The jump probabilities may also be associated with continuum functions, sampled at discrete points,

$$p_{r,\Delta}(i\Delta x,n\Delta t)=P_{r}(i,n),$$

(4)

and similarly for $P_{l}(i,n)$. At each point the continuum jump probabilities sum to one,

$$p_{r,\Delta}(i\Delta x,n\Delta t)+p_{l,\Delta}(i\Delta x,n\Delta t)=1.$$

(5)

At the point $x=i\Delta x$ and $t=n\Delta t$ ,the discrete master equation, Eq. (1) then becomes,

$$u_{\Delta}(x,t)=p_{r,\Delta}(x-\Delta x,t-\Delta t)u_{\Delta}(x-\Delta x,t-\Delta t)+p_{l,\Delta}(x+\Delta x,t-\Delta t)u_{\Delta}(x+\Delta x,t-\Delta t).$$

(6)

Defining a force function,

	$\displaystyle f_{\Delta}(x,t)$	$\displaystyle=p_{r,\Delta}(x,t)-p_{l,\Delta}(x,t),$		(7)
		$\displaystyle=2p_{r,\Delta}(x,t)-1,$		(8)

the master equation can be expressed as,

$$\begin{split}u_{\Delta}(x,t)=&\frac{1}{2}\left(u_{\Delta}(x-\Delta x,t-\Delta t)+u_{\Delta}(x+\Delta x,t-\Delta t)\right)\\ &+\frac{f_{\Delta}(x-\Delta x,t-\Delta t)}{2}\left(u_{\Delta}(x-\Delta x,t-\Delta t)\right)-\frac{f_{\Delta}(x+\Delta x,t-\Delta t)}{2}\left(u_{\Delta}(x+\Delta x,t-\Delta t)\right)\end{split}$$

(9)

Further manipulations give,

$$\begin{split}&\frac{u_{\Delta}(x,t)-u_{\Delta}(x,t-\Delta t)}{\Delta t}=\frac{\Delta x^{2}}{2\Delta t}\left(\frac{u_{\Delta}(x-\Delta x,t-\Delta t)-2u_{\Delta}(x,t-\Delta t)+u_{\Delta}(x+\Delta x,t-\Delta t)}{\Delta x^{2}}\right)\\ &+\frac{\Delta x^{2}}{2\Delta t}\left(\frac{f_{\Delta}(x-\Delta x,t-\Delta t)\left(u_{\Delta}(x-\Delta x,t-\Delta t)\right)-f_{\Delta}(x+\Delta x,t-\Delta t)\left(u_{\Delta}(x+\Delta x,t-\Delta t)\right)}{\Delta x^{2}}\right)\end{split}$$

(10)

We are now set up to take the diffusion limit, that is the limit $\Delta x\to 0$ and $\Delta t\to 0$ such that the limit, Eq. (2) exists. This limit will take the random walk stochastic process to a limit process provided that the limit process exists. In this case the limit process will exist as the original stochastic process is a simple random walk [23]. The particle distribution from the limit process will be given by the limit of the continuum embedding of the discrete process distribution, i.e.,

$$u(x,t)=\lim_{\Delta x,\Delta t\to 0}u_{\Delta}(x,t).$$

(11)

Taking the limit of Eq. (LABEL:eq_nonlim) gives,

$$\begin{split}\frac{\partial u(x,t)}{\partial t}=&D\frac{\partial^{2}u(x,t)}{\partial x^{2}}-2\beta D\frac{\partial}{\partial x}\left(F(x,t)u(x,t)\right),\end{split}$$

(12)

where,

$$\beta F(x,t)=\lim_{\Delta x\to 0}\frac{f_{\Delta}(x,t)}{\Delta x}=\lim_{\Delta x\to 0}\frac{p_{r,\Delta}(x,t)-p_{l,\Delta}(x,t)}{\Delta x}.$$

(13)

This shows that the diffusion limit of the random walk master equation is an advection-diffusion PDE. It should be noted that there is no condition on the form of the force function and this PDE may be a non-linear advection-diffusion equation by considering $F(x,t)$ as an explicit function of $u(x,t)$.

The scheme has been derived from the evolution of a probability density function. In general the integral of the initial condition over the domain will not equal one, and hence cannot be a probability density. As such we note that it is always possible to rescale any initial condition and we can consider the evolution equation to evolve a general non-normalised distribution. Another consequence of the distributional nature of the derivation is that we have a guarantee that the solutions will conserve mass, provided that the boundary conditions preserve mass, and be bound non-negative. Whilst this is often advantageous, a further complication arises if the initial condition is not strictly non-negative. This can not be interpreted as a distribution and hence may not be evolved forward by the master equation. However it is possible to write any generalised function as the difference between two distributions. We simply define two strictly non-negative distributions $u^{+}(x,t)$ and $u^{-}(x,t)$ such that,

$$u(x,t)=u^{+}(x,t)-u^{-}(x,t).$$

(32)

If the strictly non-negative distributions evolve according to the coupled non-linear advection diffusion equations,

	$\displaystyle\frac{\partial u^{+}(x,t)}{\partial t}$	$\displaystyle=\frac{\partial^{2}u^{+}(x,t)}{\partial x^{2}}-\nu\frac{\partial F(x,t,u^{+}-u^{-})u^{+}(x,t)}{\partial x},$		(33)
	$\displaystyle\frac{\partial u^{-}(x,t)}{\partial t}$	$\displaystyle=\frac{\partial^{2}u^{-}(x,t)}{\partial x^{2}}-\nu\frac{\partial F(x,t,u^{+}-u^{-})u^{-}(x,t)}{\partial x},$		(34)

then $u(x,t)$ will evolve according to the non-linear advection diffusion equation,

$$\frac{\partial u(x,t)}{\partial t}=\frac{\partial^{2}u(x,t)}{\partial x^{2}}-\nu\frac{\partial F(x,t,u)u(x,t)}{\partial x}.$$

(35)

Note that the coupling occurs through the non-linear force term. The evolution of the distributions $u^{+}$ and $u^{-}$ may be approximated as coupled DTRWs. The resulting numerical scheme is,

	$\displaystyle U^{+}(i,n)=$	$\displaystyle\frac{U^{+}(i-1,n-1)}{1+\exp(-\frac{\beta\Delta x}{2}\left(F(i-2,n-1)+2F(i-1,n-1)+F(i,n-1)\right))}$		(36)
		$\displaystyle+\frac{U^{+}(i+1,n-1)}{1+\exp(\frac{\beta\Delta x}{2}\left(F(i+2,n-1)+2F(i+1,n-1)+F(i,n-1)\right))}.$
	$\displaystyle U^{-}(i,n)=$	$\displaystyle\frac{U^{-}(i-1,n-1)}{1+\exp(-\frac{\beta\Delta x}{2}\left(F(i-2,n-1)+2F(i-1,n-1)+F(i,n-1)\right))}$		(37)
		$\displaystyle+\frac{U^{-}(i+1,n-1)}{1+\exp(\frac{\beta\Delta x}{2}\left(F(i+2,n-1)+2F(i+1,n-1)+F(i,n-1)\right))}.$

where $U(i,n)=U^{+}(i,n)-U^{-}(i,n)$. Notice that if we take the difference of Eq. (36) and Eq. (37) then we recover the original numerical scheme, Eq. (31) with no appearance of $U^{+}$ or $U^{-}$. Thus we can use the original numerical scheme for the case of initial conditions with mixed sign, provided that the boundary conditions also do not contain $U^{+}$ or $U^{-}$.

Care must be taken with the implementation of boundary conditions so that the underlying random walk remains valid. The boundary conditions thus need to be consistent with both the underlying stochastic process and the given PDE boundary conditions. In general, this excludes some conditions, but does still allow for the implementation of a wide variety. The scheme can accommodate non-conservative boundary conditions, but care needs to be taken that the boundaries do not force a positive distribution to become negative.

If the distance between the approximation from the numerical scheme and the solution to the PDE is bounded for all time steps then the numerical scheme is said to be stable. To show that the scheme is stable it is then sufficient to show that, using the $l_{1}$-norm,

$$\sum_{i}|U(i,n)-u(i\Delta x,n\Delta t)|\leq C,$$

(50)

for all $n$ where $C\in\mathbb{R}^{+}$. It is trivial to see that,

$$\sum_{i}|U(i,n)-u(i\Delta x,n\Delta t)|\leq\sum_{i}|U(i,n)|+\sum_{i}|u(i\Delta x,n\Delta t)|.$$

(51)

From this, provided that the exact solution, $u(x,t)$, is bounded, the stability of the numerical method is assured if the approximation $U(i,n)$ remains bounded.

As $U(i,n)$ is derived from an underlying stochastic process and as such is a density we know that $U(i,n)\geq 0\;\forall i$ and hence,

$$\sum_{i}|U(i,n)|=\sum_{i}U(i,n)$$

(52)

In general stability needs to be shown individually for the problem under consideration. However there are two general classes of problems where we can show unconditional stability. The first class we can consider is a periodic domain with a commensurate periodic potential. In such a case the boundary points are set to the interior point at the opposite boundary. Considering the case of the domain $x\in[0,l]$ and taking a lattice of $M+1$ points with boundary points at $i=0$ and $i=M$, then

$$\begin{split}\sum_{i=1}^{M}U(i,n)=&\sum_{i=1}^{M}U(i,n-1)+P_{r}(0,n-1)U(0,n-1)-P_{l}(1,n-1)U(1,n-1)\\ &+P_{l}(M,n-1)U(M,n-1)-P_{r}(M-1,n-1)U(M-1,n-1).\end{split}$$

(53)

The periodic boundary conditions then imply that $U(0,n)=U(M-1,n)$, and $U(M,n)=U(1,n)$. The periodic potential ensures that the jump probabilities follow the same relation, $P_{r}(0,n)=P_{r}(M-1,n)$, and $P_{l}(M,n)=P_{l}(1,n)$. Thus the total mass is conserved as $\sum_{i=1}^{M}U(i,n)=\sum_{i=1}^{M}U(i,n-1)$ from Eq. (53). As the initial mass is finite the conservation of mass guarantees the stability of the scheme.

The next class of unconditionally stable problems is the case of zero-flux boundary conditions. Here we have the condition that no mass is entering the system at the boundary. For this to be true the continuum equations must satisfy the following conditions,

$$\frac{\partial u(x,t)}{\partial x}\Big{|}_{x=0}=2\beta F(0,t)u(0,t),$$

(54)

and

$$\frac{\partial u(x,t)}{\partial x}\Big{|}_{x=l}=2\beta F(l,t)u(l,t).$$

(55)

Considering the case of the domain $x\in[0,l]$ and taking lattice of $M$ points the boundary conditions are enforced by setting ghost points at $i=0$ and $i=M+1$, then

$$\begin{split}\sum_{i=1}^{M}U(i,n)=&\sum_{i=1}^{M}U(i,n-1)+P_{r}(0,n-1)U(0,n-1)-P_{l}(1,n-1)U(1,n-1)\\ &+P_{l}(M+1,n-1)U(M+1,n-1)-P_{r}(M,n-1)U(M,n-1).\end{split}$$

(56)

Zero-flux boundary conditions are implemented by setting the ghost points such that the total mass will be conserved,

	$\displaystyle U(0,n)$	$\displaystyle=\frac{P_{l}(1,n)}{P_{r}(0,n)}U(1,n),$		(57)
	$\displaystyle U(M+1,n)$	$\displaystyle=\frac{P_{r}(M,n)}{P_{l}(M+1,n)}U(M,n).$		(58)

Using the two point quadrature Boltzmann weights, Eq. (30) and taking the limit $\Delta x\to 0$, Eq. (57) recovers Eq. (54) as required. Again as the total mass is conserved the numerical scheme will be stable.

In order to estimate the accuracy of the DTRW numerical scheme we consider the scheme as if it were derived from a Taylor series. Considering a non-linear advection diffusion equation of the form,

$$\frac{\partial u(x,t)}{\partial t}=D\frac{\partial^{2}u(x,t)}{\partial x^{2}}-2\beta D\frac{\partial F(x,t,u)u(x,t)}{\partial x}.$$

(59)

The numerical scheme that approximates this equation is then,

$$u(x,t+\Delta t)=p_{l}(x+\Delta x,t)u(x+\Delta x,t)+p_{r}(x-\Delta x,t)u(x-\Delta x,t),$$

(60)

where, as before, the jump probabilities are given by Boltzmann weights,

$$p_{r}(x,t)=\frac{\exp\left({\beta\frac{\Delta x}{2}(F(x+\Delta x,t)+F(x,t))}\right)}{\exp{\left(\beta\frac{\Delta x}{2}(F(x+\Delta x,t)+F(x,t))\right)}+\exp{\left(-\beta\frac{\Delta x}{2}(F(x-\Delta x,t)+F(x,t))\right)}},$$

(61)

and

$$p_{l}(x,t)=\frac{\exp\left({-\beta\frac{\Delta x}{2}(F(x-\Delta x,t)+F(x,t))}\right)}{\exp{\left(\beta\frac{\Delta x}{2}(F(x+\Delta x,t)+F(x,t))\right)}+\exp{\left(-\beta\frac{\Delta x}{2}(F(x-\Delta x,t)+F(x,t))\right)}}.$$

(62)

Taking the Taylor expansion with respect to $\Delta x$, of the terms on the right hand side of Eq. (60), around $\Delta x=0$, to order $\mathcal{O}(\Delta x^{3})$ we have

$$\begin{split}p_{r}(x-\Delta x,t)u(x-\Delta x,t)=&\frac{1}{2}u(x,t)+\frac{1}{2}\left(\beta F(x,t)u(x,t)-\frac{\partial u(x,t)}{\partial x}\right)\Delta x\\ &-\frac{1}{4}\left(2\beta u(x,t)\frac{\partial F(x,t)}{\partial x}+2\beta F(x,t)\frac{\partial u(x,t)}{\partial x}-\frac{\partial^{2}u(x,t)}{\partial x^{2}}\right)\Delta x^{2}\\ &-\frac{1}{24}\left(4\beta^{3}F(x,t)^{3}u(x,t)-12\beta\frac{\partial F(x,t)}{\partial x}\frac{\partial u(x,t)}{\partial x}\right.\\ &\left.-9\beta u(x,t)\frac{\partial^{2}F(x,t)}{\partial x^{2}}-6\beta F(x,t)\frac{\partial^{2}u(x,t)}{\partial x^{2}}+2\frac{\partial^{3}u(x,t)}{\partial x^{3}}\right)\Delta x^{3}+\mathcal{O}(\Delta x^{4}),\end{split}$$

(63)

and

$$\begin{split}p_{l}(x+\Delta x,t)u(x+\Delta x,t)=&\frac{1}{2}u(x,t)-\frac{1}{2}\left(\beta F(x,t)u(x,t)-\frac{\partial u(x,t)}{\partial x}\right)\Delta x\\ &-\frac{1}{4}\left(2\beta u(x,t)\frac{\partial F(x,t)}{\partial x}+2\beta F(x,t)\frac{\partial u(x,t)}{\partial x}-\frac{\partial^{2}u(x,t)}{\partial x^{2}}\right)\Delta x^{2}\\ &+\frac{1}{24}\left(4\beta^{3}F(x,t)^{3}u(x,t)-12\beta\frac{\partial F(x,t)}{\partial x}\frac{\partial u(x,t)}{\partial x}\right.\\ &\left.-9\beta u(x,t)\frac{\partial^{2}F(x,t)}{\partial x^{2}}-6\beta F(x,t)\frac{\partial^{2}u(x,t)}{\partial x^{2}}+2\frac{\partial^{3}u(x,t)}{\partial x^{3}}\right)\Delta x^{3}+\mathcal{O}(\Delta x^{4}).\end{split}$$

(64)

Substituting equations (63) and (64) into equation (60), dividing through by $\Delta t$ and collecting term appropriately we obtain,

$$\frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}=\frac{\Delta x^{2}}{2\Delta t}\frac{\partial^{2}u(x,t)}{\partial x^{2}}-\beta\frac{\Delta x^{2}}{\Delta t}\frac{\partial F(x,t)u(x,t)}{\partial x}+\mathcal{O}(\Delta x^{4}).$$

(65)

Taking the Taylor expansion of the left hand side with respect to $\Delta t$ we have,

$$\frac{\partial u(x,t)}{\partial t}+\mathcal{O}(\Delta t)=\frac{\Delta x^{2}}{2\Delta t}\frac{\partial^{2}u(x,t)}{\partial x^{2}}-\beta\frac{\Delta x^{2}}{\Delta t}\frac{\partial F(x,t)u(x,t)}{\partial x}+\mathcal{O}(\Delta x^{4}),$$

(66)

with $D=\Delta x^{2}/2\Delta t$ this reduces to,

$$\frac{\partial u(x,t)}{\partial t}=D\frac{\partial^{2}u(x,t)}{\partial x^{2}}-2\beta D\frac{\partial F(x,t)u(x,t)}{\partial x}+\mathcal{O}(\Delta t)+\mathcal{O}(\Delta x^{4}).$$

(67)

Finally noting that $\mathcal{O}(\Delta t)=\mathcal{O}(\Delta x^{2})$ we see that the scheme will be convergent with order $\Delta x^{2}$.