Extended block Krylov subspace approaches for solving large-scale linear system of fractional DEs

This study addresses the numerical solution of large-scale linear systems of fractional differential equations (LSFDEs) featuring a low-rank constant term a class of problems not previously investigated in the literature. We propose two novel numerical approaches for solving such systems. The first method exploits the integral representation of the exact solution and employs a Krylov-based approximation to compute the action of the matrix Mittag–Leffler function on a block of vectors. The second approach projects the original high-dimensional fractional system onto an extended block Krylov subspace, reducing it to a significantly smaller fractional differential system. This reduced system is then solved using either a tailored implementation of the Grünwald-Letnikov scheme or a fractional backward differentiation formula, both adapted to the projected setting. The resulting low-rank approximate solution is iteratively refined by expanding the projection subspace until a prescribed tolerance is achieved. We derive explicit expressions for the residual and error norms and establish associated convergence estimates. To validate the computational efficiency and accuracy of the proposed methods, we conduct extensive numerical experiments on several benchmark problems. The results demonstrate that our approaches substantially reduce computational time while maintaining high numerical precision, outperforming existing conventional solvers for large-scale fractional systems.