Adaptive Pseudo-Time Stepping

• Adaptive pseudo-time stepping is a method that advances an evolution variable using nonuniform increments, thereby controlling stability, accuracy, and computational cost.
• It applies across diverse areas, from residual neural networks and Hamiltonian systems to stiff ODEs, SPDEs, and energy-dissipative gradient flows.
• The approach leverages adaptive controllers, state-based reparametrizations, and error-indicator strategies to preserve structure and optimize simulation efficiency.

Searching arXiv for recent and foundational papers on adaptive pseudo-time stepping and closely related adaptive time-stepping formulations. Adaptive pseudo-time stepping denotes a class of numerical and algorithmic strategies in which an evolution variable that is not necessarily physical time is advanced with nonuniform increments chosen to balance stability, robustness, accuracy, and computational cost. Across the literature, the “pseudo-time” interpretation appears in several distinct but structurally related settings: residual networks viewed as explicit ODE discretizations along depth, geometric time reparametrizations of Hamiltonian or particle systems, adaptive multistep integration for dissipative gradient flows, explicit stabilization strategies for stiff ODEs and SPDEs, and goal-oriented timestep control for weakly instationary hyperbolic PDEs. A unifying theme is that the step size is treated as a control variable for the evolution itself rather than as a passive discretization parameter. In residual networks, this control acts along network depth and is learned jointly with model parameters (Yang et al., 2019). In geometric integration, it appears through a reparametrization $dt/d\tau=g(z)$, so that a constant step in pseudo-time $\tau$ induces a variable step in physical time $t$ (Richardson et al., 2011, Toggweiler et al., 2012). In PDEs and rate-independent systems, it is tied either to error indicators, energy laws, or complementarity residuals (Liao et al., 2020, Andreia et al., 2022, Steiner et al., 2014).

1. Depth as pseudo-time in residual and neural-ODE-like models

A direct pseudo-time interpretation arises in residual architectures when a residual block is written as

$$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$

which is exactly the forward Euler discretization of

$$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$

with $t_j=jh$, $T=Dh$, and $\mathbf{W}(t_j)\approx \mathbf{W}_j$ (Yang et al., 2019). In this interpretation, $t$ is not physical time but a pseudo-time along depth, and $\mathbf{y}(t)$ is the feature representation as it evolves through the network.

Adaptive pseudo-time stepping in this setting replaces the fixed step $\tau$0 by learned, layer-dependent increments $\tau$1:

$\tau$2

The pseudo-time grid and the effective evolution time are therefore learned jointly with the network weights (Yang et al., 2019). The motivation is explicitly compared with adaptive Runge–Kutta–Fehlberg (RKF), which updates the next step by

$\tau$3

with $\tau$4, but the network setting does not estimate local truncation error explicitly (Yang et al., 2019).

The central stability bound established for variable-step ResNets is

$\tau$5

where $\tau$6 is the spectral norm of the weight matrix in block $\tau$7 (Yang et al., 2019). This identifies the product “step size $\tau$8 spectral norm” as the relevant local amplification mechanism. The same work also shows that very small $\tau$9 can make layerwise gradients nearly identical,

$t$0

thereby smoothing the feature evolution too strongly and reducing representational power (Yang et al., 2019). This yields the characteristic trade-off: small pseudo-time steps improve stability and robustness, while larger steps improve discriminative transformation but can destabilize training.

The implemented controller is a data-independent, parameter-driven LSTM. At each residual block, the convolution kernels $t$1 are compressed, mapped through an input transform, processed recurrently, and converted into a channel-wise step-size vector

$t$2

with sigmoid output so that each channel’s step lies in $t$3 (Yang et al., 2019). The residual update becomes

$t$4

Because the controller is data-independent, it is discarded after training and only the learned step-size vectors are retained at inference, giving “no additional overhead in inference phase” (Yang et al., 2019). On ImageNet, the method improves top-1 accuracy across ResNet, ResNeXt, and SENet variants while leaving inference GFLOPs essentially unchanged; for example, on ResNet-50 the baseline has 25.56M parameters and 3.86 GFLOPs, while TSC$t$5-ResNet-50 has 27.83M parameters and 3.89 GFLOPs during training, but 25.57M parameters and 3.86 GFLOPs at inference (Yang et al., 2019). The learned step-size profiles are small in earlier layers and larger in late layers, suggesting a staged pseudo-time allocation rather than a uniform depth discretization (Yang et al., 2019).

2. Geometric reparametrization and structure-preserving pseudo-time

