Variable-Step DLN Time Integrator

• The variable-step DLN integrator is a two-step time-stepping scheme that guarantees unconditional energy stability (G-stability) and second-order accuracy even with nonuniform time grids.
• It uses weighted combinations of past solution values and adaptive time-stepping strategies to effectively simulate stiff systems, fluid models, phase-field problems, and fourth-order active flows.
• Numerical tests validate its robust performance, efficiency in handling multi-scale dynamics, and seamless integration with mixed finite element spatial discretization.

The Variable-Step Dahlquist–Liniger–Nevanlinna (DLN) Time Integrator is a nonlinear, two-step time-stepping scheme for the numerical integration of ordinary and partial differential equations, characterized by rigorous energy stability (G-stability) and provable second-order accuracy under arbitrary, nonuniform time grids. Defined by a family of weighted combinations of multiple past solution values, the DLN integrator is widely used for stiff systems, fluid models, phase-field problems, and fourth-order active flows, often coupled with mixed finite element spatial discretization and adaptive time-stepping strategies that enhance efficiency without sacrificing long-term stability or accuracy.

1. Mathematical Formulation of the DLN Integrator

The DLN method is fundamentally a multistep, one-leg, two-step scheme for initial value problems of the form $y'(t) = f(t, y(t))$. At time levels $t_n$, with variable steps $k_n = t_{n+1} - t_n$, the DLN discretization forms a weighted combination of three consecutive values: $$\sum_{\ell=0}^2 \alpha_\ell y_{n-1+\ell} = \widehat{k}_n \cdot f\left(\sum_{\ell=0}^2 \beta_\ell^{(n)} t_{n-1+\ell}, \sum_{\ell=0}^2 \beta_\ell^{(n)} y_{n-1+\ell}\right),$$ where $(\alpha_0, \alpha_1, \alpha_2)$ and $(\beta_0^{(n)}, \beta_1^{(n)}, \beta_2^{(n)})$ are coefficients designed to maintain second-order accuracy and are functions of the time-step variability parameter $\varepsilon_n = (k_n - k_{n-1})/(k_n + k_{n-1})$ and a free parameter $\theta \in [0,1]$. The effective time-step is $\widehat{k}_n = \alpha_2 k_n - \alpha_0 k_{n-1}$. This approach generalizes backward differentiation formulas while guaranteeing strong nonlinear stability properties unachievable in BDF2 under nonuniform grids (Layton et al., 2020).

In the context of partial differential equations (PDEs), especially those with mixed finite element spatial discretization, the DLN temporal update for a discrete solution $u^h$ can be stated as: $$\frac{1}{\widehat{k}_n}(u_{n,\alpha}^h, v^h) + \mathcal{A}(u_{n,\beta}^h, v^h) + \cdots = \text{(forcing terms)},$$ where $\mathcal{A}$ is typically a bilinear operator (e.g., Laplacian or convection-diffusion operator), and the evaluation points $u_{n,\alpha}^h$ and $u_{n,\beta}^h$ are defined by: $$u_{n,\alpha}^h = \alpha_2 u_{n+1}^h + \alpha_1 u_n^h + \alpha_0 u_{n-1}^h,\quad u_{n,\beta}^h = \beta_2^{(n)} u_{n+1}^h + \beta_1^{(n)} u_n^h + \beta_0^{(n)} u_{n-1}^h.$$

2. Stability Theory: G-Stability and Energy Dissipation

The cornerstone of the DLN method is its unconditional nonlinear stability, established via G-stability theory. For any time-step sequence, the DLN scheme preserves a discrete energy estimate analogous to the physical energy dissipation of the underlying model. For discrete solution vectors $u_n$, the method obeys: $$\|[u_{n+1}, u_n]^T\|^2_{G(\theta)} - \|[u_n, u_{n-1}]^T\|^2_{G(\theta)} + \|\sum_{\ell=0}^2 a_\ell^{(n)} u_{n-1+\ell}\|^2 = \text{(dissipation terms)},$$ where the G-norm is defined as $$\|[v, w]\|^2_{G(\theta)} = \frac{1}{2}[(1+\theta)\|v\|^2 + (1-\theta)\|w\|^2]$$ and the $a_\ell^{(n)}$ are scheme-specific weights. This ensures strict control over numerical energy regardless of time-step variability (Layton et al., 2020, Qin et al., 2020, Siddiqua et al., 2023, Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).

Applied to PDEs, such as the Navier–Stokes equations or phase-field models, the DLN method maintains unconditional stability of the discrete kinetic or model energy,

$$E(t_n) = \frac{1}{2}(\|u(t_n)\|^2 + \|p(t_n)\|^2 + \cdots),$$

for all $n$, independent of step size or nonlinear model coupling. This property is particularly significant for stiff problems and long-term simulations.

3. Error Analysis and Convergence Rates

The DLN integrator is second-order accurate in time even on arbitrary, nonuniform grids, as proven in multiple works (Layton et al., 2020, Qin et al., 2020, Pei, 2023, Pei, 2024, Chen et al., 2024, Zheng et al., 19 Oct 2025). The local truncation error for solutions $u$ of sufficient smoothness is

$$|(a_2 u(t_{n+1}) + a_1 u(t_n) + a_0 u(t_{n-1})) - u(t_{n,\beta})| \leq C (k_n + k_{n-1})^3 \|u_{tt}\|,$$

