Variable Time-Step Linear Relaxation Schemes

• Variable time-step linear relaxation schemes are nonuniform-in-time discretizations that handle nonlinear terms by linearizing or using auxiliary-variable updates for efficient linear solves.
• They preserve structural properties such as energy dissipation and stability through innovative design of discrete convolution kernels and auxiliary variable recalibration.
• Adaptive time-stepping in these schemes enhances simulation efficiency while maintaining second-order accuracy in complex time-fractional phase-field and related models.

Variable time-step linear relaxation schemes are nonuniform-in-time discretizations in which nonlinear terms are handled by linearization or by an auxiliary-variable relaxation so that each time level is obtained from linear algebraic solves rather than nonlinear solves. In the time-fractional phase-field setting, a representative construction uses the $L1^{+}$-CN formula to discretize the Caputo derivative and introduces an auxiliary variable to approximate the nonlinear term by directly solving algebraic equations rather than differential-algebraic equations as in the invariant energy quadratization (IEQ) and the scalar auxiliary variable (SAV) approaches. The resulting semi-discrete scheme is second-order accurate in time, and the discrete energy dissipation law is asymptotically compatible with the classical one when the fractional-order parameter $\alpha \rightarrow 1^{-}$ (Yu et al., 3 Sep 2025).

1. Conceptual setting and related formulations

The phrase linear relaxation scheme appears in several adjacent numerical traditions. In hyperbolic conservation laws, the relaxation framework replaces a nonlinear system by a nearby linear, strongly hyperbolic system plus a stiff source, with variable relaxed schemes choosing local subcharacteristic speeds based on local estimates of the wave speeds in order to reduce numerical diffusion relative to constant-relaxation formulations (Krishnamurthy et al., 2013). In nonlinear parabolic problems, a linearized variable-time-step BDF2 scheme combines variable-step BDF2 for the linear term with a Newton linearized method for the nonlinear term in time and a Galerkin finite element method in space, yielding a method in which each step requires only one linear solve (Zhao et al., 2022).

A second lineage is provided by relaxation-type Crank–Nicolson methods. For the Schrödinger–Poisson system, a Besse-style relaxation Crank–Nicolson scheme uses an auxiliary variable $\phi \approx |u|^2-\mu$, is explicit with respect to the nonlinearity, requires the solution of a linear system for each time-step, and remains compatible with variable time-steps $k_n=t_n-t_{n-1}$ (Athanassoulis et al., 2021). For stiff dissipative problems, the Dahlquist–Liniger–Nevanlinna (DLN) one-leg, two-step scheme provides unconditional $G$-stability for variable time-steps and second-order accuracy, and it can be refactorized as backward Euler with pre- and post arithmetic steps added (Layton et al., 2021). In the Allen–Cahn setting, variable time-stepping DLN schemes have been combined with a partially implicit modified algorithm and a scalar auxiliary variable algorithm, with finite element spatial discretization and proofs of unconditional, long-term stability of the model energy under any arbitrary time step sequence (Chen et al., 2024).

These formulations do not define a single unique algorithm. Rather, they delimit a class of methods characterized by three recurring features: nonuniform temporal meshes, linear treatment of the nonlinear contribution, and analysis framed around discrete structural properties such as stability, energy dissipation, or conservation.

2. Core construction on nonuniform meshes

A standard starting point is a nonuniform partition

$0=t_0<t_1<\cdots<t_N=T,\qquad \tau_k=t_k-t_{k-1},$

together with adjacent step ratios such as $r_k=\tau_k/\tau_{k-1}$. In variable-step BDF2 formulations for nonlinear parabolic equations, the derivative approximation is

$$\mathcal{D}_2 U^n = \frac{1+2r_n}{\tau_n(1+r_n)} (U^n-U^{n-1}) - \frac{r_n^2}{\tau_n(1+r_n)} (U^{n-1}-U^{n-2}),$$

and the nonlinearity is handled by Newton linearization,

$f(u^n)\approx f(U^{n-1})+f'(U^{n-1})(U^n-U^{n-1}),$

