On the numerical instability of the Riemann–Liouville fractional derivative: analysis of errors and convergence behaviour

Abstract

This paper develops a second-order numerical approximation for the Riemann–Liouville fractional derivative using a modified product trapezoidal rule. The method is derived from the fractional integral formulation via piecewise linear interpolation, yielding a discrete convolution depending only on function values. A rigorous error analysis establishes  accuracy under smoothness assumptions. However, numerical experiments on a test problem reveal significant discrepancies between theory and computation. The results exhibit instability, rapid error growth, and negative convergence orders for fractional orders. The divergence intensifies as the mesh is refined, indicating a breakdown of numerical consistency. These findings reveal notable discrepancies between theoretical predictions and computational performance, highlighting inherent limitations of the Riemann–Liouville derivative operator in applications.