with the overall error bound: $$\max_n \|u(t_n) - u_n\| \leq C [h^{r+1} + k_{\max}^2].$$ Here $h$ denotes the spatial mesh size and $k_{\max}$ the largest time-step. For semi-discrete PDE solutions (with FEM in space), analogous estimates hold for velocity, pressure, auxiliary variables, and phase-field variables (Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).

4. Adaptive Time-Stepping Strategies

DLN’s native compatibility with nonuniform grids permits adaptive time-stepping mechanisms that optimize computational efficiency. Two widely-adopted criteria are:

• Local truncation error estimators: The time-step $k_{n+1}$ is adjusted by comparing the DLN update against an independent second-order estimate (e.g., AB2-type), ensuring the error remains below a prescribed tolerance (Pei, 2023, Pei, 2024, Chen et al., 2024).
• Minimum-dissipation criterion: The ratio of numerical dissipation (ND) to physical (viscous or model) dissipation (VD, PD) is monitored:

$$\chi_u = \frac{\varepsilon_{\text{ND}}^u}{\varepsilon_{\text{VD}}^u},\quad \chi_\phi = \frac{\varepsilon_{\text{ND}}^\phi}{\varepsilon_{\text{PD}}^\phi},$$

where $$\varepsilon_{\text{ND}}^u = (1/\widehat{k}_n)\|u_{n,\alpha}^h\|^2$$, $$\varepsilon_{\text{VD}}^u = \mu\|\nabla u_{n,\beta}^h\|^2$$, etc. If $\max(\chi_u, \chi_\phi)$ is less than a tolerance $\delta$, the step is increased; if larger, it is reduced (Pei, 2023, Siddiqua et al., 2023, Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).

Adaptive DLN algorithms consistently reduce the number of time steps in stiff and multiscale problems while retaining stability and accuracy.

5. Application Domains

DLN variable-step integrators have broad applicability across computational science:

• Fluid Dynamics: Proven for incompressible Navier–Stokes and Stokes/Darcy flows (Layton et al., 2020, Qin et al., 2020, Pei, 2023, Pei, 2024), DLN schemes offer robust energy stability, second-order convergence, and superior performance over BDF2 when the solution involves rapid transitions or turbulence. Ensemble DLN methods can efficiently handle multiple simultaneous NSE systems (Pei, 2024).
• Phase-field Models: DLN schemes have been employed in Allen–Cahn equations (Chen et al., 2024) and in coupled Allen–Cahn active fluid problems (Zheng et al., 19 Oct 2025), delivering unconditional stability for stiff nonlinearities and gradient flows.
• Fourth-order Active Fluid Models: Mixed FEM discretizations, coupled with variable-step DLN, have facilitated efficient and accurate simulation of models involving higher-order operators and complex nonlinearities (Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).
• Turbulence and Smagorinsky Models: DLN integrators combined with corrected Smagorinsky closures manage adaptive dissipation and capture energy backscatter (Siddiqua et al., 2023).

6. Implementation: Refactorization and Filtering

Direct implementation of DLN can be intricate due to variable-step-dependent coefficients. Refactorization simplifies deployment, especially in legacy codes (Layton et al., 2021, Pei, 2023, Pei, 2024, Chen et al., 2024). The procedure consists of:

1. Pre-processing: Compute a weighted combination of previous solution approximations, $y_{\text{old}} = a_1 y_n + a_2 y_{n-1}$.
2. Backward Euler Step: Solve $$y_{\text{new}} - y_{\text{old}} = \Delta t_{BE} f(z_{\text{new}}, y_{\text{new}})$$, with modified stepsize.
3. Post-processing: Recover the DLN solution with $$y_{n+1} = c_2 y_{\text{new}} + c_1 y_n + c_0 y_{n-1}$$.

This method preserves the full variable-step DLN accuracy and G-stability properties while integrating efficiently with backward Euler infrastructure.

7. Numerical Verification and Practical Impact

A large body of numerical experiments confirms the theoretical properties of DLN integrators:

• Convergence tests consistently report second-order accuracy in time and optimal spatial orders for velocity, pressure, auxiliary variables, and phase fields (Layton et al., 2020, Qin et al., 2020, Pei, 2023, Pei, 2024, Chen et al., 2024, Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).
• Long-term energy stability is observed even in highly turbulent and chaotic regimes, as well as stiff dissipative systems.
• Adaptive DLN methods substantially reduce computational effort for problems with widely varying timescales, without loss of accuracy or stability (Pei, 2023, Pei, 2024, Chen et al., 2024, Zheng et al., 28 Jul 2025, Zheng et al., 23 Sep 2025, Zheng et al., 19 Oct 2025).

Summarized numerical evidence supports the DLN integrator as a robust, efficient solution for stiff, multiscale, and adaptive time-dependent PDEs.

---

In conclusion, the Variable-Step DLN Time Integrator is a rigorously stable, second-order accurate multistep scheme for nonlinear ODEs and PDEs on nonuniform time grids. Its architecture—grounded in Dahlquist theory and realized efficiently through refactorization—yields adaptability, energy conservation, and accuracy across a wide range of complex applications. Adaptive control of numerical dissipation further enhances its efficiency for modern computational simulations.