Do we need decay-preserving error estimate for solving parabolic equations with initial singularity?

Solutions exhibiting weak initial singularities arise in various equations, including diffusion and subdiffusion equations. When employing the well-known L1 scheme to solve subdiffusion equations with weak singularities, numerical simulations reveal that this scheme exhibits varying convergence rates for different choices of model parameters (i.e., domain size, final time T, and reaction coefficient κ). This elusive phenomenon is not unique to the L1 scheme but is also observed in other numerical methods for reaction-diffusion equations such as the backward Euler (IE) scheme, Crank-Nicolson (C-N) scheme, and two-step backward differentiation formula (BDF2) scheme. The existing literature lacks an explanation for the existence of two different convergence regimes, which has puzzled us for a long while and motivated us to study this inconsistency between the standard convergence theory and numerical experiences. In this paper, we provide a general methodology to systematically obtain error estimates that incorporate the exponential decaying feature of the solution. We term this novel error estimate the `decay-preserving error estimate' and apply it to the aforementioned IE, C-N, and BDF2 schemes. Our decay-preserving error estimate consists of a low-order term with an exponential coefficient and a high-order term with an algebraic coefficient, both of which depend on the model parameters. Our estimates reveal that the varying convergence rates are caused by a trade-off between these two components in different model parameter regimes. By considering the model parameters, we capture different states of the convergence rate that traditional error estimates fail to explain. This approach retains more properties of the continuous solution. We validate our analysis with numerical results.