so that the semi-discrete scheme becomes

$$\mathcal{D}_2 U^n = \Delta U^n+f(U^{n-1})+f'(U^{n-1})(U^n-U^{n-1}) .$$

Its fully discrete Galerkin FEM version seeks $\alpha \rightarrow 1^{-}$0 such that, for all $\alpha \rightarrow 1^{-}$1,

$\alpha \rightarrow 1^{-}$2

The first time step is computed by the BDF1/backward Euler scheme (Zhao et al., 2022).

In the time-fractional phase-field context, the Caputo derivative is discretized by the $\alpha \rightarrow 1^{-}$3-CN formula,

$\alpha \rightarrow 1^{-}$4

with

$\alpha \rightarrow 1^{-}$5

The defining relaxation step is the auxiliary-variable update. For a quartic potential,

$\alpha \rightarrow 1^{-}$6

the auxiliary variable is introduced by

$\alpha \rightarrow 1^{-}$7

where $\alpha \rightarrow 1^{-}$8 is a stabilization parameter. The key distinction from IEQ and SAV is that the auxiliary variable is recomputed from an algebraic equation at each new time level, rather than evolved by a differential equation (Yu et al., 3 Sep 2025).

For the time-fractional volume-conserved Allen–Cahn model, the resulting variable time-step scheme takes the form

$\alpha \rightarrow 1^{-}$9

$\phi \approx |u|^2-\mu$0

$\phi \approx |u|^2-\mu$1

Analogous constructions are given for the time-fractional Cahn–Hilliard equation and the time-fractional Swift–Hohenberg equation (Yu et al., 3 Sep 2025).

3. Discrete kernels, step-ratio constraints, and variable-step analysis

Variable time-stepping complicates the convolution structure that is available on uniform meshes. A major analytical device is therefore the introduction of auxiliary kernel systems. For the variable-step BDF2 method, the discrete orthogonal convolution (DOC) kernels are defined by

$\phi \approx |u|^2-\mu$2

and satisfy

$\phi \approx |u|^2-\mu$3

The discrete complementary convolution (DCC) kernels are defined by

$\phi \approx |u|^2-\mu$4

with the identity

$\phi \approx |u|^2-\mu$5

These kernels facilitate summation-by-parts, energy arguments, and Grönwall inequalities under variable time-steps, and their positivity and summability under a step-ratio restriction are central to the proof of stability and error bounds (Zhao et al., 2022).

Related nonuniform convolution machinery appears in variable-step L1 schemes for time-fractional equations. For the time-fractional Allen–Cahn equation, the L1 weights

$\phi \approx |u|^2-\mu$6

are used together with DOC kernels $\phi \approx |u|^2-\mu$7 and DCC kernels

$\phi \approx |u|^2-\mu$8

which are the discrete analogues of the fractional integral kernel and are used in the energy functional and in the discrete fractional Grönwall lemma (Liao et al., 2021).

A recurrent misconception is that there is a single universal admissible step-ratio bound for variable-step schemes. The published bounds are method-dependent and norm-dependent. For the linearized variable-time-step BDF2 scheme for nonlinear parabolic equations, the unconditionally optimal error estimate is proved under

$\phi \approx |u|^2-\mu$9

and a mild maximum step size condition

$k_n=t_n-t_{n-1}$0

with no CFL restriction (Zhao et al., 2022). For variable step-size BDF2 in the $k_n=t_n-t_{n-1}$1 setting for linear and semilinear parabolic equations, the sharp zero-stability bound is

$k_n=t_n-t_{n-1}$2

while $k_n=t_n-t_{n-1}$3 and $k_n=t_n-t_{n-1}$4 stability are established under the more relaxed bound

$k_n=t_n-t_{n-1}$5

(Wang et al., 2020). For the time-fractional Cahn–Hilliard equation with a refined L2-type approximation, the step-ratio restriction is relaxed to

$k_n=t_n-t_{n-1}$6

and a specialized nonuniform mesh is constructed so that

