Pseudo-Time Stepping in Computational Methods

• Pseudo-time stepping is a strategy that reformulates difficult nonlinear updates as an artificial time evolution, driving residuals or state corrections toward convergence.
• It finds applications in steady-state CFD, dual-time stepping for unsteady flow, optimization, and neural PDE solvers by decoupling physical limitations from convergence efficiency.
• The method leverages techniques such as adaptive step sizing, modified Jacobian approximations, and underrelaxation to improve stability and convergence without eliminating inherent CFL restrictions.

Pseudo-time stepping denotes a family of formulations in which a difficult update, nonlinear solve, or large-step evolution is recast as an evolution in an artificial time variable, or as an inner sequence of auxiliary substeps, until a target residual, state, or macro-step condition is satisfied. In computational fluid dynamics this often means marching a steady residual to zero or enforcing an implicit physical-time discretization through inner pseudo-time iterations; in optimization and inverse problems it can mean following a solution manifold in an auxiliary parameter; in meshfree and neural methods it can mean wrapping explicit microdynamics or training dynamics inside a pseudo-implicit outer map (Li et al., 2016, Moujaes et al., 2024, Cao et al., 2023, Hart et al., 28 Aug 2025).

1. Canonical formulations and scope

Across the literature, pseudo-time stepping appears in several mathematically distinct but structurally related forms. The common pattern is that a target equation is not solved directly; instead, one constructs an auxiliary evolution whose stationary or terminal state realizes the desired update.

Context	Representative formulation	Representative paper
Steady CFD residual relaxation	$\frac{\partial x}{\partial \tau}+r(x)=0$	(Wang et al., 8 Apr 2025)
Dual-time stepping for unsteady flow	$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$	(Li et al., 2016)
Pseudo-transient continuation	$$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$	(Zandbergen et al., 2023)
Pseudo-implicit wrapper for explicit MPM	$\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$	(Jiang, 15 Aug 2025)
Continuation of inverse-problem minimizers	$$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$	(Hart et al., 28 Aug 2025)
PINN pseudo-time relaxation	$\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$	(Wang et al., 26 Apr 2026)

For steady-state CFD, pseudo-time is explicitly artificial: the objective is not faithful time integration, but convergence to a discrete steady solution. In one standard formulation, the nonlinear algebraic problem $\mathcal N(x)=0$ is replaced by $\partial x/\partial\tau + r(x)=0$, and pseudo-time iterations are judged by residual decay rather than by physical-time accuracy (Wang et al., 8 Apr 2025). In finite-element incompressible Navier–Stokes, pseudo-transient continuation plays the role of a globalization device for Newton’s method by shifting the Jacobian with a mass-like term, so that small pseudo-time steps produce damped corrections and large pseudo-time steps recover Newton-like behavior (Zandbergen et al., 2023).

For unsteady problems, pseudo-time stepping usually appears in the dual-time form. There the physical-time discretization is retained, but each physical step is solved through an inner pseudo-steady problem. In the gas-kinetic scheme, the BDF2-discretized update is rewritten as a steady problem in $t^*$, and convergence of the inner loop enforces the desired second-order physical-time advance (Li et al., 2016). This dual role—artificial inner time, physical outer time—is one of the defining features of modern pseudo-time methods.

The term is broader than residual relaxation alone. In Material Point Method substepping, many explicit CFL-sized substeps are run internally and then compressed into a macro-step grid field on the original configuration, so that the outer solver sees a large-step “pseudo-implicit” map even though the internal dynamics remain explicit (Jiang, 15 Aug 2025). In parameterized inverse problems, pseudo-time continuation is not a stabilizing relaxation of one residual equation; it is the ODE governing how a local minimizer moves along an auxiliary-parameter path (Hart et al., 28 Aug 2025). In neural PDE solvers, pseudo-time may be introduced in training space rather than state space, replacing direct residual minimization by a sequence of pseudo-time subproblems (Cao et al., 2023).

2. Dual-time stepping for unsteady simulation

A canonical high-fidelity use of pseudo-time stepping is the enforcement of implicit physical-time discretizations for unsteady flow. In the gas-kinetic scheme for compressible Navier–Stokes-level flow variables, the semi-discrete update is written in finite-volume form as

$$\bm{w}_I^{n+1}=\bm{w}_I^{n} - \frac{1}{\Omega_I} \int_{t_n}^{t_{n+1}} \sum_{J=1}^{N_{IF}}\bm{F}(t)_J S_J\,dt,$$

