An Efficient Solver for Fractional Diffusion on Unbounded Combs with Exact Absorbing Boundary Conditions

Despite its importance in modeling subdiffusion in fractal and heterogeneous media, a rigorous and computational scheme for solving the fractional diffusion equation on generalized comb structures over unbounded domains has remained elusive, mainly due to the nonlocal memory effect and slow spatial decay of solutions. To the best of our knowledge, we address this long-standing gap by presenting a fully integrated framework that simultaneously resolves both challenges. We derive the governing equation from constitutive relations and establish exact absorbing boundary conditions (ABCs) for the multi-skeleton comb model, a result absent in prior work. A transparent Dirichlet-to-Neumann (DtN) map, constructed via Laplace analysis, rigorously handles skeletal Dirac delta singularities and eliminates spurious reflections without empirical parameters. Furthermore, we propose a novel structure-preserving finite difference scheme that applies the sum-of-exponentials (SOE) approximation not only to the interior Caputo derivative but also to the convolution kernels arising from the ABCs. This yields a dramatic reduction in computational complexity, from quadratic O(Nt2) to quasi-linear O(NtlogNt), while preserving the physics of anomalous transport. We prove the well-posedness, unconditional stability, and convergence of the method. Numerical results confirm theoretical error estimates and show excellent agreement between simulated particle distributions, mean square displacement profiles, and exact asymptotics, validating both accuracy and robustness. The speedup (CPU time ratio Direct/Fast) is about 1.00×–1.23× for Nt=5000 in our tests. Our approach sets a new benchmark for simulating anomalous dynamics in fractal-inspired media.