Phase and amplitude responses for delay equations using harmonic balance

Robust delay induced oscillations, common in nature, are often modeled by delay-differential equations (DDEs). Motivated by the success of phase-amplitude reductions for ordinary differential equations with limit cycle oscillations, there is now a growing interest in the development of analogous approaches for DDEs to understand their response to external forcing. When combined with Floquet theory, the fundamental quantities for this reduction are phase and amplitude response functions. Here, we develop a framework for their construction that utilises the method of harmonic balance.

Time delays can lead to oscillatory behaviour in many real world systems, as exemplified by laser networks [1], machine dynamics [2], and neural systems [3].These can often be described by delay-differential equations (DDEs) with a stable limit cycle.Such systems do not occur in isolation and are perturbed by external forcing, or interactions with other oscillators.It is therefore vital to be able to quantify the effect these perturbations have on the system behaviour.However, the understanding of the response of DDEs to external forces is challenging due to their infinite dimensionality For weak perturbations a popular method of oscillator reduction for ordinary differential equations (ODEs) is based on formulating phase dynamics along a cycle, and this has recently been extended to account for some notion of distance from cycle using isostable coordinates [4][5][6].For a recent overview see [7].It has been shown that phase-amplitude reduction retaining the phase and slowest decaying amplitude (isostable) of oscillators can capture qualitative changes in stability of phase-locked states occurring under increasing interaction strength in networks of coupled ODE oscillators [8].Extending these results to networks of delayed systems requires the computation of phase and amplitude response to perturbations of limit cycle solutions of DDEs, as well as associated Floquet exponents and eigenfunctions.The phase only reduction has previously been generalised to DDEs in [9,10], and only recently has a phase-amplitude formulation been proposed using a functional analytic perspective by Kotani et al. [11].In both these approaches the practical application of the theory is via the numerical solution of DDEs using time evolution methods.Here, we propose an alternative approach, based upon the method of harmonic balance, that can side-step the numerical challenges associated with time evolution.
In this paper we describe how to arrive at the linear (adjoint) equations for phase and amplitude responses (with appropriate normalisations) using a generalisation of the approach of Novičenko and Pyragas [9].This involves a discretisation of the (infinite dimensional) DDE system as a system of (high dimensional) ODEs and the use of the ODE theory for phase-amplitude reduction to obtain the equations for the DDE phase and amplitude responses in the continuum limit.We further introduce a practical methodology, based upon harmonic balance, to approximate the periodic solutions of the (linear) DDEs describing the response functions in addition to determining the delay induced orbit (from a nonlinear DDE) and its Floquet exponents and eigenfunctions.To illustrate the utility of our approach we compare against two nonlinear DDE models for which analytical results are known.
For limit-cycle oscillators described by ODEs of the form where y ∈ R m and ϵP (t) is a small time-dependent perturbation, we assume that when ϵ = 0 the system has a stable limit cycle y γ (t) with period T .The phase θ = θ(y) and amplitudes ψ i = ψ i (y) can be defined near the limit cycle such that θ = ω = 2π/T and ψi = µ i ψ i where µ i , i = 1, • • • m−1 are Floquet exponents for which Re(µ i ) < 0. The exponent with real part closest to zero is µ 1 := µ which here we assume to be real and small so that perturbations in the direction of the corresponding eigenfunction g 1 (t) := g(t) decay slowly.We assume that all other exponents have Re(µ i ) large and negative so we may retain only ψ 1 := ψ and all other amplitudes may be neglected as they decay much faster.
In the presence of forcing, the standard first order phase-amplitude reduction is given by [12], where T denotes transpose, but higher order corrections may also be included [8,12,13].Here Z(t) and I(t) are the T -periodic phase and amplitude response functions respectively, evaluated on the limit cycle.They quantify the linear response of the phase and amplitude of the oscillator to the perturbation and can be computed as the T periodic solutions of the adjoint equations arXiv:2404.17356v1 [math.DS] 26 Apr 2024 normalised according to where I m is the m × m identity matrix and J := DG(y γ (t)) is the Jacobian of the vector field G evaluated on the limit cycle.Here g(t) is the T -periodic Floquet eigenfunction associated with the exponent µ which satisfies the linear equation [14] ġ(t) = (J − µI m )g(t), (7) and we are free to specify the normalisation, usually making |g(0)| = 1.
In the current work we focus on the computation of phase and amplitude response functions for perturbed limit cycle oscillators described by DDEs of the form where x ∈ R m and τ is a constant delay.We assume that when ϵ = 0, (8) admits a limit cycle solution x γ (t) with period T .Following [9], equation ( 8) is equivalent to where s ∈ [0, τ ] and therefore ξ(τ, t) = x(t − τ ).Next, discretise this as a system of m(N + 1) ODEs by defining x 0 (t) = x(t) and x i (t) = ξ(iτ /N, t) ≈ x(t − iτ /N ) for i = 1, . . ., N .Then where .
The phase adjoint equation for system (11) is therefore of the form (4) where and the normalisation condition is for i = 1, . . ., N , Novičenko and Pyragas [9] observe that in the limit N → ∞, z(t) (the phase response function for the DDE ( 8)) satisfies the adjoint equation with the normalisation Similarly, extending the work of Novičenko and Pyragas [9], the amplitude adjoint equation for system (11) is of the form (5) where I(t) = (q T 0 (t), q T 1 (t), . . ., q T N (t)) T .Observe that for an ODE system, if we write W (t) = e −µt I(t) then W (t) satisfies (4) when I(t) solves (5).Writing W (t) = (w T 0 (t), w T 1 (t), . . ., w T N (t)) T , making the substitution w 0 (t) = w(t) and w i (t) as in (13) with w replacing z, and then taking the limit N → ∞ we conclude that w(t) = w 0 (t) = e −µt q 0 (t) = e −µt q(t) satisfies (14).Therefore the amplitude response for the DDE (8), q(t) satisfies the adjoint equation The normalisation for (16) requires the Floquet eigenfunction ρ(t) for the DDE.For the discretised ODE system (11) the Floquet eigenfunction g(t) = (ρ T 0 (t), ρ T 1 (t), . . ., ρ T N (t)) T satisfies (7).Then h(t) = e µt g(t) satisfies ḣ = Jh (the linearised equation for deviations h from limit cycle).In the limit N → ∞ this has solution h i (t) = h 0 (t − iτ /N ) where and we make the normalisation choice max t |ρ(t)| = 1.We also note that ẋγ (t) is the eigenfunction of the linearised system with µ = 0. Since for i = 1, . . ., N , q i (0) = w i (0) and ρ i (0) = h i (0) we see that Taking the limit N → ∞ the normalisation condition for ( 16) is given by [15] q T (0)ρ(0) + e −µτ 0 We now turn to the problem of how to efficiently numerically calculate the quantities z(t), q(t), ρ(t) and µ as well as the limit cycle.Improvements in accuracy and computational speed can be made over numerical DDE solvers by using the method of harmonic balance (also known as the Fourier-Galerkin method).The harmonic balance method is widely used for analysing the periodic solutions of nonlinear differential equations, often in a mechanical setting, see e.g., [16].The periodic solution in R m is expressed as a Fourier series truncated at M Fourier modes.To determine the m(2M + 1) unknown Fourier coefficients, an algebraic zero problem of size m(2M + 1) is formulated by considering the solution sampled at 2M + 1 time points.The method can also be used for determining periodic solutions of delay differential equations, see e.g., [17].Simmendinger et al. [18] also use Fourier expansions to approximately calculate Floquet exponents and eigenfunctions for delay differential equations.
For the DDE (8) first consider the T -periodic orbit x γ (t) ∈ R m with truncated Fourier series representation where a p ∈ C m with a −p = a * p and * denotes complex conjugation.There are m(2M + 1) + 1 real unknowns to determine, namely the Fourier coefficients a 0 , a 1 , . . ., a M and the period T .Sample x γ (t) at 2M + 1 time instants t n = nT /(2M + 1), n = −M, . . ., 0, . . ., M and introduce the notation T , (21) Then X = (S ⊗ I m )A where S is the symmetric Vandermonde matrix with [S] np = e 2πinp/(2M +1) , n, p ∈ {−M, . . ., 0, . . .M }, and ⊗ is the Kronecker product.Sampling (8) on the limit cycle results in the system of m(2M + 1) nonlinear algebraic equations for the components of X and we note that S −1 = S * /(2M + 1).To fix an origin of the one parameter family of periodic orbits we set the nth component of x γ to have a vanishing derivative at the origin, or equivalently: where e n ∈ R m is a canonical vector for the nth direction.
The zero problem given by ( 23)-( 24) can then be solved numerically for X, for example using an optimisation approach based upon the Levenberg-Marquardt algorithm.
To determine stability we consider again (18) describing a small T -periodic perturbation ρ(t) such that x(t) = x γ (t) + e µt ρ(t) for µ ̸ = 0. Sampling (18) at times t n and using a truncated Fourier series representation of ρ(t) as ρ(t) = M p=−M b p e iωpt yields the linear algebraic system where , For nontrivial solutions of (25) we require that E(µ) = det(M(µ)) = 0.The approximated Floquet exponents µ j are solutions of E(µ) = 0 with corresponding sampled eigenfunctions R j .
For stable limit cycles Re(µ j ) < 0 for all nonzero Floquet exponents.Here we consider the case where the largest nontrivial Floquet exponent µ 1 := µ is real.
The phase and amplitude responses can be computed as appropriately normalised solutions of the adjoint equation ( 16) (with µ = 0 giving phase response).Sampling ( 16) at times t n and using a truncated Fourier series representation of q(t) as q(t) = M p=−M c p e iωpt yields another linear algebraic system where , and the relationship to the Fourier coefficients is given by Q The systems of linear equations ( 25) and ( 26) can be solved faster than using DDE solvers to determine the eigenfunctions and response functions, also recalling the normalisations (15) has a limit cycle x γ (t) = cos(t).For small values of δ the forms for µ, ρ(t), z(t), and q(t) can be found analytically [11].For ϵ = 0.05 we plot these analytical expressions against the result of the harmonic balance method with truncations at M = 20 in Figure 1.The harmonic balance approach faithfully reproduces the eigenfunction and both phase and amplitude response curves.We also find that harmonic balance improves on the accuracy in addition to the speed of these computations over standard DDE solvers (not shown).As a second (nonscalar) example we consider a simplified model for cortico-thalamic EEG rhythms [19], dx(t) dt = y(t), For α = −0.039,β = −0.4,γ = −2.0,δ = −10.0 and τ = 8.0, centre manifold reduction can be used to analytically approximate the y component of the phase response [20].Kotani et al. [10] match this approximation to the PRC calculated by numerically integrating the adjoint equation and by direct perturbation.Using the harmonic balance approach we calculate the limit cycle and the phase response to perturbations in both the x and y components as in Figure 2(A) and (B) respectively, showing that we also obtain agreement with the analytical expression from [20].We further use harmonic balance to determine the Floquet exponent as µ = −0.00296and to compute the corresponding eigenfunction and amplitude response function as shown in Figure 2(C) and (D) respectively.In this case, since the Floquet exponent is so close to zero one would expect only a very slow decay to the limit cycle, which would require long computation times using a direct numerical method based upon time evolution.This issue does not arise when using the harmonic balance approach, whose accuracy can be further improved with increasingly larger choices for the truncation parameter M .For simplicity we have restricted to systems with a single time delay, however it is straightforward to extend the approaches given here to systems with multiple time delays.Extensions to consider systems where the Floquet exponent with largest real part is complex are natural and allow for the treatment for a wider range of DDEs.Furthermore, higher order correction terms to the phase and isostable response curves can be computed for ordinary differential equation systems with limit cycle solutions [12,14,21,22] by solving inhomogeneous linear equations.Analogues of these equations for delay differential equations can also be computed and solved using harmonic balance to obtain more accurate equations for the dynamics of the phase and amplitude of solutions of delay systems.This allows for the exploration of synchronisation and phase-locked state dynamics in networks of delay systems, with and without additional interaction delays, extending previous work [8] to delay systems.These extensions will be reported on at length elsewhere.

FIG. 1 .
FIG. 1.The (A) limit cycle, (B) phase response curve (C) Floquet eigenfunction, and (D) amplitude response curve for model (27) plotted against time.The analytical curves from [11] are shown with red circles and numerical results using the harmonic balance method are shown with blue lines.Here δ = 0.05 and M = 20.

FIG. 2 .
FIG. 2. The (A) limit cycle, (B) phase response curve (C) Floquet eigenfunction, and (D) amplitude response curve for model (28) plotted against phase.The analytical curve from [20] is shown with yellow circles and the x and y components from the harmonic balance method are shown with a blue and red line respectively.Parameter values are α = −0.039,β = −0.4,γ = −2.0,δ = −10.0 and τ = 8.0 and M = 20.