$k_n=t_n-t_{n-1}$7

(Li et al., 24 Aug 2025). By contrast, the DLN framework is explicitly described as supporting arbitrary variation in time step and, in the Allen–Cahn application, the model energy is shown to be stable under any arbitrary time step sequence (Layton et al., 2021, Chen et al., 2024).

4. Energy laws, conservation, and auxiliary-variable alignment

Structure preservation is one of the principal reasons variable time-step linear relaxation schemes are studied. In the 2025 time-fractional phase-field formulation, each model satisfies a discrete modified energy law,

$k_n=t_n-t_{n-1}$8

and, under mild step ratio constraints, a discrete variational dissipation law. For the Allen–Cahn model,

$k_n=t_n-t_{n-1}$9

As $G$0, this law converges to the classical second-order Crank–Nicolson dissipation law, which is the precise sense in which the scheme is asymptotically compatible with the classical one (Yu et al., 3 Sep 2025).

The same structural theme appears in earlier variable-step fractional schemes. For the time-fractional Allen–Cahn equation, the discrete variational energy is

$G$1

and the scheme satisfies

$G$2

on arbitrary nonuniform meshes (Liao et al., 2021). For the time fractional Swift–Hohenberg model, the modified discrete energy

$G$3

is non-increasing at each step, and as $G$4 the law reduces to the classical case (Zhao et al., 2023).

Conservation properties are likewise retained in several formulations. The refined variable-step L2 scheme for the time-fractional Cahn–Hilliard equation is proved to have exact discrete volume conservation and proper energy dissipation laws (Li et al., 24 Aug 2025). The relaxation Crank–Nicolson method for Schrödinger–Poisson preserves a discrete mass law exactly even for variable time-steps, while its discrete energy balance acquires a remainder term proportional to the change in time-step; when the time-steps are constant, the extra term vanishes (Athanassoulis et al., 2021).

A central point of the algebraic relaxation strategy is the auxiliary-variable error. In the 2025 phase-field paper, numerical results confirm that the auxiliary variable remains well aligned with the original variable, and the error between them does not continue to increase over time before the system reaches steady state. The explicit contrast drawn in that work is that, in IEQ and SAV, the error between the auxiliary variable and its defining nonlinear function can accumulate continuously in long simulations because the auxiliary variable is propagated dynamically rather than corrected at each step (Yu et al., 3 Sep 2025).

5. Accuracy, initialization, and implementation pathways

The principal accuracy claims are high-order but not uniform across all formulations. For the variable-time-step BDF2 discretization of nonlinear parabolic equations, the main result is the unconditional optimal estimate

$G$5

and the analysis also shows that the first level solution $G$6 obtained by BDF1 does not cause the loss of global accuracy of second order (Zhao et al., 2022). For the time-fractional phase-field relaxation scheme based on $G$7-CN, the semi-discrete scheme is second-order accurate in time (Yu et al., 3 Sep 2025). For the Schrödinger–Poisson relaxation Crank–Nicolson method, the scheme is second order in time and explicit with respect to the nonlinearity (Athanassoulis et al., 2021).

The interaction between startup procedures and global order is more delicate in BDF-type methods than a schematic description sometimes suggests. In the analysis of variable step-size BDF2 for linear and semilinear parabolic equations, a trapezoidal start retains second-order accuracy in all norms, whereas a backward Euler start leads to order reduction in several norms unless the initial step is chosen sufficiently small, for example $G$8 (Wang et al., 2020). This should be read together with the different BDF2 initialization result above: the conclusions depend on the exact scheme, norm, and proof framework.

Implementation complexity is also method-specific. The DLN scheme has been underutilized due, partially, to its complexity of direct implementation, but its refactorization shows that it is equivalent to the backward Euler method with pre- and post arithmetic steps added, making it possible to add DLN capability with only scalar algebraic pre/post steps wrapped around an existing backward Euler solver (Layton et al., 2021). In the Allen–Cahn DLN-SAV algorithm, the refactorization process similarly reduces the main solve to standard linear solves, with only a few steps added to a backward Euler SAV code (Chen et al., 2024).