A second major line of work interprets adaptive stepping as a time reparametrization. For Hamiltonian systems, naive variable-step symplectic integration typically fails for two distinct reasons. If the step depends on time alone, $t$6, backward error analysis shows that the method can remain symplectic while developing parametric instabilities caused by resonance between oscillations of $t$7 and orbital motion (Richardson et al., 2011). If the step depends on phase-space position, $t$8, the transformation

$t$9

leads to equations in pseudo-time $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$0 that are generally no longer in canonical Hamiltonian form (Richardson et al., 2011). The proposed remedy is to integrate with constant step in $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$1 while preserving the transformed geometric structure either by an extended phase-space method or by a non-canonical mixed-variable generating function (Richardson et al., 2011).

In the extended phase-space construction, time is promoted to a coordinate $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$2 with conjugate momentum $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$3, and an extended Hamiltonian

$$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$4

is introduced, with $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$5 (Richardson et al., 2011). Constant pseudo-time steps $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$6 then induce adaptive physical steps $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$7. The alternative non-canonical method directly preserves the two-form $$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$8 in one degree of freedom by constructing a Poisson integrator via generating functions (Richardson et al., 2011). In both cases, pseudo-time stepping is not an auxiliary implementation trick but the mathematically correct mechanism by which adaptivity is reconciled with geometric structure.

The same reparametrization principle is used in accelerator simulation with the Boris–Buneman integrator (Toggweiler et al., 2012). There, adaptive pseudo-time stepping is realized through a Sundman transformation

$$\mathbf{y}_{j+1} = \mathbf{y}_j + h\,\mathcal{F}(\mathbf{y}_j,\mathbf{W}_j),\quad j=0,\ldots,D-1,$$9

with $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$0, so that a constant step $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$1 in pseudo-time gives a variable physical step $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$2 (Toggweiler et al., 2012). The choice of $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$3 is driven not by a local truncation error estimate but by a physically motivated state function proportional to a power of the maximum self-field acceleration: $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$4 with $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$5, hence $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$6 (Toggweiler et al., 2012). Strong space charge therefore produces small physical steps, while weak space charge produces large steps. This is embedded into an adaptive multiple-time-stepping Boris–Buneman method in which the outer step $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$7 controls the expensive self-field updates, while an inner step $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$8 remains close to a prescribed $$\dot{\mathbf{y}}(t)=\mathcal{F}(\mathbf{y}(t),\mathbf{W}(t)),\quad \mathbf{y}(0)=\mathbf{y}_0,\ t\in[0,T],$$9 for external fields (Toggweiler et al., 2012). In a photoinjector test, the adaptive method achieves a given emittance error with dramatically fewer self-field solves than a fixed-step MTS method; in a cyclotron test, where space charge varies only mildly, the gain disappears, indicating that the usefulness of adaptive pseudo-time depends strongly on scale separation and force variability (Toggweiler et al., 2012).

These geometric and force-based formulations suggest a broad interpretation: pseudo-time stepping is a controlled reparametrization of evolution. The step size is encoded in $t_j=jh$0, and the numerical problem becomes one of integrating the $t_j=jh$1-system with a structure-preserving method.

3. Stability-driven adaptivity in stiff deterministic and stochastic systems

A third strand uses adaptive pseudo-time steps to stabilize explicit or semi-explicit discretizations in stiff settings. In stiff ODEs, the basic difficulty is the classical explicit stability restriction. For $t_j=jh$2, explicit Euler is stable only if

$t_j=jh$3

The proposed explicit stabilization strategy is to take a large unstable step $t_j=jh$4, followed by $t_j=jh$5 small stable steps $t_j=jh$6, so that the composite amplification

$t_j=jh$7

remains stable (Eriksson et al., 2012). For $t_j=jh$8, the number of stabilizing small steps scales only logarithmically,

$t_j=jh$9

which implies that a small amount of local pseudo-time refinement can stabilize a much larger effective step (Eriksson et al., 2012). For systems with a spectral gap this yields large speedups; for parabolic spectra, Chebyshev and dyadic damping strategies extend the real-axis stability interval with cost reductions of order $T=Dh$0 and $T=Dh$1, respectively (Eriksson et al., 2012). This suggests a pseudo-time interpretation in which large outer steps advance the slow dynamics and short bursts of small steps quench stiff modes.

In SDEs with multiplicative, non-globally Lipschitz coefficients,

$T=Dh$2

the adaptive rule

$T=Dh$3

is designed so that

$T=Dh$4

