Delayed Matrix Differential Equations

• Delayed matrix differential equations are functional differential equations whose state evolution depends on past values, often with noncommutative matrix coefficients and multiple delays.
• They model dynamic systems with internal memory in fields like control theory, signal processing, and biology, employing recursive and explicit solution frameworks.
• Recent advancements include explicit fundamental matrix solutions, fractional and distributed delay formulations, and high-accuracy numerical methods for robust stability analysis.

Delayed matrix differential equations (DMDEs) are functional differential equations whose unknown is a matrix- or vector-valued function, in which the derivative at each time depends explicitly on delayed (history) values of the function—possibly coupled through noncommuting matrix coefficients, multiple delays, distributed delays, fractional derivatives, or additional structure. DMDEs are pervasive in modeling dynamical systems with internal memory or communication lags, and constitute a crucial class for control, signal processing, mathematical biology, and engineering. This article details key structural paradigms, analytical and computational methodologies, and their extensions to fractional, distributed, and noncommutative frameworks, drawing directly on recent advances and rigorous research results.

1. Fundamental Forms and Noncommutativity in DMDEs

A prototypical first-order DMDE is

$$\dot{X}(t) = A_0 X(t-\sigma) + X(t-\sigma)A_1 + G(t), \qquad t \geq 0,$$

with prescribed initial matrix-valued history $X(t)=\Psi(t)$ for $t \in [-\sigma,0]$, where $A_0,\,A_1 \in \mathbb{R}^{d \times d}$ and $A_0 A_1 \neq A_1 A_0$ (the noncommutative case) (Asadzade et al., 21 Sep 2025). The discrete-time analogue is

$\Delta X(u) = A_0 X(u-m) + X(u-m) A_1 + G(u).$

The noncommutativity of system matrices fundamentally alters analytical solution strategies: classic approaches (e.g., via matrix exponentials or binomial identities) no longer apply, necessitating recursive auxiliary matrices and specialized fundamental solutions.

For $A_0, A_1$ commutative, binomial expansions yield closed formulas for the $Z$-function fundamental solution (e.g., $$Z(u) = \sum_{j=0}^{n} \binom{u-(j-1)m}{j} \sum_{\ell=0}^j \binom{j}{\ell} A_0^{j-\ell}A_1^{\ell}$$), but when $A_0, A_1$ do not commute, the solution construction follows an inductive table where each $Q_{u+1}(\ell \delta)$ is built from $A_0 Q_u((\ell-1)\delta) + Q_u((\ell-1)\delta)A_1$.

Explicit solutions for both homogeneous and inhomogeneous systems directly incorporate the structure of the initial history and the noncommutative convolution over delay intervals (Asadzade et al., 21 Sep 2025). This extends foundational works on delayed exponential and binomial expansions (Mahmudov, 2018, Mahmudov, 27 Feb 2025), fundamentally broadening the class of systems that are analytically tractable.

2. Explicit Representation: Fundamental Matrix Solutions and Solution Formulas

Solution formulas for linear DMDEs with delays rest on the construction of a matrix-valued fundamental solution or 'propagator.' In the continuous case, the fundamental function $Z(\vartheta)$ is defined piecewise over $[(u-1)\sigma,\,u\sigma)$ as

$$Z(\vartheta) = Q_1(0) + Q_2(\sigma) \frac{\vartheta}{1!} + \cdots + Q_{u+1}(u \sigma) \frac{(\vartheta -(u-1)\sigma)^u}{u!}$$

where each $Q_{u+1}(l\delta)$ is computed recursively. The solution is then

$$X(\vartheta) = Z(\vartheta) \Psi(-\sigma) + \int_{-\sigma}^{0} Z(\vartheta - \sigma - s)\Psi'(s)\,ds + \int_{0}^{\vartheta} Z(\vartheta - \sigma - s) G(s)\,ds.$$

In discrete time,

$$X(u) = Z(u)\Psi(-m) + \sum_{r=-m+1}^{0} Z(u-m-r)\Delta\Psi(r-1) + \sum_{r=1}^{u} Z(u-m-r)G(r-1).$$

When $A_0,A_1$ commute, $Q_{u+1}(l\delta)$ reduces to binomial expressions, but otherwise each level must be built up recursively.

