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

TL;DR

The paper tackles the puzzling observation that time-stepping schemes for diffusion and subdiffusion with weak initial singularities exhibit parameter-dependent convergence orders. It introduces decay-preserving error estimates that couple an exponential-in-time component with a higher-order algebraic component, deriving explicit bounds for implicit Euler, Crank-Nicolson, and BDF2 on ODEs, and extending them to PDEs via an eigenvalue framework. The results reveal two competing scales controlled by model parameters (e.g., κ, λ_1, and the final time $T$) that determine whether convergence at the final time is governed by $O(τ^{α})$ or $O(τ^{k})$ (with $k=1$ or $2$). The work further discusses a conjectured decay-preserving estimate for the L1 scheme in subdiffusion and supports it with numerical evidence, thereby offering a more faithful reflection of the continuous problem’s behavior and guiding practical time-step selection.

Abstract

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.

Figures (1)

• Figure 1: Errors with fixed model parameters $T, \kappa, \lambda_1$: absolute error for ODEs and PDEs.

Theorems & Definitions (17)

• Lemma 2.1
• Remark 1
• Lemma 2.2
• Theorem 2.1
• Lemma 2.3
• Theorem 2.2
• Theorem 2.3
• Theorem 3.1
• Theorem 3.2
• Theorem 3.3