thereby uniformly bounding the drift and diffusion responses per step (Kelly et al., 2017). The analysis is delicate because the Brownian increments over random steps are no longer independent, so the step times are constructed as stopping times and the proof uses a discrete Itô formula with conditional moments together with martingale convergence lemmas (Kelly et al., 2017). Under $T=Dh$5, the adaptive Euler–Maruyama scheme satisfies

$T=Dh$6

for sufficiently small $T=Dh$7, while if $T=Dh$8, then $T=Dh$9 (Kelly et al., 2017). The same rule also yields positivity with arbitrarily high probability on any finite horizon (Kelly et al., 2017). Here pseudo-time stepping is not guided by a classical error estimate but by a Lyapunov-driven control of qualitative dynamical properties.

The corresponding SPDE literature uses state-dependent, admissible time-step rules to control explicit exponential integrators in the presence of one-sided Lipschitz but non-globally Lipschitz drift (Campbell et al., 2018). For semilinear SPDEs

$\mathbf{W}(t_j)\approx \mathbf{W}_j$0

the adaptive explicit exponential schemes

$\mathbf{W}(t_j)\approx \mathbf{W}_j$1

and

$\mathbf{W}(t_j)\approx \mathbf{W}_j$2

use adaptive $\mathbf{W}(t_j)\approx \mathbf{W}_j$3 chosen so that along the discrete trajectory

$\mathbf{W}(t_j)\approx \mathbf{W}_j$4

holds (Campbell et al., 2018). Admissible step controls include

$\mathbf{W}(t_j)\approx \mathbf{W}_j$5

among others (Campbell et al., 2018). Under these rules, the full discretization achieves strong convergence

$\mathbf{W}(t_j)\approx \mathbf{W}_j$6

which is essentially the same temporal order as in the globally Lipschitz case (Campbell et al., 2018). Numerically, these adaptive schemes outperform fixed-step tamed or stopped alternatives at equal computational effort in Allen–Cahn and Swift–Hohenberg SPDEs (Campbell et al., 2018).

A closely related but distinct phenomenon appears in the Poisson–Nernst–Planck equations with generalized Frumkin–Butler–Volmer boundary conditions (Yan et al., 2017). The paper develops adaptive VSBDF2 and VSSBDF2 schemes using a coarse–fine local truncation error estimator. For the semi-implicit scheme, when the underlying dynamics approach a steady state, the adaptive time step does not grow indefinitely but stabilizes at a limiting size $\mathbf{W}(t_j)\approx \mathbf{W}_j$7, and the numerical solution only “nearly” converges to the steady state (Yan et al., 2017). The paper states that in the companion article the linearized SBDF2 scheme about the steady-state solution is shown to be conditionally stable, and that this conditional stability is the cause of the adaptive time-stepper’s behaviour (Yan et al., 2017). By contrast, the fully implicit BDF2 adaptive scheme is not subject to such time-step constraints, but the required nonlinear solve makes it significantly slower (Yan et al., 2017). This suggests that in pseudo-time marching to steady states, adaptivity alone does not remove the fundamental stability ceiling of a conditionally stable base integrator.

4. Energy-dissipative pseudo-time stepping in gradient flows and rate-independent systems

In dissipative PDEs, adaptive pseudo-time stepping is often governed by discrete energy laws rather than by classical stability regions. For the phase field crystal (PFC) equation,

$\mathbf{W}(t_j)\approx \mathbf{W}_j$8

the adaptive variable-step BDF2 discretization is

$\mathbf{W}(t_j)\approx \mathbf{W}_j$9

with variable-step BDF2 operator

$t$0

and step ratios $t$1 (Liao et al., 2020). The key structural condition is the zero-stability bound

$t$2

under which the BDF2 convolution kernels remain positive definite (Liao et al., 2020). The discrete modified energy

$t$3

is then nonincreasing under an additional mild time-step restriction (Liao et al., 2020). The adaptive controller uses the cheap relative-change indicator

$t$4

together with a safety factor $t$5 and the cap $t$6, to accept, reject, and update the next step (Liao et al., 2020). The method captures fast transients with small steps and slow coarsening with large steps while preserving the modified energy law and mass conservation (Liao et al., 2020). In a large-domain polycrystal growth simulation, the adaptive strategy uses about 393 steps versus 1000 uniform steps to reach $t$7 while producing essentially the same energy curve (Liao et al., 2020).

For rate-independent systems with non-convex energy, continuous physical time is often not the appropriate evolution parameter because the solution can jump. The paper on $t$8-parametrized balanced viscosity solutions introduces an explicit curve parameter $t$9 and interprets the evolution as a pseudo-time path

$\mathbf{y}(t)$0

