Variable-Step-Size Linear Multistep Methods

• Variable-Step-Size Linear Multistep (VLM) methods are adaptive numerical integration schemes that modify stepsizes and coefficients based on historical ratios to solve ODEs and related problems.
• They dynamically compute method coefficients to satisfy order conditions and maintain zero-stability on nonuniform grids, ensuring robust convergence.
• VLM frameworks enable adaptive order and stepsize control, supporting applications from stiff ODEs and PDEs to optimization and inverse problem regularization.

A variable-step-size linear multistep method (VLM) refers to any family of $k$-step linear multistep methods applied to ODEs or related problems in which the grid $\{t_n\}$ is nonuniform, with the steps $h_n = t_{n+1} - t_n$ determined adaptively or by external criteria, and the method's coefficients $\alpha_j, \beta_j$ (and related operator structure) are explicit functions of the current and prior stepsizes or their ratios. This framework underlies the majority of modern ODE and PDE time integrators used in scientific computing, controls the accuracy-order/stability tradeoff, and forms the analytic basis for contemporary adaptive and embedded-variable-order schemes.

1. General Formulation and Core Properties

A $k$-step variable-step-size LMM for an initial value problem $y' = f(t, y)$ assumes the recurrence

$$\sum_{j=0}^k \alpha_{j,n} y_{n+j} = h_{n+k-1} \sum_{j=0}^k \beta_{j,n} f(t_{n+j}, y_{n+j}),$$

where $h_n = t_{n+1} - t_n$ are generally nonuniform and the grid $\{t_n\}$ is arbitrary or adaptively determined. The coefficients $\alpha_{j,n}, \beta_{j,n}$ become functions of stepsize ratios, such as $r_n = h_n / h_{n-1}$ and related sequences, breaking the translation invariance and Toeplitz structure of classical, uniform-grid theory.

Zero-stability—expressed via boundedness of solutions to the homogeneous recurrence $\sum_{j=0}^k \alpha_{j, n} e_{n+j} = 0$—is necessary for convergence and, in the variable-step context, is substantially more intricate, involving analysis of products of varying companion matrices rather than a single operator. For strongly zero-stable methods on uniform grids, bounded smoothness in the step-ratio sequence ($h_n/h_{n-1} = 1 + O(1/N)$ as $N \to \infty$) suffices to transfer zero-stability to the variable-step case (Söderlind et al., 2018).

The general VLM form extends to implicit, explicit, and linearly-implicit settings, to quadrature for Volterra equations, and to convex optimization interpretation as in the analysis of Nesterov acceleration (Nozawa et al., 2024).

2. Zero-Stability on Nonuniform Grids

Zero-stability for variable-step LMMs depends critically on the manner in which the grid is constructed. If the grid is the image of an equidistant grid under a $C^2$-diffeomorphism $\Phi: [0,1] \to [0,1]$ with $\Phi(0)=0, \Phi(1)=1$, and $\Phi' > 0$, then for strongly zero-stable schemes (i.e., all extraneous characteristic roots strictly inside the unit circle), there exists a finite threshold $N^*$ so that for all $N>N^*$, the variable-step LMM remains zero-stable (Söderlind et al., 2018).

This critical result is realized by analyzing the extraneous operator $R_N(\Phi)$, which is a perturbation of the constant-step Toeplitz operator $T_0$. The norm of the perturbation scales as $O(1/N)$, with explicit Neumann-series control. Thus,

$$\|R_N(\Phi)^{-1}\| \leq \frac{C_0}{1 - \|E_N\| C_0}, \quad \text{if } \|E_N\| C_0 < 1,$$

where $C_0$ is the uniform Toeplitz inverse bound and $E_N$ the perturbation norm.

The practical implication is that strong zero-stability is preserved under step sequences with $h_n/h_{n-1} = 1 + O(1/N)$, which is enforced in standard smooth-step controllers using digital filters or bounded-variation updates to $\log h$ (Söderlind et al., 2018).

3. Order Conditions, Consistency, and SSP Theory

Order conditions for variable-step LMMs generalize the classical approach by enforcing consistency equations in terms of cumulative sums of stepsize ratios. For an explicit $k$-step method, the coefficients depend on historical step ratios $\omega_{j,n} = h_{n-k+j}/h_n$ and their sums $\Omega_{j,n} = \sum_{i=1}^j \omega_{i,n}$. The $p$th-order conditions read (Hadjimichael et al., 2015): $$\sum_{j=0}^{k-1} \alpha_j = 1, \qquad \sum_{j=0}^{k-1} \left( \Omega_j^m \alpha_j + m \Omega_j^{m-1} \beta_j \right) = \Omega_k^m \quad (1 \le m \le p).$$ Strong Stability Preserving (SSP) variable-step LMMs require nonnegativity constraints on the combinations $\alpha_{j,n} - r \beta_{j,n} \ge 0$ and lead to an explicit optimal step-size control: $$h_n = \mathcal{C}_n \mu_n, \quad \mathcal{C}_n = \max \{ r \ge 0 : \alpha_{j,n} - r \beta_{j,n} \ge 0, \forall j \},$$ with $\mu_n$ a local problem-dependent bound such as the Forward Euler step restriction. Tight sharp upper bounds on SSP coefficients are given, and the existence of arbitrarily high-order methods for bounded step-ratios is established (Hadjimichael et al., 2015).