Scheme family	Linearization or relaxation device	Reported property
Variable-step BDF2 for nonlinear parabolic equations	Newton linearization	Unconditionally optimal $G$9 error estimate (Zhao et al., 2022)
DLN and DLN2BE	Pre-/post-processed backward Euler refactorization	Unconditional $0=t_0<t_1<\cdots<t_N=T,\qquad \tau_k=t_k-t_{k-1},$0-stability for variable time-steps (Layton et al., 2021)
Variable-step DLN for Allen–Cahn	Partially implicit modified algorithm or SAV algorithm	Unconditional, long-term stability under any arbitrary time step sequence (Chen et al., 2024)
Relaxation Crank–Nicolson for Schrödinger–Poisson	Auxiliary variable $0=t_0<t_1<\cdots<t_N=T,\qquad \tau_k=t_k-t_{k-1},$1	Linear system for each time-step; exact mass preservation (Athanassoulis et al., 2021)
$0=t_0<t_1<\cdots<t_N=T,\qquad \tau_k=t_k-t_{k-1},$2-CN linear relaxation for fractional phase fields	Algebraic auxiliary-variable update	Second-order accuracy; asymptotically compatible energy law (Yu et al., 3 Sep 2025)

The common implementation consequence is straightforward: once the nonlinearity is shifted to extrapolation, Newton linearization, or an algebraic auxiliary-variable relation, the time step is organized around linear algebraic subproblems rather than nonlinear coupled solves.

6. Applications, adaptivity, and numerical behavior

Variable time-step linear relaxation schemes are used across dissipative, dispersive, and conservation-law models. In the 2025 phase-field setting, the time-fractional volume-conserved Allen–Cahn equation, the time-fractional Cahn–Hilliard equation, and the time-fractional Swift–Hohenberg equation serve as model problems for the algebraic relaxation strategy (Yu et al., 3 Sep 2025). In the Allen–Cahn DLN study, one- and two-dimensional numerical tests are used to verify unconditional non-linear stability, long-term energy behavior, and the performance of time-adaptive algorithms (Chen et al., 2024). For the time fractional Swift–Hohenberg equation, adaptive time-stepping is used to capture the multi-scale time behavior efficiently, and the reported tests show about 5–20 times fewer steps in some cases while still preserving energy decay (Zhao et al., 2023).

The adaptive mechanisms are typically solution-driven. In the 2025 phase-field relaxation work, the control law is

$0=t_0<t_1<\cdots<t_N=T,\qquad \tau_k=t_k-t_{k-1},$3

so that rapid transients trigger refinement and slow evolution permits coarsening (Yu et al., 3 Sep 2025). Closely related feedback formulas are used in variable-step L1 methods for the Allen–Cahn and Swift–Hohenberg equations, where they are described as tracking multi-scale behaviors and drastically reducing computational cost in slow-evolving regimes (Liao et al., 2021, Zhao et al., 2023). This suggests that variable time-stepping is not merely a consistency device for initial weak singularities; it is also a practical mechanism for long-time simulation in memory-driven gradient flows.

Outside phase-field models, the numerical behavior is similarly tied to local adaptation. The variable relaxed schemes for multidimensional hyperbolic conservation laws use local maximum and minimum speeds, are free of nonlinear Riemann solvers and independent of the underlying eigenstructure, and produce sharper fronts than constant relaxation in the weakly hyperbolic gas injection displacement examples discussed in the source material (Krishnamurthy et al., 2013). For Schrödinger–Poisson dynamics, the variable time-step relaxation Crank–Nicolson scheme remains robust in stiff regimes and with highly oscillatory solutions, and the energy deviation induced by time-step changes is described as predictable and small (Athanassoulis et al., 2021).

The broader significance is therefore methodological rather than tied to a single equation. Across these works, variable time-step linear relaxation schemes are used to reconcile nonuniform temporal resolution with linear per-step solvers, while preserving, in the discrete setting, the qualitative structures that govern the underlying continuous models.