In this paper, we propose a numerical approach for the fractional reaction sub-diffusion equation with a Caputo-Hadamard derivative of fractional order \(\alpha \in (0,1)\) , whose solutions typically exhibit singular behavior at the initial time. The approach involves the use of an L2-type discrete operator of order \(3-\alpha \) to approximate the fractional derivative on nonuniform temporal meshes, while a standard second-order difference method is employed for uniform spatial meshes. Detailed discussions are presented on the truncation errors and coefficients of the proposed discrete fractional derivative operator. The stability as well as the accuracy of the resulting numerical scheme is rigorously analyzed on a unique nonuniform mesh. Several numerical examples are provided to validate the theoretical results and to demonstrate the efficiency of the proposed method.
