Random walk models of anisotropic diffusion on rectangular and hexagonal lattices

TL;DR

This work addresses simulating anisotropic diffusion in two dimensions by constructing equivalent discrete-space discrete-time random walk models for a diffusion equation with a spatially invariant tensor $\mathbf{D}$ under no-flux boundaries. A vertex-centered finite-volume discretisation on rectangular and hexagonal lattices, combined with forward Euler time-stepping, yields a simple Markov-chain representation with transition matrix $\mathbf{P}=\mathbf{I}+\tau\mathbf{C}$ and explicit transition probabilities. Rectangular lattices impose tensor-dependent constraints (e.g., $\det(\mathbf{D})\ge 3D_{xy}^2$) for valid probabilities, whereas flat-top and pointy-top hexagonal lattices accommodate any SPD diffusion tensor, ensuring a stochastic mapping. Across multiple test tensors and initial conditions, the random-walk models show good agreement with the deterministic PDE solutions, quantified by low mean-squared error, and the authors provide MATLAB code for replication. This framework offers a practical, analytically tractable approach to simulate anisotropic diffusion with intrinsic randomness on common lattice geometries.

Abstract

The diffusive transport of particles in anisotropic media is a fundamental phenomenon in computational, medical and biological disciplines. While deterministic models (partial differential equations) of such processes are well established, their inability to capture inherent randomness, and the assumption of a large number of particles, hinders their applicability. To address these issues, we present several equivalent (discrete-space discrete-time) random walk models of diffusion described by a spatially-invariant tensor on a two-dimensional domain with no-flux boundary conditions. Our approach involves discretising the deterministic model in space and time to give a homogeneous Markov chain governing particle movement between (spatial) lattice sites over time. The spatial discretisation is carried out using a vertex-centred element-based finite volume method on rectangular and hexagonal lattices, and a forward Euler discretisation in time yields a nearest-neighbour random walk model with simple analytical expressions for the transition probabilities. For each lattice configuration, analysis of these expressions yields constraints on the time step duration, spatial steps and diffusion tensor to ensure the probabilities are between zero and one. We find that model implementation on a rectangular lattice can be achieved with a constraint on the diffusion tensor, whereas a hexagonal lattice overcomes this limitation (no restrictions on the diffusion tensor). Overall, the results demonstrate good visual and quantitative (mean-squared error) agreement between the deterministic model and random walk simulations for several test cases. All results are obtained using MATLAB code available on GitHub (https://github.com/lukefilippini/Filippini2025).

Figures (7)

• Figure 1: Representations of diffusion in anisotropic media. (a) Particle diffusion in two-dimensional anisotropic media. The movement of particles (orange circles) occurs primarily along aligned structures (solid grey boxes) of the surrounding environment (indicated by dashed arrows). (b) Representation of the tensor describing diffusion in (a) using an ellipse. The eigenvalues ($\lambda_{1}$ and $\lambda_{2}$) and eigenvectors ($\mathbf{v}_{1}$ and $\mathbf{v}_{2}$) of the eigendecomposition $\mathbf{D}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\mathsf{T}}$ of the diffusion tensor $\mathbf{D}$ correspond to the length and magnitude of the ellipse axes (large arrows), respectively.
• Figure 2: Rectangular and hexagonal lattice configurations. (a) Rectangular and (b)--(c) hexagonal lattices for the spatial discretisation of the deterministic model (\ref{['eq:diffusion_eq']})--(\ref{['eq:D_tensor']}). The problem domain $\Omega = [0,L_{x}] \times [0,L_{y}]$ is partitioned exactly by the (a) rectangular lattice (defined by structured rectangular elements), whereas the (b) flat-top and (c) pointy-top hexagonal lattices (defined by structured triangular elements) approximate $\Omega$. For each lattice configuration, $\delta_{x}$ and $\delta_{y}$ represent the spatial steps between lattice sites in the $x$ and $y$ directions, respectively. Around each lattice site $\mathbf{x}_{i,j}=(x_{i},y_{j})$ (black circles), we construct a control volume $\Omega_{i,j}$ (shaded orange) using a vertex-centred approach, implying the centroids of the elements $E_{m}$ (shaded grey), which have the lattice site $\mathbf{x}_{i,j}$ as a vertex, correspond to the vertices of the control volume boundary $\partial\Omega_{i,j}$ (shaded purple). An approximation for the flux across the control volume boundary is obtained by separating $\partial\Omega_{i,j}$ into line segments $\sigma \in \mathcal{E}_{i,j}$, where $\mathcal{E}_{i,j}$ denotes the set of line segments which constitute $\partial\Omega_{i,j}$, computing the outward (relative to $\sigma$) flux at each segment midpoint $\bar{\mathbf{x}}_{\sigma}$ (black crosses) in the direction of the unit normals $\hat{\mathbf{n}}_{\sigma}$ and multiplying by the segment lengths $\ell_{\sigma}$.
• Figure 3: Prototype tensor constraints for particle diffusion on a rectangular lattice. (a) Proportion (between zero and one) of the unit semicircle within which the angle of rotation $0 < \theta < \pi$ satisfies (\ref{['eq:rec_theta']}) for implementation of the deterministic and random walk models (with the prototype diffusion tensor (\ref{['eq:prototype_tensor']})) on a rectangular lattice. If $\lambda_{y}/3 \leq \lambda_{x} \leq 3\lambda_{y}$, where $\lambda_{x}$ and $\lambda_{y}$ are the tensor eigenvalues, there is no restriction on $\theta$ (labelled region with bounds denoted by dashed black lines). Otherwise, the proportion of the unit semicircle (range of appropriate angles) decreases as $\lambda_{x}$ and $\lambda_{y}$ diverge away from each other. (b) Angular regions (shaded) of the unit semicircle within which $\theta$ satisfies (\ref{['eq:rec_theta']}) for eigenvalues $\lambda_{x} = 0.08$ and $\lambda_{y} = 0.02$ (labelled on (a)). The three prototype diffusion tensors used in this work are generated using these eigenvalues and the respective angles $\theta = \pi/24$, $\theta = 5\pi/12$ and $\theta=3\pi/4$ (dashed red lines). The latter angle does not satisfy the constraint (\ref{['eq:rec_theta']}).
• Figure 4: Simple and complex initial conditions. Two-dimensional (a) simple and (b) complex initial particle densities (at time $t = 0$) for the deterministic and random walk models implemented on a flat-top hexagonal lattice. The simple and complex initial conditions are defined in (\ref{['eq:f_simple']}) and (\ref{['eq:f_complex']}), respectively, where the latter utilises the density aggregation algorithm (\ref{['eq:flat-top_complex']}). Simulation results (two-dimensional model comparisons) using these initial conditions are presented in Figures \ref{['fig:2D_flat-top_simple']} and \ref{['fig:2D_flat-top_complex']}, respectively. Parameters: $L_{x} = 1$, $L_{y} = 0.8$, $N_{x} = 202$, $N_{y} = 89$, $\delta_{x} = 0.01$, $\delta_{y} \approx 0.0182$, $\kappa = 25$, $\gamma = 0.5$.
• Figure 5: Two-dimensional comparisons of anisotropic diffusion models (simple initial condition). Two-dimensional comparisons of particle density obtained from the deterministic (left column of panels) and random walk (right column of panels) models implemented on a flat-top hexagonal lattice. The solutions are plotted for times $t = 0.03, 0.13, 0.25$ (corresponding to each row of panels in descending order) where the deterministic particle densities are obtained from the solution of the spatially-discretised model (\ref{['eq:density_ODEs']}), and the stochastic particle densities are averaged across five simulations of the random walk model (see Algorithm \ref{['alg:random_walk']}). A simple (rectangular) initial condition (\ref{['eq:f_simple']}) is used: $f(\mathbf{x}_{i,j}) = 1$ if $3L_{x}/10 \leq x_{i} \leq 7L_{x}/10$ and $3L_{y}/10 \leq y_{j} \leq 7L_{y}/10$, otherwise $f(\mathbf{x}_{i,j}) = 0$. Parameters: $L_{x} = 1$, $L_{y} = 0.8$, $N_{x}=202$, $N_{y}=89$, $\delta_{x} = 0.01$, $\delta_{y} \approx 0.0182$, $N_{p}^{*}=449189$, $N_{s}=5$, $N_{t}=2332$, $T=1$, $\tau \approx 0.00043$, $\lambda_{x} = 0.08$, $\lambda_{y}=0.02$, $\theta = 5\pi/12$, $D_{xx} \approx 0.024$, $D_{xy}=0.015$ and $D_{yy} \approx 0.076$. Colour map sourced from https://matplotlib.org/cmocean/thyng_et_al_2016.