Fast Sine-Transform Preconditioning for Global-in-Time Fractional Diffusion

Time-fractional diffusion equations describe subdiffusive transport in heterogeneous media, but their numerical treatment is complicated by the nonlocal Caputo derivative and by the weak singularity that the solution develops at the initial time. We study a global-in-time discretization that combines spectral collocation in time—on the fractional power basis {tℓα}ℓ=0N, evaluated at Chebyshev–Gauss–Lobatto nodes, which reproduces the leading terms of the singular expansion of the solution—with a second-order conservative finite-difference stencil in space that uses harmonic averaging of the diffusivity at the cell faces and therefore remains accurate across discontinuous media. The resulting fully discrete problem is a large, nonsymmetric, dense-in-time linear system whose two-norm condition number grows like the inverse square of the spatial mesh size, so that Krylov subspace iteration without preconditioning stalls under refinement. Exploiting the Kronecker sum structure of the discrete operator, we build a preconditioner by fast diagonalization of the spatial factor through the discrete sine transform. For constant diffusivity the preconditioner reproduces the operator exactly and yields a direct solver; for variable diffusivity it is spectrally equivalent to the operator, and we prove that the eigenvalues of the preconditioned system cluster in a disk centered at one whose radius depends only on the coefficient contrast, and not on the mesh, the number of temporal degrees of freedom, or the fractional order. Numerical experiments in one and two space dimensions confirm second-order spatial accuracy and a preconditioned iteration count that stays flat—twelve iterations from M=32 up to M=1024 in one dimension and eleven up to M=256 per direction in two—while the unpreconditioned count grows by more than two orders of magnitude. In time, the accuracy is spectral until round-off in the ill-conditioned Vandermonde matrix of the power basis takes over: the barrier is reached at N=9,10,13 for α=0.3,0.5,0.7, where the attainable error is about 10−6. A benchmark against the L1 scheme on uniform and graded meshes, the Alikhanov L2-1σ scheme and Grünwald–Letnikov convolution quadrature quantifies when the global approach pays: on forced problems and on modes with κλTα≲2 it reaches a prescribed accuracy one to two orders of magnitude faster and with several times less memory, while for strongly damped modes the fractional power basis converges only algebraically and graded time marching is preferable below a relative error of 10−2.