Delay Ordinary Differential Equations

• Delay ODEs are differential equations where the derivative at time t depends on past states, requiring a history function for well-posedness.
• They appear with constant, time-dependent, or state-dependent delays, and are studied using tools like functional analysis, bifurcation theory, and numerical schemes.
• Advanced methods such as Lie symmetry, variational principles, and neural approaches enable insights into stability, linearization, and diffusion-like large-delay behaviors.

Delay ordinary differential equations (DODEs) are ordinary differential equations whose evolution at time $t$ depends on one or more past states rather than only on the instantaneous value. In the retarded form, a first-order DODE is commonly written as

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$

so the state at time $t$ is the history segment $x_t$ on a delay interval rather than a finite-dimensional point alone (Gselmann et al., 2016). In Lie-symmetry and variational treatments, one often works with a delay ordinary differential system (DODS), where the differential equation is paired with an explicit delay relation such as $x_-=g(\cdot)$ or $t_-=g(\cdot)$; this reflects the fact that, unlike an ODE, a delayed equation is not fully specified until the delayed argument is fixed (Dorodnitsyn et al., 2017). DODEs occur with constant, time-dependent, piecewise-constant, and state-dependent delays, and their analysis draws on functional analysis, operator theory, asymptotics, symmetry methods, bifurcation theory, numerical approximation, and, more recently, data-driven identification (Aigner et al., 2024).

1. Formulations, histories, and delay relations

A standard first-order retarded problem consists of

$x'(t)=F(t,x_t), \qquad t\ge t_0,$

together with an initial history

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$

where the delayed variables may be represented by continuous functions $g_1,\dots,g_m$ satisfying $\gamma\le g_j(t)\le t$ and

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$0

(Gselmann et al., 2016). This formulation makes explicit that the natural initial datum is a function on an interval, not a point value. In discrete-delay form, one often writes

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$1

with constant delays $x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$2 (Pecile et al., 2024).

In symmetry-based treatments, the equation and the delay relation are kept on equal footing. A first-order DODS may be written as

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$3

while a second-order DODS takes the form

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$4

(Dorodnitsyn et al., 2017). This formulation rules out a common simplification: the delayed argument is not always externally prescribed, and in state-dependent problems it is itself part of the dynamics.

The distinction from classical ODEs is structural, not merely notational. Delay equations evolve on spaces of histories, and even scalar equations are infinite-dimensional in this sense. Large-delay analysis makes this especially vivid: for

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$5

the delay interval can act as a pseudo-spatial domain, so that the dynamics acquire a continuum-limit description more akin to a PDE than to a finite-dimensional ODE (Kozyreff, 2023). A related misconception is that continuity of the prehistory should suffice for a well-posed initial-value problem. For ordinary state-dependent delay equations, this is not generally enough for uniqueness; an $x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$6-based framework with controlled prehistories is used precisely because low regularity can cause nonuniqueness (Aigner et al., 2024).

2. Well-posedness, continuous dependence, and stability

For first-order retarded systems, a central well-posedness hypothesis is a Lipschitz bound in the delayed variables,

$x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$7

with continuous $x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$8. Under this assumption, the initial-value problem has at most one solution on any subinterval, and solutions depend continuously on the initial function $x'(t)=F(t,x_t), \qquad x_t(\theta)=x(t+\theta),$9 (Gselmann et al., 2016). If $t$0 and $t$1 solve the same equation with different initial data, then

$t$2

which is the delay analogue of classical continuous dependence (Gselmann et al., 2016).

The same framework yields Hyers-Ulam stability. If a continuous function satisfies the DODE only approximately,

$t$3

and shares the same initial history, then there exists a unique exact solution $t$4 such that

$t$5

This quantifies the principle that a uniformly small residual in the delayed differential law forces proximity to a genuine solution (Gselmann et al., 2016).

For linear retarded functional differential equations, the variation-of-constants formula is subtler than for ODEs because the natural fundamental object is a principal fundamental matrix built from a discontinuous initial history. A mild-solution framework on $t$6 resolves this by allowing histories that are $t$7 on $t$8 but possess a value at $t$9. The resulting principal fundamental matrix $x_t$0 satisfies a delay analogue of the matrix-exponential role in ODE theory, and the variation-of-constants representation then supports linearized stability and a Poincaré-Lyapunov theorem for DDEs without assuming uniqueness of nonlinear solutions (Nishiguchi, 2022).

