Matrix-Free Stabilized BDF Schemes for Semilinear Parabolic Equations with Unconditional Maximum Bound Principle Preservation and Energy Stability

We develop a family of stabilized backward differentiation formula (sBDF) schemes of orders one through four for semilinear parabolic equations. The proposed methods are designed to achieve three properties that are rarely available simultaneously in high-order time discretizations: unconditional preservation of the maximum bound principle (MBP), unconditional discrete energy stability, and practical matrix-free implementation. The construction integrates carefully designed stabilization terms, fixed-point iterations, and a pointwise cut-off strategy. The nonlinear algebraic systems arising from the implicit sBDF discretizations are solved by fixed-point iteration, resulting in fully matrix-free algorithms. This makes the approach particularly attractive for practical computations on general domains and under mixed boundary conditions, where FFT-based exponential time differencing methods are often unavailable or inefficient. We further present a unified analysis for the fully implemented schemes, explicitly incorporating the interplay among time discretization, nonlinear iteration, and cut-off. Unconditional contractivity of the fixed-point iterations and error estimates are established. For the Allen-Cahn equation, we additionally prove an unconditional discrete energy dissipation law. Numerical experiments confirm the theoretical convergence rates and demonstrate the robustness and efficiency of the proposed methods, particularly relative to ETD-based approaches for problems with mixed boundary conditions.

💡 Research Summary

The paper addresses a long‑standing challenge in the numerical simulation of semilinear parabolic equations, namely the simultaneous preservation of three desirable properties in a high‑order time discretization: (i) unconditional adherence to the maximum bound principle (MBP), (ii) unconditional discrete energy stability, and (iii) a matrix‑free implementation that scales efficiently on general domains and under mixed boundary conditions.

Starting from the classical backward differentiation formula (BDF) of orders one through four, the authors introduce carefully calibrated stabilization terms of the form β_k B Δt ∇_k ϕ^{n+1}, where ∇_k denotes the k‑th order backward finite‑difference operator and β_k = {1,2,6,12} for k = 1,…,4. These terms are O(Δt^{k+1}) and therefore do not degrade the nominal k‑th order accuracy of the underlying BDF scheme, yet they render the associated fixed‑point mapping a strict contraction for any choice of time step Δt and mesh size h.

The nonlinear algebraic systems generated by the implicit sBDF_k schemes are solved by a simple fixed‑point iteration. Because each iteration only requires applying the discrete Laplacian and the linear stabilization term, no global matrix assembly or factorization is needed; the algorithm is fully matrix‑free and costs O(N) per iteration (N = number of spatial degrees of freedom). For the first‑order scheme (sBDF1) the fixed‑point map is monotone and already respects the MBP, so no additional treatment is required. For higher orders (k ≥ 2) a pointwise cut‑off operator Π_{