This framework covers pure-delay as well as advanced or functional difference-advance equations, allowing one to systematically track how noncommutative algebraic effects—such as the precise ordering of $A_0$ and $A_1$—influence the state evolution (Asadzade et al., 21 Sep 2025, Mahmudov, 2018, Mahmudov, 27 Feb 2025).

3. Fractional and Distributed-Delay DMDEs

Recent advances address DMDEs with both fractional order derivatives and distributed delays. For Caputo/index-$\alpha$ systems,

$(CD_h^\alpha y)(t) = A y(t) + B y(t-h) + f(t),$

solution formulas involve a delayed Mittag-Leffler matrix function, constructed recursively (in index $k$) by $Q_{k+1}(s) = A Q_k(s) + B Q_k(s-h)$, with $Q_1(0)=I$ (Mahmudov, 2018). The solution is then expressed as

$$y(t) = \int_0^t X_{A,B}^{h,\alpha,\alpha}(t-s) f(s)ds,$$

where $X_{A,B}^{h,\alpha,\beta}(t)$ aggregates the effect of A, B, and fractional order $\alpha$ in a convolutive series over delay intervals.

For multi-delay or distributed-delay systems (with polynomially weighted kernels or general density functions), explicit reduction formulas are derived: either by introducing auxiliary 'moment' variables for each term in the kernel expansion, leading to an extended system with only discrete delays (Pulch, 19 Aug 2024), or by quadrature-based reduction (e.g., Gaussian nodes) that approximates the distributed delay as a sum of weighted discrete delays (Cassidy, 12 Oct 2024). Both cases yield high-dimensional but finite systems amenable to standard DDE practices.

In the case of $\mu$-neutral fractional multi-delay systems with noncommutative matrices, solutions in terms of generalized Mittag-Leffler functions (with recursive matrix arguments) are established, together with conditions for existence, uniqueness, and Ulam–Hyers-type stability (Aydin et al., 2022, Mahmudov, 2020).

4. Oscillatory and Stability Behaviors

Oscillatory criteria for delayed DMDEs are derived by casting the original differential system into associated difference equations; for instance, for a system with argument $[t \pm k]$, explicit conditions for oscillation and non-oscillation are tied to the asymptotic behavior of integral transforms involving the system coefficients (Naranjo et al., 17 Jul 2025). Sufficient criteria generalize the classical discrete oscillation theory (e.g., via Erbe–Zhang, Ladas–Philos–Sficas, Győri–Ladas conditions) to the continuous–discrete hybrid case with impulsive jumps.

Stability analysis utilizes fundamental matrix solutions and Lyapunov functionals—Lyapunov matrix constructions for delay-difference equations characterize stability through symmetry, dynamic, and algebraic properties; generalized relations are established even when delays are noncommensurate, by rational approximation (Rocha et al., 2016). For fractional and Hadamard-type delayed systems, Ulam–Hyers stability and unique solvability under Lipschitz (and growth) conditions are proved via contraction mappings in appropriate function spaces (Mahmudov, 2020, Aydin et al., 2022).

Principles of linearized stability for DMDEs (including cases where uniqueness is not a priori guaranteed) are established via mild solutions and variation-of-constants representations, leveraging the structure of the fundamental solution and fixed-point iteration (Nishiguchi, 2022).

5. Numerical Frameworks: Matrix and Spectral Discretization

Discretizing DMDEs for simulation or numerical analysis employs several strategies:

• Matrix Approach: Fractional derivatives and delay operators are discretized using (upper/lower) triangular strip matrices; on a space-time grid, the Kronecker product merges differential operators into a unified matrix system, to be solved as a linear algebraic system (0811.1355). Delay terms are accounted for via shifters (for fixed delays aligning with the grid) and eliminator matrices (to impose boundary or initial conditions). MATLAB routines are developed for automated assembly and solution.
• Spectral Methods and Magnus Expansion: For numerical integration of non-autonomous linear/semilinear DMDEs, the delayed terms are discretized via Chebyshev collocation (yielding high-accuracy representations), mapping the problem to a large but finite-dimensional ODE which is subsequently integrated with a Magnus expansion-based scheme (Arnal et al., 2022). Such methods are effective in handling characteristic multipliers and preserving spectral properties.
• Functional Continuous Runge-Kutta (FCRK): To simulate distributed or multi-delay systems, FCRK methods combine time-stepping with quadrature-based discretization for convolution integrals; the dominant error arises from the less accurate of the two components (either quadrature or time integration). This framework guarantees global error of order $\min(p,q)$, where $p$ and $q$ are the orders of the Runge–Kutta and quadrature components (Cassidy, 12 Oct 2024).
• Reduction to Discrete-Delay DDEs: Distributed kernels are often discretized using quadrature, resulting in equivalent multi-delay DDEs, for which existing DDE solvers can be applied with error control directly related to the quadrature scheme (Pulch, 19 Aug 2024, Cassidy, 12 Oct 2024).

