Variable-Step DLN Time Integrator

• The DLN integrator is a one-leg, two-step method providing unconditional G-stability and second-order accuracy for variable or adaptive time steps.
• It uses an exact algebraic refactorization akin to backward Euler with linear filtering, reducing computational overhead and coupling seamlessly with mixed FEM.
• Numerical studies demonstrate its efficiency and controlled energy dissipation in applications such as incompressible flows, phase-field models, and nonlinear active fluids.

The variable-step Dahlquist–Liniger–Nevanlinna (DLN) time integrator is a one-leg, two-step multistep method designed for unconditionally stable and second-order accurate time discretization of ordinary and partial differential equations with variable or adaptive time steps. The method is distinguished by its uniform $G$-stability, robust handling of non-uniform time grids, and ability to avoid nonphysical energy growth even under aggressive adaptivity. High-performance applications include incompressible flow, phase-field models, and nonlinear active fluid systems, especially when coupled with divergence-free mixed finite element spatial discretizations or minimum-dissipation adaptivity strategies (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Zheng et al., 28 Jul 2025, Chen et al., 28 Sep 2024, Pei, 2023, Pei, 26 Jul 2024, Qin et al., 2020, Layton et al., 2020, Layton et al., 2021).

1. General Formulation and Coefficient Structure

The DLN time integrator discretizes the ODE system $y'(t) = f(t, y(t))$ using an arbitrary, possibly adaptive partition $0 = t_0 < t_1 < \cdots < t_M = T$. Time steps are denoted $k_n = t_n - t_{n-1}$, with the step-ratio parameter

$$\epsilon_n = \frac{k_n - k_{n-1}}{k_n + k_{n-1}},\quad r_n = \frac{k_n}{k_{n-1}},\quad \epsilon_n = \frac{r_n - 1}{r_n + 1}.$$

The method uses the recurrence

$$\alpha_0 y_{n-1} + \alpha_1 y_n + \alpha_2 y_{n+1} = \widehat k_n\, f(t_{n,\beta}, y_{n,\beta}),$$

where the coefficients depend only on a user-prescribed parameter $\theta \in [0,1]$: $$\alpha_2 = \frac{1+\theta}{2},\quad \alpha_1 = -\theta,\quad \alpha_0 = \frac{\theta-1}{2},\qquad \widehat k_n = \alpha_2 k_n - \alpha_0 k_{n-1}.$$ The stage values $(t_{n,\beta}, y_{n,\beta})$ are convex combinations of $(t_{n-1}, t_n, t_{n+1})$ and $(y_{n-1}, y_n, y_{n+1})$ with weights $(\beta_0^{(n)}, \beta_1^{(n)}, \beta_2^{(n)})$ determined by $\theta$ and $\epsilon_n$. These guarantee second-order accuracy and $G$-stability for arbitrary step sequences. Uniform time steps recover classical schemes: $\theta=1/2$ yields BDF2, $\theta=1$ the midpoint rule, and $\theta=0$ a two-step midpoint scheme (Zheng et al., 19 Oct 2025, Chen et al., 28 Sep 2024, Layton et al., 2021).

2. Nonlinear $G$-Stability and Energy-Dissipation Properties

The nonlinear and unconditional $G$-stability of DLN is central to its robust behavior on variable grids. Defining the $G$-norm for a pair $(u, v)$ by

$$\| (u, v) \|^2_{G(\theta)} = \tfrac{1}{4}(1+\theta)\|u\|^2 + \tfrac{1}{4}(1-\theta)\|v\|^2,$$

the method satisfies a telescoping energy identity: $$(y_{n,\alpha}, y_{n,\beta}) = \| (y_{n+1}, y_n)\|_{G(\theta)}^2 - \| (y_n, y_{n-1})\|_{G(\theta)}^2 + \|\sum_{\ell=0}^2 a_\ell^{(n)} y_{n-1+\ell} \|^2,$$ for explicit $a_\ell^{(n)}$. For nonlinear dissipative problems, this structure yields a monotonic decay of the discrete energy, preventing nonphysical growth due to variation in step size, with or without viscosity (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Pei, 26 Jul 2024, Layton et al., 2020, Layton et al., 2021). In the linear test equation $y' = \lambda y$ with $\mathrm{Re}\,\lambda \le 0$, $G$-stability yields unconditional $A$-stability.

3. Consistency, Local Truncation Error, and Accuracy

The DLN method is globally second-order accurate for any time-grid with step ratios bounded away from zero. A Taylor expansion of the solution reveals that the local truncation error at each step $n$ is