State-dependent delay equations require a different analytic mechanism. In a Hilbert-space setting, one studies

$x_t$1

in exponentially weighted Sobolev spaces $x_t$2, where the history operator becomes Lipschitz with norm at most $x_t$3. After splitting off the prehistory and projecting onto the closed convex set

$x_t$4

the problem is reduced to a contraction mapping. This yields local existence and uniqueness, continuous dependence on initial prehistories and on the right-hand side, and a continuation criterion in which either the maximal existence time is the whole interval or

$x_t$5

(Aigner et al., 2024).

A broader operator-theoretic perspective treats delay and memory equations on weighted Hilbert spaces over the whole real line via an invertible time-derivative operator $x_t$6. Within this setting, delays, neutral terms, and convolutions are handled by functional calculus, and the resulting solution operators are causal for sufficiently large $x_t$7 (Kalauch et al., 2012). This suggests that causality, rather than merely existence, is a fundamental organizing principle in abstract DODE theory.

3. Symmetry classification, linearization, and variational structure

Lie symmetry methods reveal a sharp divide between genuinely nonlinear DODEs and linear or linearizable ones. For first-order DODSs,

$x_t$8

the admitted Lie symmetry algebra can have dimension $x_t$9. Infinite-dimensional algebras occur for linear DODSs with solution-independent delay relations and for nonlinear systems that are linearizable by invertible point transformations. A decisive criterion is the presence of two linearly connected symmetries: if a DODS admits such a two-dimensional subalgebra, it can be transformed into a linear DODE with solution-independent delay (Dorodnitsyn et al., 2017). For linear first-order DODSs, the superposition principle produces an infinite-dimensional symmetry algebra generated by vector fields of the form

$x_-=g(\cdot)$0

with $x_-=g(\cdot)$1 a solution of the homogeneous equation (Dorodnitsyn et al., 2017).

For second-order DODSs with one delay point,

$x_-=g(\cdot)$2

the symmetry threshold is higher. Genuinely nonlinear systems have symmetry algebras of dimension $x_-=g(\cdot)$3, whereas $x_-=g(\cdot)$4 characterizes linear or linearizable systems. A four-dimensional symmetry algebra realized by linearly connected vector fields implies linearizability, and linear homogeneous second-order DODSs again admit infinite-dimensional superposition symmetries $x_-=g(\cdot)$5 (Dorodnitsyn et al., 2019).

Variational DODE theory departs even more strongly from the ODE pattern. A first-order delay functional with one delay,

$x_-=g(\cdot)$6

generates second-order DODEs with two delays under variation (Dorodnitsyn et al., 2023). Two distinct variational equations arise: the vertical variation yields the Elsgolts equation,

$x_-=g(\cdot)$7

while horizontal variation produces a separate equation in $x_-=g(\cdot)$8. Unlike the non-delayed case, these equations are not generally equivalent (Dorodnitsyn et al., 2023).

The corresponding Noether theory is a mixed differential-difference theory. For the generator $x_-=g(\cdot)$9, the delay Noether operator identity has the form

$t_-=g(\cdot)$0

linking Lagrangian invariance, the appropriate variational equations, and conserved quantities (Dorodnitsyn et al., 2023). Consequently, first integrals split into differential first integrals $t_-=g(\cdot)$1 and difference first integrals $t_-=g(\cdot)$2, a distinction without an ODE analogue (Dorodnitsyn et al., 2023).

The Hamiltonian side can also be reconstructed. A delay analogue of the Legendre transform relates a delay Lagrangian to a Hamiltonian $t_-=g(\cdot)$3, subject to a compatibility condition $t_-=g(\cdot)$4, and produces first-order Hamiltonian DODEs with two delays (Dorodnitsyn et al., 2024). The Hamiltonian Noether identity then yields a conservation law of the same mixed type,

$t_-=g(\cdot)$5