4. Adaptive Control of Stepsize and Order

Adaptive step-size and variable-order strategies are central in VLM practice. Embedded variable-step variable-order (VSVO) families, such as the MOOSE234 scheme, implement simultaneous solutions of multiple order (e.g., via a BDF3 solve and filtered linear combinations for order 2, 3, and 4) with negligible overhead. Step and order selection employ local error estimators computed from embedded lower- and higher-order solutions, accepting the order yielding the largest feasible new step: $$j^* = \arg\max_{j: |\mathrm{Est}_j| < \varepsilon} \left( \frac{\varepsilon}{|\mathrm{Est}_j|} \right)^{1/(j+1)}, \quad h_{n+1} = \gamma h_n \left( \frac{\varepsilon}{| \mathrm{Est}_{j^*}| } \right)^{1/(j^*+1)}.$$ This approach enables robust and computationally efficient adaptive methods that require only one nonlinear solve per step, yet provide competitive accuracy and efficiency with minimal added complexity (DeCaria et al., 2018).

5. Special Classes: Linearly-Implicit, SSP, and Optimization-VLMs

Variable-step linearly implicit multistep methods (LIMMs) are characterized by requiring only one linear system solve per step, formulated as

$$\sum_{i=-1}^{k-1} \alpha_i y_{n-i} = h_n \sum_{i=0}^{k-1} \beta_i f_{n-i} + h_n J_n \left( \sum_{i=-1}^{k-1} \mu_i y_{n-i} + h_n \sum_{i=0}^{k-1} \nu_i f_{n-i} \right),$$

where $J_n$ is the local Jacobian, and the method, including its order and stability maintenance, is adapted to variable steps (Glandon et al., 2020). LIMM, BDF, and explicit VLM classes are competitive, with LIMM offering lower per-step cost in many large-scale and stiff regimes.

In convex optimization, the Nesterov accelerated gradient method for $L$-smooth convex functions has been reinterpreted as an explicit two-step variable-step LMM, with linearly increasing steps $h_n$. Zero-stability and absolute stability are addressed via variable-step characteristic polynomials. NAG-c is proved optimal among consistent, stable explicit two-step VLMs with a certain Lyapunov-based $O(1/n^2)$ guarantee, and a broader class of coefficient-tuned VLMs provides superior performance in ill-conditioned settings (Nozawa et al., 2024).

6. Regularization via VLM for Volterra Integral Equations

In the context of first-kind Volterra integral equations with smooth kernels, VLM-based quadrature combines a linear multistep discretization on an equidistant mesh with a variable global stepsize $h$, determined a priori or via an adaptive balancing principle. Error bounds of the form

$$\max_n |u_n^\delta(h) - u(x_n)| = O(h^p + \delta/h),$$

lead to the optimal regularizing choice $h \sim \delta^{1/(p+1)}$ and error $O(\delta^{p/(p+1)})$. The balancing principle adaptively selects $h$ from a geometric mesh family, comparing reconstructions for stability and consistency without requiring knowledge of the solution smoothness (Plato, 2016).

7. Implementation Considerations and Practical Implications

The main practical implications of VLM theory are:

• Zero-stability on nonuniform grids critically depends on the smoothness of the stepsize sequence. Digital-filter-based step controllers ensure $h_n/h_{n-1} = 1 + O(1/N)$ without restricting the global step range (Söderlind et al., 2018).
• Embedded VSVO VLMs deliver adaptive order at negligible marginal cost and are especially effective for stiff ODEs and PDEs, as demonstrated in MOOSE234 implementations (DeCaria et al., 2018).
• Linearly implicit and explicit VLMs offer computational savings and robust stability in large-scale simulations, provided that suitable stability regions and error controls are enforced (Glandon et al., 2020).
• In optimization, recognizing accelerated gradient flows as variable-step LMMs expands the theoretical understanding of acceleration mechanisms and enables construction of new schemes adapted to spectral or conditioning considerations (Nozawa et al., 2024).
• In the regularization of inverse problems, step-size adaptation via balancing principles achieves optimal error rates for noisy data with modest computational overhead, even when underlying smoothness is unknown (Plato, 2016).

A plausible implication is that the mathematical flexibility of VLM frameworks enables seamless integration of modern step-size and order control, stability theory, efficiency optimization, and interpretive links to advanced algorithms across numerics and optimization.