My PhD was on the topic of delay differential equations (DDEs) and, at the time, I wrote lots of utilities to do things like calculate the stability of linearised equations. All that code has been lost to the mists of time but I still have call for it every now and again (e.g., for real-time dynamic substructuring where I need to model the delays in the control loop).

Here is some code in Julia that calculates the stability of a DDE in the form

$$ x'(t) = Ax(t) + Bx(t-\tau) $$

where $A$ and $B$ can be (square) matrices, i.e., systems of DDEs. If you’ve got a nonlinear model, you just need to linearise around an equilibrium first.

The algorithm implemented is that of *Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential Equations*; D. Breda, S. Maset, and R. Vermiglio, SIAM Journal on Scientific Computing 2005 27:2, 482-495. It has the property that it approximates the characteristic roots of the DDE closest to the origin best.

```
using LinearAlgebra
"""
cheb_mesh(T::Type, N::Integer)
cheb_mesh(N::Integer)
Return the Chebyshev nodes of degree `N` on the interval [0, 1].
"""
cheb_mesh(T::Type, N::Integer) = (-cos.((T(π) * (0:N)) ./ N) .+ 1) ./ 2
cheb_mesh(N::Integer) = cheb_mesh(Float64, N)
"""
cheb(T::Type, N::Integer)
cheb(N::Integer)
Return a tuple `(x, D)` where `x` contains the Chebyshev nodes of degree `N` on the interval
[0, 1] using the numerical type `T` and `D` is the Chebyshev differentiation matrix on the
same interval.
## Example - differentiate a sinusoid
julia> (x, D) = cheb(10) ;
julia> D * sinpi.(x) ≈ π * cospi.(x)
true
"""
function cheb(T::Type, N::Integer)
if N == 0
return (ones(T, 1), zeros(T, 1, 1))
else
x = cheb_mesh(T, N)
c = [2; ones(T, N - 1, 1); 2] .* (-1) .^ (0:N)
dx = x .- x'
D = (c * (1 ./ c)') ./ (dx + I)
D[diagind(D)] .-= vec(sum(D; dims=2))
return (x, D)
end
end
cheb(N::Integer) = cheb(Float64, N)
"""
dde_stability(A, B, τ; N=10)
Calculate (a finite number of) the characteristic roots of the linear delay differential
equation
x'(t) = Ax(t) + Bx(t-τ)
where `A` and `B` can be (square) matrices. `N` controls the order of the Chebyshev
polynomial approximation to the infinitesimal generator; higher values of `N` result in a
more accurate approximation. A property of the algorithm used is that characteristic roots
closest to the origin are approximated best.
## Reference
Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential
Equations; D. Breda, S. Maset, and R. Vermiglio, SIAM Journal on Scientific Computing 2005
27:2, 482-495
"""
function dde_stability(A, B, τ; N=10)
@assert size(A) == size(B)
@assert size(A, 1) == size(A, 2)
T = eltype(A) == Int ? Float64 : eltype(A)
n = size(A, 1)
(x, D) = cheb(T, N)
M = kron(-D, Matrix(I / convert(T, τ), n, n))
M[1:n, 1:n] .= A
M[1:n, (n + 1):(end - n)] .= 0
M[1:n, (end - n + 1):end] .= B
return eigen(M)
end
# Tests
# 1. x'(t) = 2x(t) - exp(1)*x(t - 1) has a double root at λ = 1
# 2. x'(t) = -x(t - τ) goes unstable at τ=π/2
# 3. x'(t) = α*A*x(t) + (1-α)*A*x(t-τ)
# A = [
# 50 284 41 23 50 32
# -280 -46 -19 -37 -10 -28
# 35 -1 26 143 35 17
# 5 -31 -139 -22 5 -13
# 20 -16 11 -7 -115 137
# -10 -46 -19 -37 -145 -163
# ]
# with α=0.9 is stable for 3π/(270*(2*0.9 - 1)) ≤ τ ≤ 4π/270
# that is 0.04363323129985824 ≤ τ ≤ 0.046542113386515455
```