for invariant Hamiltonian functionals (Dorodnitsyn et al., 2024). This suggests that classical Lagrangian-Hamiltonian duality survives in the delayed setting only after the delay structure itself is built into the transform.

4. Spectral structure, large-delay asymptotics, and PDE embeddings

Large delay changes the effective geometry of a DODE. For

$t_-=g(\cdot)$6

multiple-scale analysis introduces a strained fast variable

$t_-=g(\cdot)$7

and a slow variable

$t_-=g(\cdot)$8

with $t_-=g(\cdot)$9. The leading-order problem is a periodic map in $x'(t)=F(t,x_t), \qquad t\ge t_0,$0, not a differential equation, and the solvability condition at the next nontrivial order yields

$x'(t)=F(t,x_t), \qquad t\ge t_0,$1

Thus the large-delay limit of this linear DODE is governed by a diffusion equation on a unit-periodic pseudo-space (Kozyreff, 2023). The associated spectrum becomes dense as $x'(t)=F(t,x_t), \qquad t\ge t_0,$2, and the real parts scale like $x'(t)=F(t,x_t), \qquad t\ge t_0,$3, matching a diffusive dispersion relation (Kozyreff, 2023).

A different route from DDEs to spatial dynamics starts from the history transformation

$x'(t)=F(t,x_t), \qquad t\ge t_0,$4

For

$x'(t)=F(t,x_t), \qquad t\ge t_0,$5

this yields the advection equation

$x'(t)=F(t,x_t), \qquad t\ge t_0,$6

with a nonlinear boundary condition at $x'(t)=F(t,x_t), \qquad t\ge t_0,$7. A Galerkin expansion in shifted Legendre polynomials then converts the infinite-dimensional problem into a finite-dimensional ODE system. For smooth solutions, the approximation error decreases exponentially as the number of retained modes increases, and the reduced ODEs can approximate strange attractors and Lyapunov exponents of chaotic DDEs with fewer modes than a standard method of lines (Sadath et al., 2018).

Pseudospectral collocation yields a closely related reduction tailored to bifurcation analysis. Approximating the history by polynomial interpolation at nodes $x'(t)=F(t,x_t), \qquad t\ge t_0,$8 produces an ODE whose universal part encodes translation along the history segment and whose extension equation imposes the original DDE at $x'(t)=F(t,x_t), \qquad t\ge t_0,$9. The discrete characteristic function $x(t)=\phi(t), \qquad t\in[\gamma,t_0],$0 converges rapidly to the DDE characteristic function $x(t)=\phi(t), \qquad t\in[\gamma,t_0],$1, and Hopf points, frequencies, and direction coefficients of the ODE approximation converge to those of the original DDE (Wolff et al., 2020). A frequent misunderstanding is that such reductions are merely pragmatic discretizations; the convergence results show that they are also structure-preserving approximations of the generator of translation along solution histories (Wolff et al., 2020).

5. Numerical continuation, bifurcation computation, and solver architectures

The infinite-dimensional phase space of DODEs complicates direct bifurcation computation, but it can sometimes be reduced exactly at critical spectral configurations. For a one-delay, two-parameter equation

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$2

Takens-Bogdanov points are defined in the abstract phase space $x(t)=\phi(t), \qquad t\in[\gamma,t_0],$3 by a double zero eigenvalue with a Jordan chain. By exploiting explicit bases of the eigenspace and adjoint eigenspace, the Banach-space defining system can be reduced to a finite-dimensional algebraic system

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$4

At a quadratic Takens-Bogdanov point, this reduced system is regular, so standard Newton iteration applies directly (Xu et al., 2012).

For Hopf bifurcation, pseudospectral reduction makes standard ODE continuation software available. The discretized ODE can be continued in tools such as MatCont, while DDE-specific packages such as DDE-BIFTOOL remain available for the original problem. The key caveat is that the first Lyapunov coefficient reported by generic ODE software is not directly comparable to the DDE direction coefficient unless the normalization conventions of eigenvectors are matched (Wolff et al., 2020).

General-purpose time integration has also shifted toward reuse of ODE infrastructure. DelayDiffEq constructs DDE solvers by recursively embedding an ODE solver within itself: a continuous ODE interpolant supplies the history, and the outer solver advances the delayed problem while sharing caches so that interpolation and stage values remain consistent (Widmann et al., 2022). This design supports the classical method of steps,

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$5