satisfying

$\mathbf{y}(t)$1

with complementarity

$\mathbf{y}(t)$2

and an energy-dissipation balance (Andreia et al., 2022). The discrete pseudo-time steps $\mathbf{y}(t)$3 advance a local incremental stationary scheme

$\mathbf{y}(t)$4

and the physical time is updated by

$\mathbf{y}(t)$5

Hence $\mathbf{y}(t)$6 during viscous jumps and $\mathbf{y}(t)$7 in sticking or rate-independent slip (Andreia et al., 2022). Adaptivity is driven by two residuals: one for local stability/complementarity,

$\mathbf{y}(t)$8

and one for the energy-dissipation balance, $\mathbf{y}(t)$9, built from the discrete energy identity (Andreia et al., 2022). Steps are accepted if both residuals are below a tolerance, doubled if both are below half the tolerance, and halved otherwise (Andreia et al., 2022). The resulting scheme increases the pseudo-time step during sticking and in viscous jumps, and numerical experiments show a significant reduction in the number of steps relative to a uniform pseudo-time discretization (Andreia et al., 2022). This is one of the clearest cases where “pseudo-time stepping” is literal: the step size is attached to an auxiliary curve parameter rather than to physical time.

A recent filtered implicit framework also fits this energy-stability-oriented perspective. A variable-step pre-filter/post-filter method based on Implicit Euler constructs a second/third-order embedded pair,

$\tau$00

with

$\tau$01

and an explicit variable-step formula for $\tau$02 (McGovern, 1 Apr 2025). The difference between the second- and third-order outputs gives the estimator

$\tau$03

and a simple halving/doubling controller adjusts $\tau$04 (McGovern, 1 Apr 2025). The method is proposed in the standard ODE initial-value setting but is described as directly relevant to pseudo-time stepping toward steady states because it retains the Implicit Euler solve while allowing larger, adaptively selected steps once the solution changes slowly (McGovern, 1 Apr 2025).

5. Goal-oriented and locality-aware adaptivity in hyperbolic and particle systems

For weakly instationary hyperbolic conservation laws, pseudo-time adaptivity is often best understood as goal-oriented global time-step control. In scalar conservation laws and the 2D Euler equations, the error in a target functional is represented as

$\tau$05

where $\tau$06 is the residual and $\tau$07 solves the adjoint problem (Steiner et al., 2014, Steiner et al., 2014). Using space–time $\tau$08 projections, the error splits into a temporal part $\tau$09 and a spatial part $\tau$10: $\tau$11 up to boundary and linearization terms (Steiner et al., 2014, Steiner et al., 2014). The main innovation is that, because only $\tau$12 matters by Galerkin orthogonality, it is sufficient to compute the dual gradient $\tau$13, which satisfies a conservation law instead of a transport equation: $\tau$14 in 1D, or more generally a conservative dual PDE in multiple dimensions (Steiner et al., 2014, Noelle et al., 2014, Steiner et al., 2014). This allows the same finite-volume or DG machinery as in the primal solver.

The time indicator $\tau$15 is then used to design new step sizes by an equidistribution principle. In the scalar case, the update is written as

$\tau$16

so that stationary regions, where $\tau$17 is tiny, receive very large timesteps, while transient regions receive small ones (Steiner et al., 2014). In the Euler bump problem, the combined space–time adaptive method takes timesteps with CFL numbers of order $\tau$18 in stationary intervals and small timesteps when a perturbation enters the flow field (Noelle et al., 2014, Steiner et al., 2014). The method becomes even more efficient when coupled with a hybrid explicit–implicit strategy: use explicit steps when the equidistributed timestep would imply a small CFL and implicit steps when the timestep is large (Noelle et al., 2014). This is directly analogous to adaptive pseudo-time continuation, except that the controller is driven by an adjoint-weighted functional error rather than by a residual norm.

A more localized multirate interpretation appears in CFD–DEM particle integration. The method partitions particles into boxes using orthogonal recursive bisection and assigns each box a local timestep based on whether particles are close to collision partners (Sitaraman et al., 2018). For particle $\tau$19, the contact time is

$\tau$20

and if a potential collision partner lies within $\tau$21, the local timestep is set to

$\tau$22

Otherwise the particle uses a much larger sub-iteration step $\tau$23 (Sitaraman et al., 2018). Each box takes its own local timestep $\tau$24, which is the minimum over the particles in that box, and advances independently within each fluid time slab (Sitaraman et al., 2018). This is a clear example of local pseudo-time allocation in a multiphysics setting: the global physical time is fixed by the fluid step $\tau$25, while the DEM subsystem evolves in locally adaptive pseudo-time substeps whose density reflects local collisional stiffness.

