Fractional Partial Differential Equations (FPDEs) are essential for modeling complex systems across various scientific and engineering areas. However, efficiently solving FPDEs presents significant computational challenges due to their inherent memory effects, often leading to increased execution times for numerical solutions. This study proposes a highly parallelizable hybrid computational approach that combines the Finite Element Method (FEM) for spatial discretization with Numerical Inversion of the Laplace Transform (NILT) for time-domain solutions, optimized for execution on Graphics Processing Units (GPUs). The NILT method’s high parallelizability,
stemming from the independence of its series terms, combined with the robust spatial discretization provided by FEM, enables the efficient and accurate solution of FPDEs on GPUs, demonstrating substantial performance improvements over traditional CPU-based implementations. We observe a generalized pattern in execution time behavior that accounts for both the number of nodes and
the number of NILT terms. Specifically, execution time scales quadratically with the number of nodes, while showing only a logarithmic marginal increase with the number of NILT terms. These behaviors not only enables consistent performance assessment but also highlights potential areas for algorithm optimization. Validation against exact solutions of fractional diffusion and wave equations, employing Caputo, modified Caputo-Fabrizio, and modified Atangana-Baleanu derivatives, demonstrates the accuracy and convergence of the hybrid FEM-NILT method. Notably, the exact solutions of wave equation are novel in literature. The results highlight the method’s potential for enabling high-precision, large-scale simulations in fractional calculus applications, thereby advancing
computational capabilities and efficiency in the field.