but also permits unconstrained stepping via fixed-point iterations and handles discontinuity propagation explicitly (Widmann et al., 2022). Because it inherits the solver ecosystem of OrdinaryDiffEq, the framework can use explicit Runge-Kutta methods for non-stiff problems and stiff integrators such as KenCarp5, Rosenbrock methods, SDIRK methods, BDF methods, and FIRK methods for stiff delayed systems (Widmann et al., 2022). This suggests that, numerically, the boundary between ODE and DODE solvers is increasingly architectural rather than categorical.

6. Data-driven discovery and neural delay models

Recent data-driven work treats DODE identification as a sparse-regression or neural-system-identification problem on delayed feature spaces. In the SINDy extension to DDEs with discrete delays, one augments the library by delayed samples,

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$6

and solves

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$7

by sparse regression, typically with STLSQ or LASSO (Pecile et al., 2024). Unknown delays are then identified by minimizing the reconstruction residual

$x(t)=\phi(t), \qquad t\in[\gamma,t_0],$8

with Bayesian optimization used as an outer loop to reduce the number of SINDy calls relative to brute-force search (Pecile et al., 2024). Across the reported test cases, the average reduction in error-function evaluations or SINDy calls is about $x(t)=\phi(t), \qquad t\in[\gamma,t_0],$9 (Pecile et al., 2024).

The methodology extends to multiple delays, systems of DDEs, and unknown non-multiplicative parameters. Reported examples include the delayed logistic equation, a delayed SIR model, the Mackey-Glass equation, and multiple-delay polynomial examples (Pecile et al., 2024). A practical limitation is that the approach is designed for autonomous DDEs with constant discrete delays and is not yet designed for distributed delays (Pecile et al., 2024).

Neural delay models make the delay explicit in the learned vector field. Neural State-Dependent DDEs parameterize

$g_1,\dots,g_m$0

with a history function $g_1,\dots,g_m$1, thereby covering multiple delays, constant delays, piecewise constant delays, time-dependent delays, and state-dependent delays, though not continuous or integral delays (Monsel et al., 2023). In reported benchmarks, Neural SDDDE attains the lowest test MSE on the time-dependent delayed logistic problem and on a state-dependent Mackey-Glass-type problem, while Neural Laplace is slightly better on the delayed diffusion PDE benchmark (Monsel et al., 2023).

A distinct neural branch replaces continuous delays by piecewise-constant ones. Neural Piecewise-Constant Delay Differential Equations use terms such as

$g_1,\dots,g_m$2

and can incorporate multiple previous frozen states without augmenting the state dimension (Zhu et al., 2022). The framework includes shared-parameter and unshared-parameter versions and is proved to have universal approximation capability in the forms stated in the paper (Zhu et al., 2022). Reported image-classification results show consistent improvements over NODE and NDDE baselines on MNIST, CIFAR10, and SVHN, with the best values typically obtained by UNPCDDE5 in the reported tables (Zhu et al., 2022).

A broader synthesis distinguishes direct sparse identification of the DDE itself from identification of an ODE approximation obtained by pseudospectral collocation. In this terminology, E-SINDy preserves interpretability of the delayed law, whereas P-SINDy is more pragmatic when delays are numerous or only the maximal delay is known (Breda et al., 4 Dec 2025). The same synthesis also notes an important caution for neural delay models: good trajectory fitting does not by itself guarantee correct recovery of the underlying delay structure, especially when the data are not sufficiently rich (Breda et al., 4 Dec 2025).

DODEs therefore occupy a technically distinctive position among dynamical systems: they are ordinary differential equations in derivative order, but functional-differential equations in state space; they admit ODE-like notions of uniqueness, stability, bifurcation, symmetry, and conservation law, but each of these notions is modified by the presence of history, delay relations, and delayed variational structure. The modern theory shows that these modifications are not peripheral. They determine the correct state space, the right form of well-posedness, the relevant symmetry algebra, the nature of conserved quantities, the architecture of numerical solvers, and even the design of data-driven models.