6. Functional and Spectral Properties: Characteristic Matrix Functions and Symmetry

Characteristic matrix functions provide a finite-dimensional analytical method for extracting the spectrum of the (infinite-dimensional) solution operator for periodic DMDEs. By conjugating the monodromy operator with group symmetries (when the DDE is equivariant under a compact group), the characteristic matrix function $\Delta(z)$ encodes eigenvalues (Floquet multipliers) as zeros of $\det \Delta(z)$, reducing stability analysis to root-finding in matrix functions (Wolff, 2022). This is particularly effective for systems with discrete wave solutions, or in delayed feedback stabilization problems with spatio-temporal symmetry.

7. Extensions: Nonlinear, Fractional, and Large-Delay Regimes

The frameworks above extend naturally to:

• Nonlinear DMDEs: Existence, uniqueness, and stability under almost automorphic forcing functions and nonlinearity, using generalizations of classical Massera and bi-almost automorphicity concepts (Coronel et al., 2014).
• Fractional/Neutral-type: Analysis covers Caputo, Hadamard, and $\mu$-neutral cases, using generalized Mittag-Leffler functions for both commutative and noncommutative coefficient matrices (Mahmudov, 2018, Mahmudov, 2020, Aydin et al., 2022).
• Large-delay and multiple-scale analysis: Asymptotic reduction of strongly delayed DMDEs shows convergence to effective partial differential equations (e.g., diffusion equations governing the slow evolution of the envelope), with secular divergences eliminated by solvability conditions (Kozyreff, 2023).

8. Applications and Broader Impacts

Delayed matrix differential equations encompass models in:

• Control theory (delayed feedback, robust stabilization, iterative learning)
• Signal processing (filters with memory)
• Networked systems (information propagation, consensus under delays)
• Mathematical biology (compartmental dynamics, maturation, distributed-latency regulatory feedback)
• Finance (delayed stochastic processes with matrix-valued noise) Interpretation and simulation of these systems benefit from analytical explicit representations, robust stability criteria, and spectrally accurate numerical schemes. The integration of noncommutativity, fractional order, and distributed delay structures substantially widens the scope of real-world systems to which the modern theory of DMDEs applies.

References Table

Key Technique/Concept	Representative Reference(s)	Role/Utility
Recursive fundamental solution $Z$ and $Q_u$	(Asadzade et al., 21 Sep 2025, Mahmudov, 2018)	Explicit solution for noncommutative delays
Matrix approach for discretization	(0811.1355)	Converts fractional/DMDEs to algebraic system
Delayed Mittag-Leffler matrix function	(Mahmudov, 2018, Mahmudov, 2020)	Explicit solution for fractional/neutral DDEs
Characteristic matrix functions & symmetry	(Wolff, 2022)	Finite-dimensional spectral reduction
Spectral/Magnus and FCRK-based numerics	(Arnal et al., 2022, Cassidy, 12 Oct 2024)	High-accuracy simulation of delayed systems
Multi-scale/large-delay reduction	(Kozyreff, 2023)	Asymptotic PDE limits for strong delay
Oscillation and hybrid discrete-continuous	(Naranjo et al., 17 Jul 2025)	Criteria for oscillation in advanced/delayed systems
Lyapunov matrices for delay-difference eqns	(Rocha et al., 2016)	Construction/stability for delay-difference
Polynomial/distributed delay reduction	(Pulch, 19 Aug 2024)	Reduction to equivalent multi-discrete DDE

The current state of the art in delayed matrix differential equations offers a systematic foundation for explicit analysis, robust stability criteria (Lyapunov, spectral, automorphicity-based), high-order simulation, and advanced model reduction. The extension to noncommutative, high-dimensional, and distributed-delay regimes marks a significant generalization, accommodating modern applications in engineering, biology, and networked dynamical systems.