A related space–time adaptive strategy appears in convection–diffusion finite differences (Malheiro et al., 2021). There, the scheme starts with the highest-order centered spatial discretization and RK$\tau$26 in time, using a global time step given by the largest CFL allowed by linear stability, then detects local instabilities a posteriori and selectively switches troubled nodes to an upwind spatial scheme and a more dissipative RK$\tau$27 time scheme while keeping the same global $\tau$28 (Malheiro et al., 2021). This suggests a different kind of adaptive pseudo-time stepping: instead of shrinking the global step, one modifies the local scheme so that the large step remains admissible.

6. Cross-cutting principles and recurrent trade-offs

Across these disparate formulations, several recurrent principles define adaptive pseudo-time stepping.

The first is that step-size control is tied to a structural quantity rather than to time alone. In residual networks, the decisive combination is $\tau$29 (Yang et al., 2019). In Boris–Buneman particle acceleration, it is the maximum self-field acceleration $\tau$30 through $\tau$31 (Toggweiler et al., 2012). In non-Lipschitz SDEs and SPDEs, it is the size of the drift and diffusion responses through bounds such as

$\tau$32

or

$\tau$33

(Kelly et al., 2017, Campbell et al., 2018). In rate-independent systems, it is the complementarity residual and the energy-dissipation residual (Andreia et al., 2022). In hyperbolic PDEs, it is the adjoint-weighted time contribution to a goal functional (Steiner et al., 2014, Noelle et al., 2014, Steiner et al., 2014).

The second principle is that pseudo-time is often more fundamental than physical time. In Hamiltonian and accelerator problems, the evolution is explicitly rewritten in terms of $\tau$34 with $\tau$35 (Richardson et al., 2011, Toggweiler et al., 2012). In rate-independent systems, $\tau$36 is the primary evolution variable and physical time is recovered as part of the solution (Andreia et al., 2022). In residual networks, depth itself plays the role of pseudo-time (Yang et al., 2019). Even when the papers work in physical time, as in PFC or weakly instationary Euler flows, the adaptive mechanism behaves as if it were allocating pseudo-time density to where the solution or the target functional changes most (Liao et al., 2020, Steiner et al., 2014).

The third principle is that adaptivity can preserve or destroy structure depending on how it is introduced. Variable steps inserted naively into symplectic schemes can cause resonant instabilities or destroy canonical Hamiltonian form (Richardson et al., 2011). Semi-implicit BDF2 can exhibit a limiting step and only “nearly” converge to steady state because of conditional stability (Yan et al., 2017). By contrast, adaptive variable-step BDF2 for the PFC model preserves a modified energy law under the sharp ratio bound $\tau$37 (Liao et al., 2020), and the rate-independent pseudo-time algorithm converges to parametrized balanced viscosity solutions as the tolerance tends to zero (Andreia et al., 2022).

The fourth principle is that training-time or offline adaptation can be separated from runtime cost. The ResNet controller is used only during training; after training, the learned pseudo-time schedule is frozen and incurs negligible inference overhead (Yang et al., 2019). In accelerator simulation, the adaptive rule is algebraic in the current state and cheaper than repeated local error estimates or rejected steps (Toggweiler et al., 2012). In adjoint-based Euler timestep control, the dual solve is performed on a coarse grid, and the expensive fine-grid run uses the resulting step sequence (Noelle et al., 2014, Steiner et al., 2014). This suggests a broad practical pattern: adaptive pseudo-time stepping is often most useful when the adaptation mechanism shapes the evolution law while leaving the dominant runtime kernel unchanged.

A plausible implication is that “adaptive pseudo-time stepping” is best regarded not as a single algorithm but as a design paradigm. Its core ingredients are: a surrogate evolution variable; a structurally meaningful step-selection signal; and a discretization or controller that allows nonuniform stepping without destroying the qualitative behavior one wishes to preserve. The literature surveyed here shows that this paradigm can be instantiated by learned channelwise depth increments (Yang et al., 2019), state-dependent geometric time rescaling (Richardson et al., 2011, Toggweiler et al., 2012), positivity- and stability-preserving explicit step control (Kelly et al., 2017, Campbell et al., 2018, Eriksson et al., 2012), energy-law-aware variable-step multistep schemes (Liao et al., 2020, McGovern, 1 Apr 2025), and adjoint-driven goal-oriented timestep allocation (Steiner et al., 2014, Noelle et al., 2014, Steiner et al., 2014).