with conservative variables

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$0

and flux

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$1

The physical-time discretization is taken as BDF2: $$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$2 This is then converted into the pseudo-time evolution

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$3

with

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$4

When the pseudo-time residual vanishes, the inner artificial-time problem has enforced the outer BDF2 update (Li et al., 2016).

The inner solve is implicit in pseudo-time. Linearization gives

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$5

A notable feature is the Jacobian construction. Since the exact BGK-based flux Jacobian is difficult to compute, the residual Jacobian is approximated by conventional Navier–Stokes/Euler Jacobians and split into convective and viscous parts,

$$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$6

The convective part uses Roe linearization, while the viscous part is explicitly included as an addition absent from earlier implicit GKS treatments. The resulting large sparse nonsymmetric systems are solved by GMRES (Li et al., 2016).

The computational rationale is explicit. Direct explicit GKS for unsteady flow is accurate, but the step size is tied to CFL restrictions, and GKS flux evaluations are relatively expensive. Dual-time stepping decouples physical-step selection from explicit stability. In the reported laminar cylinder case at $$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$7, the explicit method used $$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$8, whereas dual time stepping used $$\frac{\partial (\Omega_I^{n+1}\bm{w}_I^{*})}{\partial t^*}=-\bm{R}_I^*(\bm{w}_I^*)$$9, $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$0 inner iterations, and pseudo-steady residual $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$1; to reach $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$2, the explicit method required $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$3 steps and the dual-time method $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$4 steps. For turbulent flow around a square cylinder, the reported cost was roughly one-tenth of the explicit scheme for the same simulated time; for transonic buffet on NACA0012, the explicit step was $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$5 while the dual-time method used $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$6 with $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$7 inner iterations and pseudo-steady residual $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$8 (Li et al., 2016).