$$\tau_{n+1} = C(\theta, r_n)\,k_n^2\, y'''(\xi) + O(k_n^3),$$

with $C(\theta, r_n)$ given explicitly, and $\xi \in (t_{n-1}, t_{n+1})$ (Zheng et al., 19 Oct 2025, Chen et al., 28 Sep 2024, Layton et al., 2021). Spatial discretizations, such as conforming FEM assemblages, combine to yield full error bounds of

$\|u_n^h - u(t_n)\| = O(k_{\max}^2 + h^{s+1}),$

under standard regularity assumptions and mild step-size moderation.

4. Refactorization and Efficient Implementation

Direct deployment of DLN requires managing non-trivial stage coefficients, but it admits exact algebraic refactorization as backward Euler plus linear time filtering. Each DLN step is realized as:

1. Linear combination of previous solutions (pre-filter),
2. Single backward Euler solve with effective time step and collocation,
3. Linear combination of new and old solutions (post-filter).

This structure dramatically reduces code complexity and enables embedding DLN into backward-Euler codebases with only minor augmentation (Pei, 2023, Layton et al., 2021, Chen et al., 28 Sep 2024). All history dependence is contained in short recurrences, and storage requirements are minimal (two solution vectors per unknown plus filters).

5. Adaptive Time-Stepping and Minimum-Dissipation Criterion

A distinctive feature of modern DLN practice is robust adaptive control, notably via the minimum-dissipation criterion. At each step, the numerical dissipation (ND) and physical dissipation (PD; often associated to viscosity) are evaluated: $$\mathrm{ND} = \| w_{n,\alpha} \|^2 / \widehat k_n, \qquad \mathrm{PD} = \mu \| \nabla w_{n,\beta} \|^2.$$ The ratio $\chi = \mathrm{ND}/\mathrm{PD}$ is compared to a tolerance $\delta$, and the step updated via: $$k_{n+1} = \begin{cases} \min(2k_n, k_{\max}), & \chi \le \delta \ \max(\tfrac12 k_n, k_{\min}), & \chi > \delta \end{cases}$$ (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Zheng et al., 28 Jul 2025, Pei, 26 Jul 2024, Layton et al., 2020). This ensures controlled energy dissipation and efficiently tracks dynamic stiffness. Alternative adaptivity strategies deploy embedded AB2-type local error estimation, but dissipation-based adaptivity aligns more directly with the $G$-stability framework for complex nonlinear systems.

6. Coupling With Mixed FEM and Application to Nonlinear Fourth-Order Systems

DLN integrators couple naturally with divergence-free preserving mixed FEM for nonlinear PDEs, such as the Allen–Cahn active fluid model or fourth-order active/dissipative flows. High-order operators (e.g., biharmonic) are handled via variable splitting (e.g., $w = \Delta u$), lowering the regularity requirements. Auxiliary, divergence-free variables (e.g., $\phi$ as Lagrange multiplier) and monolithic solvers are deployed at each time step (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Zheng et al., 28 Jul 2025). Numerical experiments confirm that minimum-dissipation-adaptive DLN can yield dramatic reductions in total time steps (e.g., 500 instead of 20,000 for comparable solution accuracy at high Reynolds numbers) and remain stable for complex long-time nonlinear dynamics.

7. Practical Recommendations and Observations From Numerical Studies

Key practical guidelines include:

• Limiting the step-ratio $r_n \in [0.5, 2.0]$ stabilizes coefficient conditioning.
• Choices of $\theta$ influence damping and error: $\theta \in (0.3, 0.7)$ often yields the best compromise; $\theta=1/2$ matches BDF2 under uniform steps.
• Startup values require two initial states; Crank–Nicolson initialization is commonly used.
• Implementation overhead is minimal compared to implicit Euler, especially under refactorization.
• Tests on Navier–Stokes, Allen–Cahn, Stokes/Darcy, and active fluid models confirm unconditional nonlinear stability, second-order temporal accuracy, and significant efficiency gains under adaptivity, compared to BDF2 or trapezoidal integrators, especially in highly stiff or multiscale regimes (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Zheng et al., 28 Jul 2025, Chen et al., 28 Sep 2024, Layton et al., 2020, Layton et al., 2021).

---

Summary Table: DLN One-Leg Two-Step Variable-Step Scheme

Step	Formula	Notes
Recurrence	$$\sum_{\ell=0}^2 \alpha_\ell y_{n-1+\ell} = \widehat k_n f(t_{n,\beta},y_{n,\beta})$$	$\alpha$, $\widehat k_n$, $\beta$ as above
Coefficient Form	$$\alpha_2 = \tfrac12(1+\theta),\ \alpha_1 = -\theta,\ \alpha_0 = \tfrac12(\theta-1)$$	$\theta\in[0,1]$
Adaptivity	If $\chi\le\delta,\ k_{n+1}=\min(2 k_n, k_{\max})$, else $k_{n+1}=\max(\tfrac12 k_n, k_{\min})$	Based on ND/PD ratio

This synthesis is directly grounded in the key recent and foundational references on the DLN integrator family (Zheng et al., 19 Oct 2025, Zheng et al., 23 Sep 2025, Zheng et al., 28 Jul 2025, Chen et al., 28 Sep 2024, Pei, 2023, Pei, 26 Jul 2024, Qin et al., 2020, Layton et al., 2020, Layton et al., 2021).