A recurrent misconception is that dual time stepping removes all timestep restrictions. It does not. The paper is explicit that pseudo-time iterations remove the explicit stability restriction from the physical advance, but the physical step $$\mathbf{v}^{n+1}=\mathbf{v}^n-(\frac{1}{\Delta t^n}M+F'(\mathbf{v}^n))^{-1}F(\mathbf{v}^n)$$9 must still be chosen according to temporal accuracy requirements. In the NACA0012 buffet case, simulations with physical time steps from $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$0 to $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$1 showed convergence of the buffet response as the physical step decreased, and $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$2 was selected (Li et al., 2016).

3. Pseudo-time marching to steady solutions

In steady-state solvers, pseudo-time stepping is used as a nonlinear iteration rather than as an unsteady integrator. For continuous finite-element Euler equations with monolithic convex limiting, the fully implicit pseudo-time step is

$\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$3

and the steady target is

$\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$4

Here backward Euler in artificial time is combined with invariant-domain-preserving low-order operators and explicitly treated limited antidiffusive terms. The practical nonlinear iteration uses

$\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$5

or equivalently

$\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$6

The stopping criterion inside each pseudo-time step is unusual but deliberate: instead of a residual tolerance, admissibility of the candidate state is used, i.e. the iteration stops when the updated state lies in the invariant set $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$7 (Moujaes et al., 2024).

That formulation is paired with adaptive explicit underrelaxation: $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$8 where $\Delta t\, v_p^\star=x_p^{n+1,S}-x_p^n$9 is chosen by approximate minimization of nodal entropy residuals over a small candidate set. The reported conclusion is strong: some kind of explicit underrelaxation is an essential prerequisite for convergence in most tests; without it, convergence was obtained only in two simple scenarios. With adaptive relaxation and large pseudo-CFL values, the method converged robustly on test cases including the GAMM channel, transonic nozzle flows, NACA 0012, and a Mach 20 bow shock, with typical pseudo-time counts such as about $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$0 steps for the GAMM channel and about $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$1 for the subsonic nozzle case (Moujaes et al., 2024).

Pseudo-time stepping also serves as a steady-state convergence accelerator. In the MMRES method, the baseline steady CFD solver is written abstractly as

$$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$2

where $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$3 represents the implicit pseudo-time marching scheme, approximate inner solve, and any local time stepping or residual smoothing. MMRES periodically collects recent pseudo-time iterates, forms a mean-based affine trial space

$$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$4

approximates the residual image by

$$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$5

and computes a correction by reduced residual minimization: $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$6 The purpose is not to reduce inner linear-solver iterations, but to reduce the number of pseudo-time steps. Reported results include $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$7 to $$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$8 times wall-clock speedup across two- and three-dimensional steady-flow problems, together with stabilization of flows that otherwise remain effectively unsteady, including cylinder wakes with vortex shedding and transonic buffet without symmetry boundary conditions (Wang et al., 8 Apr 2025).

For quasilinear elliptic PDEs with solution-dependent diffusion, pseudo-time becomes part of a broader regularization framework. The core linearized solve is

$$\dot m(t)=-H(m(t),\theta(t))^{-1}B(m(t),\theta(t))\Delta\theta$$9

with update $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$0. Here $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$1 is a pseudo-time/Tikhonov term, $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$2 gives numerical dissipation with predicted residual rate $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$3, $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$4 is a Picard-like regularization term chosen to cancel estimated linearization error, and $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$5 scales the source term while the mesh is too coarse to resolve steep layers. This framework is expressly designed for pre-asymptotic adaptive finite-element computation, where standard Newton solves are unreliable (Pollock, 2016).

4. Step control, linear algebra, and variable-step stability

Once pseudo-time stepping is viewed as a nonlinear solver, the dominant implementation questions become step selection, Jacobian approximation, and inner linear algebra. The literature shows no single answer.

At one end are explicitly assembled Jacobian solves, as in dual-time GKS, where Roe convective terms, viscous couplings, and GMRES are combined in a standard sparse nonsymmetric implicit solve (Li et al., 2016). At another are deliberately low-order Jacobian approximations, as in the MCL Euler solver, where the practical linearization uses a robust low-order Jacobian and keeps the limited antidiffusive correction explicit for robustness (Moujaes et al., 2024). MMRES occupies a different level: it does not alter the inner linear solve, but periodically wraps the outer pseudo-time sequence with a reduced least-squares correction built from recent iterates (Wang et al., 8 Apr 2025).

Pseudo-transient continuation for steady incompressible Navier–Stokes makes the dependence on the pseudo-time step especially explicit: $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$6 In the local version, one uses an elementwise field $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$7 and replaces the scalar shift by $\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$8. Classical CFD software chooses

$\frac{u^k-u^{k-1}}{\tau}+R_{\mathrm{int}}[u^k]=0$9

with a global controller for $\mathcal N(x)=0$0, but the paper argues that convergence depends critically on this choice and proposes predicting elementwise local pseudo-time steps with a neural network using local geometry, solution values, and residual features (Zandbergen et al., 2023). Reported averages include, for example, $\mathcal N(x)=0$1 nonlinear iterations for the NN method versus $\mathcal N(x)=0$2 for $\mathcal N(x)=0$3 and $\mathcal N(x)=0$4 for $\mathcal N(x)=0$5 on one backward-facing step case, and $\mathcal N(x)=0$6 versus $\mathcal N(x)=0$7 and $\mathcal N(x)=0$8 on another (Zandbergen et al., 2023).

Variable-step stability is a separate issue. For stiff dissipative systems, the refactorized Dahlquist–Liniger–Nevanlinna scheme offers second-order accuracy and unconditional $\mathcal N(x)=0$9-stability for variable time steps. Its significance for pseudo-time or pseudo-transient continuation is practical: the method is exactly equivalent to one backward Euler solve bracketed by arithmetic pre- and postprocessing steps, so a legacy backward Euler code can be upgraded without changing its nonlinear or linear solver kernel (Layton et al., 2021). This suggests that variable-step pseudo-time continuation need not be confined to first-order backward Euler, provided the underlying dynamics satisfy the contractivity assumptions used in the DLN analysis.

5. Broader generalizations

Outside classical CFD, pseudo-time stepping has been generalized in several ways that enlarge the concept itself.

In Material Point Method substepping, the objective is not to obtain unconditional stability, but to expose a large-step interface for coupling, constraints, or multiphysics while retaining fully explicit internal mechanics. A macro-step $\partial x/\partial\tau + r(x)=0$0 is decomposed into $\partial x/\partial\tau + r(x)=0$1 explicit substeps of size $\partial x/\partial\tau + r(x)=0$2. After substepping, secant particle data are formed,

$\partial x/\partial\tau + r(x)=0$3

and a macro-step grid velocity is reconstructed on the original configuration. With mass lumping and an APIC-style approximation, the reconstruction becomes

$\partial x/\partial\tau + r(x)=0$4

The method is explicitly described as a pseudo-implicit wrapper around explicit MPM, and equally explicitly stated not to remove CFL restrictions (Jiang, 15 Aug 2025).

For parameterized PDE-constrained inverse problems, pseudo-time continuation follows a path of local minimizers rather than relaxing one residual to zero. If $\partial x/\partial\tau + r(x)=0$5 and $\partial x/\partial\tau + r(x)=0$6 is the corresponding optimizer, then

$\partial x/\partial\tau + r(x)=0$7

A predictor-corrector discretization advances this ODE, with Hessian systems solved by PCG and accelerated by an adaptive quasi-Newton preconditioner. The paper’s numerical examples emphasize that continuation can be substantially cheaper and more robust than re-solving each inverse problem independently from the nominal minimizer (Hart et al., 28 Aug 2025).

Neural PDE solvers provide a different generalization. In TSONN, pseudo-time is used to convert direct PDE-residual minimization into a sequence of subproblems. For steady problems the artificial evolution is

$\partial x/\partial\tau + r(x)=0$8

and the neural network is trained either against the explicit label

$\partial x/\partial\tau + r(x)=0$9

or the implicit relation

$t^*$0

The corresponding losses replace the standard PINN residual loss by one-step pseudo-time consistency terms, with reported gains on Laplace, Burgers, lid-driven cavity, and Allen–Cahn problems (Cao et al., 2023).

A later PINN study argues for a different interpretation. There, pseudo-time stepping is not presented primarily as a conditioning trick, but as a device that makes spurious solutions less viable when collocation points are resampled. The pseudo-time-relaxed loss is based on

$t^*$1

and the adaptive step is chosen from a finite-difference surrogate of the local residual Jacobian: $t^*$2 Across ten benchmarks, the adaptive method improved relative $t^*$3 error over both a strong baseline PINN pipeline and fixed-step pseudo-time stepping (Wang et al., 26 Apr 2026). This contrast with the conditioning-based explanation in TSONN suggests that “why pseudo-time works” remains context-dependent (Cao et al., 2023, Wang et al., 26 Apr 2026).

6. Misconceptions, limitations, and recurring design trade-offs

A persistent misconception is that pseudo-time stepping is synonymous with unconditional stability or with true implicitness. Several papers explicitly reject that reading. Dual-time stepping allows the physical step to be chosen by accuracy rather than explicit CFL limits, but it does not remove the need to resolve unsteady physics (Li et al., 2016). The MPM pseudo-implicit wrapper preserves explicit substep costs and explicit CFL restrictions; its “implicitness” is architectural, not a statement about altered stability theory (Jiang, 15 Aug 2025).

A second misconception is that pseudo-time stepping is a single method. The literature instead shows a taxonomy: residual relaxation for steady problems, dual-time enforcement of outer physical-time schemes, globalization of Newton iterations, continuation along parameter manifolds, substepping wrappers, and training-space relaxations in PINNs. This suggests that pseudo-time stepping is better understood as a design pattern than as a fixed algorithm.

Step-size choice is a recurring difficulty. In CFD globalization, local pseudo-time steps strongly affect convergence and motivate learned controllers (Zandbergen et al., 2023). In PINNs, the effectiveness of pseudo-time stepping depends critically on $t^*$4, and training loss alone is not a reliable tuning signal; this motivates Jacobian-surrogate adaptivity (Wang et al., 26 Apr 2026). In variable-step stiff integration, the DLN results indicate that preserving nonlinear dissipativity under changing steps is itself a nontrivial requirement, not an automatic consequence of using an implicit scheme (Layton et al., 2021).

Finally, pseudo-time stepping often succeeds by combining multiple mechanisms rather than one. In steady Euler solves it works together with invariant-domain preservation and entropy-based underrelaxation (Moujaes et al., 2024). In nonlinear diffusion it is one part of a larger regularization apparatus involving Tikhonov-like terms, numerical dissipation, Picard-like stabilization, and source scaling (Pollock, 2016). In inverse problems it is effective because the continuation path stays near a branch of local minimizers with positive definite Hessian (Hart et al., 28 Aug 2025). In PINNs, one line of work attributes the benefit to better-conditioned subproblems, while another attributes it to exposing spurious solutions under resampling (Cao et al., 2023, Wang et al., 26 Apr 2026).

Taken together, these formulations show that pseudo-time stepping is less a single discretization than a general strategy for replacing a difficult global update by a controlled auxiliary evolution. Its central technical questions are always the same: what is being relaxed, what artificial dynamics are introduced, what linear or nonlinear operators define the inner step, how step size is chosen, and which invariant, residual, or consistency property is meant to survive the auxiliary evolution.