Bifurcation Analysis of Systems with Delays: Methods and Their Use in Applications
BBifurcation Analysis of Systems withDelays:Methods and Their Use in Applications
Bernd Krauskopf * and Jan Sieber † * Department of Mathematics, University of Auckland, Private Bag92019, Auckland 1142, New Zealand † College of Engineering, Mathematics and Physical Sciences, Uni-versity of Exeter, Exeter EX4 4QF, United Kingdom (preprint)
This chapter presents a dynamical systems point of view of the study ofsystems with delays. The focus is on how advanced tools from bifurcationtheory, as implemented for example in the package DDE-BIFTOOL, canbe applied to the study of delay differential equations (DDEs) arising inapplications, including those that feature state-dependent delays. We discussthe present capabilities of the most recent release of DDE-BIFTOOL. Theyinclude the numerical continuation of steady states, periodic orbits and theirbifurcations of codimension one, as well as the detection of certain bifurca-tions of codimension two and the calculation of their normal forms. Twolonger case studies, of a conceptual DDE model for the El Niño phenomenonand of a prototypical scalar DDE with two state-dependent feedback terms,demonstrate what kind of insights can be obtained in this way.
Systems with delays arising in applications come in many different forms. From a generalperspective a DDE is an ordinary differential equation (ODE) with a number of terms thatfeature delays. When the delays are zero, or parameters multiplying such terms are zero,then the DDE reduces to the underlying ODE, that is, a finite-dimensional dynamicalsystem. When delays are present, on the other hand, one is dealing with an actual DDEand, hence, with an infinite-dimensional dynamical system. As for ODEs, the task is todetermine the possible dynamics of a given DDE as a function of its parameters. In otherwords, what is called for is a bifurcation analysis of the DDE that unveils the divisionof parameter space into regions of different behavior. In spite of this difference in the1 a r X i v : . [ m a t h . D S ] S e p imension of the phase space, the bifurcation theory of DDEs is effectively that of ODEsin the sense that the same bifurcations of equilibria and periodic orbits arise in bothcases. The complicating issue is that equilibria, periodic orbits and their bifurcationsof a given DDE ‘live’ in an infinite-dimensional space. As is the case for ODEs, thisrequires specialized numerical tools for finding and tracking invariant objects and theirbifurcations.As we will demonstrate, such advanced tools are available today. We focus here on thecapabilities as implemented in the package DDE-BIFTOOL — a Matlab/octave (MAT-LAB, 2018; Eaton et al., 2017) compatible library for performing numerical bifurcationanalysis of DDEs of different types. DDE-BIFTOOL uses a numerical continuationapproach, originally implemented by Engelborghs et al. (2000b, 2001); Samaey et al.(2002), and is currently accessible and maintained at sourceforge.net/projects/ddebiftool (Sieber et al., 2015). Its capabilities are a subset of those of commonly usedtools for ODEs and maps, such as AUTO (Doedel et al., 1999; Doedel, 2007), MATCONT(Dhooge et al., 2003; Govaerts, 2000) or COCO (Dankowicz and Schilder, 2013). Thebifurcation analysis tool knut offers an alternative, stand-alone implementation (in C++)of many of the methods used in DDE-BIFTOOL; see (Roose and Szalai, 2007).In contrast to the tools for ODEs, the package DDE-BIFTOOL permits differentialequations with a finite number of discrete delays in their arguments. More precisely, itconsiders differential equations of the form 𝑀 𝑥 ′ ( 𝑡 ) = 𝑓 ( 𝑥 ( 𝑡 ) , 𝑥 ( 𝑡 − 𝜏 ) , … , 𝑥 ( 𝑡 − 𝜏 𝑑 ) , 𝑝 ) , where (1) 𝑀 ∈ ℝ 𝑛 × 𝑛 and 𝑓 ∶ ℝ 𝑛 ×( 𝑑 +1) × ℝ 𝑛 𝑝 → ℝ 𝑛 ,which is a DDE with 𝑑 + 1 delays (where 𝜏 = 0 is included in the list of discrete delays).Further, there are 𝑛 𝑝 system parameters 𝑝 ∈ ℝ 𝑛 𝑝 . We call the space ℝ 𝑛 , in which 𝑥 ( 𝑡 ) lives,the physical space for (1) (not to be confused with the infinite-dimensional phase space ofthe DDE as defined below). The matrix 𝑀 on the left-hand side is most commonly (andby default) equal to the identity matrix, which corresponds to the case of ‘standard’ DDEs.Different choices of 𝑀 are used to define other types of DDEs, including equations thatare neutral (featuring delayed derivatives), are of differential algebraic form, or are ofmixed type with both delayed and advances terms.DDE-BIFTOOL distinguishes two types of DDEs, depending on the nature of thediscrete delays 𝜏 , … , 𝜏 𝑑 : DDEs with constant delays and DDEs with state-dependent delays. In the case of constant-delay DDEs the delays 𝜏 , … , 𝜏 𝑑 need to be part of thevector 𝑝 of system parameters. As part of the setup, the user has to specify the list of 𝑑 indices of 𝑝 that correspond to the delays 𝜏 , … , 𝜏 𝑑 . Alternatively, when dealing withstate-dependent delays one may specify the delays 𝜏 , … , 𝜏 𝑑 as functions of current ordelayed states of 𝑥 . More precisely, the user has to specify the number 𝑑 of delays andthe functions 𝜏 𝑓𝑗 ∶ ℝ 𝑛 × 𝑗 × ℝ 𝑛 𝑝 → ℝ for 𝑗 = 1 , … , 𝑑 , where then 𝜏 𝑗 = 𝜏 𝑓𝑗 ( 𝑥 , … , 𝑥 𝑗 −1 , 𝑝 ) for 𝑗 = 1 , … , 𝑑 , with (2) 𝑥 = 𝑥 ( 𝑡 ) , 𝑥 𝑗 = 𝑥 ( 𝑡 − 𝜏 𝑗 ) for 𝑗 = 1 , … , 𝑑 , such that 𝑥 ′ ( 𝑡 ) = 𝑓 ( 𝑥 , … , 𝑥 𝑑 , 𝑝 ) . (3)2his recursive definition permits arbitrary levels of nesting such that delays may dependon delayed states.We focus here on ‘standard’ DDEs with constant and state-dependent delays. We firstbriefly discuss their relevant properties as dynamical systems. Subsequently, we presentthe tasks of bifurcation analysis and then discuss how they are performed and set up inpractice in DDE-BIFTOOL; here, we use a constant-delay DDE for the inverted pendu-lum with delayed control as the illustrating example throughout. We further illustratethe overall capabilities with two longer case studies: (1) a conceptual DDE model forthe El Niño Southern Oscillation (ENSO) system with negative delayed feedback andperiodic forcing, where the delay is initially constant and then state dependent; and (2)a prototypical scalar DDE with two state-dependent feedback terms that features onlytrivial dynamics in the absence of state dependence. While DDEs with discrete delays are the most common type of DDEs considered forpractical implementation of numerical methods, the underlying mathematical theory doesnot distinguish between discrete and, e.g., distributed delays. The general theory permitsgeneral functionals on the right-hand side, of the form ̃𝑓 ∶ 𝐶 ([− 𝜏 max , ℝ 𝑛 ) × ℝ 𝑛 𝑝 → ℝ 𝑛 . Here 𝐶 𝑘 ([− 𝜏 max , ℝ 𝑛 ) (or 𝐶 𝑘 for short) is the space of 𝑘 times continuously differentiablefunctions — the history segments — on the interval [− 𝜏 max , , where 𝜏 max is an upperbound for the delays; in particular, 𝐶 is the space of continuous functions on [− 𝜏 max , with values in ℝ 𝑛 . The general DDE (also called a functional differential equation, FDE)then has the form 𝑥 ′ ( 𝑡 ) = ̃𝑓 ( 𝑥 𝑡 , 𝑝 ) , (4)for the standard case that 𝑀 in (1) is the identity matrix. One looks for solutions 𝑥 ( 𝑡 ) ∈ ℝ 𝑛 with 𝑡 ∈ [− 𝜏 max , 𝑡 end ] of (4), and the solution 𝑥 𝑡 with subscript 𝑡 is the current historysegment in 𝐶 , that is, 𝑥 𝑡 ∶ [− 𝜏 max , → ℝ 𝑛 𝜃 ↦ 𝑥 ( 𝑡 + 𝜃 ) .For the DDE (1) with discrete delays the functional ̃𝑓 has the form ̃𝑓 ( 𝑥, 𝑝 ) = 𝑓 ( 𝑥 (0) , 𝑥 (− 𝜏 ) , … , 𝑥 (− 𝜏 𝑑 ) , 𝑝 ) for 𝑥 ∈ 𝐶 . In case the delays are state dependent, the 𝜏 𝑗 are defined as described by (2)when setting 𝑡 = 0 . For DDEs with constant delays the textbooks by Hale and Verduyn Lunel (1993) andDiekmann et al. (1995) develop the necessary theory that permits one to consider DDEs3f the general type (4) (and, hence, (1)) as regular dynamical systems on the phase space 𝐶 of continuous functions over the (maximal) delay interval; these DDEs are referredto as abstract ODEs by Diekmann et al. (1995). For 𝑡 = 0 one has to provide an initialhistory segment 𝜙 ∈ 𝐶 , and then at each time 𝑡 ≥ the current state is the function 𝑥 𝑡 ,which is also in 𝐶 ; in particular, 𝑥 = 𝜙 . The textbooks show that the map 𝑋 ∶ [0 , ∞) × 𝐶 → 𝐶 ( 𝑡, 𝜙 ) ↦ 𝑥 𝑡 , which maps time 𝑡 and initial value 𝜙 to the history segment 𝑥 𝑡 at time 𝑡 of the solutionof the DDE, is as regular with respect to its argument 𝜙 as the right-hand side ̃𝑓 of (1) iswith respect to its first argument. For example, if 𝜙 ↦ ̃𝑓 ( 𝜙, 𝑝 ) is 𝓁 times continuouslydifferentiable, then so is 𝜙 ↦ 𝑋 ( 𝑡 ; 𝜙 ) . Consequently, the general theory transfers manyresults of the bifurcation theory for ODEs to the case of DDEs with constant delays. Inparticular, the solution map 𝑋 is eventually compact, which implies that local centermanifolds in equilibria and periodic orbits of DDEs are finite-dimensional and as regularas the right-hand side ̃𝑓 . Therefore, the local bifurcation theory of DDEs with a finitenumber of constant delays is identical to the local bifurcation theory of ODEs. For DDEs with state-dependent delays the claim that their local bifurcation theory isidentical to the theory of ODEs is not fully resolved. A review of well established resultsand an exposition of the obstacles that one initially faces are described in the review byHartung et al. (2006). Even assuming that the state-dependent delays are always boundedwithin an interval [0 , 𝜏 max ] , the space of continuous 𝐶 is not a suitable phase space,since no local uniqueness of solutions to initial-value problems can be guaranteed. Thedifficulty lies in the fact that for state-dependent delays the functional ̃𝑓 ∶ 𝐶 → ℝ 𝑛 inthe right-hand side of the DDE is not continuously differentiable (or locally Lipschitzcontinuous), even if all coefficients 𝑓 and 𝜏 𝑓𝑗 are smooth. In fact, for general DDEs ofthe form (4) the assumption that ̃𝑓 is continuously differentiable with respect to its firstargument is not satisfied when the delays are state dependent. In fact, the assumptionof regularity of 𝜙 ↦ ̃𝑓 ( 𝜙, 𝑝 ) being satisfied could be considered as the general propertyunderlying and, hence, defining the constant-delay case of the theory.Walther (2003) observes that functionals ̃𝑓 involving state-dependent delays satisfy aweaker regularity condition, which could be called mild differentiability. For these typesof functionals Walther (2003) proves that for history segments 𝜙 within the manifold of compatible initial conditions, defined as 𝜙 ∈ 𝐶 ∶= { 𝜙 ∈ 𝐶 ∶ 𝜙 ′ (0) = ̃𝑓 ( 𝜙 ) } ,a unique solution 𝑋 ( 𝑡 ; 𝜙 ) of the DDE exists and is also in 𝐶 . Moreover, for each time 𝑡 ≥ , the solution map 𝑋 ∶ 𝐶 → 𝐶 𝜙 ↦ 𝑋 ( 𝑡 ; 𝜙 ) (5)4s continuously differentiable once, which means that 𝑋 meets the conditions for basicstability theory. For example, the principle of linearized stability holds for equilibria andperiodic orbits (Skubachevskii and Walther, 2006; Mallet-Paret and Nussbaum, 2011b).Moreover, local center manifolds exist, are finite-dimensional and are continuouslydifferentiable once (Stumpf, 2011).The results for the solution map 𝑋 cannot be generalized to higher degrees of continu-ous differentiability; hence, 𝑋 is not sufficiently regular to support all aspects of localbifurcation theory. However, there exist some results on higher-order differentiability forDDEs with state-dependent delays. First, solutions of periodic boundary-value problemsfor state-dependent delays can be reduced to finite-dimensional algebraic systems ofequations that are as regular as the coefficients, such as 𝑓 and 𝜏 𝑓𝑗 in (2) and (3). Thus,all computations performed during numerical bifurcation analysis of equilibria, periodicorbits or their local bifurcations can be performed as expected and depend smoothly ontheir data. This includes the standard tasks of continuation of solutions using Newtoniterations and pseudo-arclength continuation, or branching off at singularities. Similarly,all computations performed during normal form analysis are feasible (Sieber, 2012, 2017).Furthermore, Krisztin (2006) checked that the techniques used for obtaining 𝓁 > timesdifferentiable local unstable manifolds of equilibria (Krisztin, 2003) are also applicableto local center manifolds of equilibria.Thus, the results by Krisztin (2006) strongly suggest that local center manifolds inDDEs with state-dependent delays are differentiable as often as the coefficients 𝑓 and 𝜏 𝑓𝑗 , even though the solution map 𝑋 of (5) is not. For this reason, while this is notfully resolved, the local bifurcation theory of DDEs with state-dependent delays is stillexpected to be identical to the theory for ODEs. Indeed, this claim is supported by alltheoretical results this far, as well as by numerical investigations such as the one presentedin Section 5. The general theory of Section 2 implies that numerical bifurcation analysis should allowone to perform a range of tasks for equations of type (1), similar to those arising in thebifurcation analysis of ODEs. More specifically, the local bifurcations of (standard)DDEs with both constant and state-dependent delays are the same as those one findsin ODEs and can, hence, be found in standard textbooks such as (Guckenheimer andHolmes, 1983; Govaerts, 2000; Kuznetsov, 2013). We do not present or review here thisextensive theory but rather focus on the typical tasks required for the bifurcation analysisof a given DDE (or ODE), which include:1. continuation of equilibria in a single system parameter;2. linear stability analysis at equilibrium points, that is, finding the (leading) eigenval-ues of their linearization; 5. detection of codimension-one bifurcations of equilibria and their continuation ascurves in two system parameters; in generic systems, these are the saddle-node(or fold) bifurcation and the Hopf bifurcation, while in the presence of additional(symmetry) properties they include the transcritical bifurcation and the pitchforkbifurcation;4. detection of codimension-two bifurcations of equilibria; in generic systems, theseinclude the saddle-node Hopf, Hopf-Hopf, cusp and degenerate Hopf bifurcations;5. normal form analysis of generic codimension-one and codimension-two bifurca-tions of equilibria; this includes, for example, computing the Lyapunov coefficientfor the Hopf bifurcation, which determines whether the bifurcation is supercriticalor subcritical, and branching off to secondary solution or bifurcation branches;6. continuation of periodic orbits in a single system parameter (with automatic adjust-ment of the period);7. linear stability analysis of periodic orbits, that is, determining their (leading) Floquetmultipliers;8. detection of codimension-one bifurcations of periodic orbits and their continuationas curves in two system parameters; in generic systems, these are the saddle-node(or fold) bifurcation of periodic orbits, the period-doubling bifurcation and thetorus (or Neimark-Sacker) bifurcation;9. identification and continuation of connecting orbits between equilibria in a suitablenumber of system parameters (depending on the dimensions of the respecitve stableand unstable eigenspaces of the involved equilibria);10. computation of unstable manifolds of equilibria and periodic orbits with a singleunstable direction to find, for example, certain invariant tori and global bifurcations.Continuation tasks require formulating an algebraic system of equations of the form 𝐺 ( 𝑦 ) = 0 , where 𝐺 ∶ dom( 𝐺 ) → range( 𝐺 ) (6)is differentiable and the nullspace ker 𝐺 ′ ( 𝑦 ) is one-dimensional in solutions 𝑦 ∈ dom 𝐺 .For the contination of equilibria and their bifurcations we will have dom( 𝐺 ) = ℝ 𝑁 +1 and range( 𝐺 ) = ℝ 𝑁 for some problem dependent 𝑁 ∈ ℕ . For the continuation of periodicorbits and their bifurcations 𝐺 will map between infinite-dimensional spaces. In thesecases a discretized problem 𝐺 d ∶ ℝ 𝑁 +1 → ℝ 𝑁 has to be constructed, where 𝑁 dependson the number of mesh points, and may be increased to improve accuracy. Solution lociof system (6) are generically curves (branches), which are tracked by the continuationalgorithm (Doedel et al., 1999). Within DDE-BIFTOOL, the user specifies two pointsnear each other on the respective curve of interest; the (pseudo-arclength) continuationalgorithm generates a predictor by extrapolating the secant through the (last) two pointson the curve and then uses Newton iteration to correct the prediction and, hence, find the6 a) (b) Figure 1:
Sketch of an inverted pendulum on a cart subject to a control force 𝐹 with its equationin rescaled time units √ 𝐿 ∕ 𝑔 (a), and the linear stability chart in the ( 𝑎, 𝑏 ) -plane for PD control 𝐹 ( 𝑡 ) = 𝑎𝜃 ( 𝑡 − 𝜏 ) + 𝑏𝜃 ′ ( 𝑡 − 𝜏 ) and delay 𝜏 = √ (b). next point along the curve; this predictor-corrector step is repeated until a sufficient pieceof the curve defined by (6) has been computed.Stability computation tasks for DDEs require the creation of a matrix eigenvalueproblem of an approximating, sufficiently large system of ODEs. Normal form analysis forequilibrium bifurcations in DDE-BIFTOOL computes explicit normal form coefficientsfor codimension-one and codimension-two bifurcations, which determine the dynamicsclose to the bifurcation according to the textbook by Kuznetsov (2013). The normal formanalysis is also used to construct predictors for starting the continuation of secondarysolution branches that emerge from the respective bifurcation. We will proceed to explain how these different tasks are performed by DDE-BIFTOOL. Todemonstrate how this works in practice, we will use throughout this section the exampleof a simple DDE model for balancing an inverted pendulum with delayed feedback,described by Sieber and Krauskopf (2004a,b), 𝜃 ′′ ( 𝑡 ) = sin 𝜃 ( 𝑡 ) − 𝐹 ( 𝑡 ) cos( 𝑡 ) , where 𝐹 ( 𝑡 ) = 𝑎𝜃 ( 𝑡 − 𝜏 ) + 𝑏𝜃 ′ ( 𝑡 − 𝜏 ) . (7)The dependent variable 𝜃 ( 𝑡 ) is the angle by which the (approximately mass-less) pendulumdeviates from the upright position. The term 𝐹 ( 𝑡 ) is the feedback force exerted by movingthe base of the pendulum, such as the cart sketched in fig. 1(a), to achieve upright balancing.In (7), the feedback is of proportional-plus-derivative type, where the correcting forcedepends linearly on the angle 𝜃 and on the angular velocity 𝜃 ′ ; one speaks of PD control.Time in (7) is in units of the intrinsic time scale √ 𝑔 ∕ 𝐿 of the pendulum, where 𝐿 is thelength of the pendulum and 𝑔 is the gravitational acceleration. The delay 𝜏 models areaction delay relative to the intrinsic time scale of the pendulum; hence, varying the delayis similar to changing the length of the pendulum, where a shorter pendulum correspondsto a longer delay after scaling. 7ewriting (7) as as a first-order system gives 𝑥 ′1 ( 𝑡 ) = 𝑥 ( 𝑡 ) , 𝑥 ′2 ( 𝑡 ) = sin 𝑥 ( 𝑡 ) − cos 𝑥 ( 𝑡 ) [ 𝑎𝑥 ( 𝑡 − 𝜏 ) + 𝑏𝑥 ( 𝑡 − 𝜏 ) ] . (8)Hence, the physical space for 𝑥 ( 𝑡 ) = ( 𝑥 ( 𝑡 ) , 𝑥 ( 𝑡 )) is ℝ , the parameter space for 𝑝 =( 𝑝 , 𝑝 , 𝑝 ) = ( 𝑎, 𝑏, 𝜏 ) is ℝ , and the right-hand side 𝑓 is 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = [ 𝑥 sin 𝑥 − cos 𝑥 ( 𝑝 𝑥 + 𝑝 𝑥 )] (9)in the formulation (1) of DDE-BIFTOOL, where 𝑀 is the identity in ℝ . We remark that,even though the right-hand side 𝑓 does not depend on 𝜏 , the delay 𝜏 must be included inthe parameter vector 𝑝 (as 𝑝 ) for 𝑓 .We observe that the right-hand side 𝑓 has the reflection symmetry 𝑓 (− 𝑥 , − 𝑥 , 𝑝 ) = − 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) . Therefore, system (8) has two types of solutions: symmetric solutions, which are invari-ant under this reflection symmetry, and symmetrically related pairs of non-symmetricsolutions. As we will see later on, distinguishing the two types is of practical relevancebecause the control force 𝐹 ( 𝑡 ) is applied by adjusting the position of the pendulum at itsbase, that is, by pushing the cart on which it is mounted in fig. 1(a). For non-symmetricsolutions the average of 𝐹 ( 𝑡 ) is non-zero, which implies that the cart accelerates away inone direction; in particular, the cart position, which we refer to as 𝛿 ( 𝑡 ) in fig. 2 below, isunbounded. For symmetric solutions, on the other hand, the long-term average of theforce 𝐹 ( 𝑡 ) on the cart is zero, which means that the cart position 𝛿 ( 𝑡 ) may be bounded.Thus, only symmetric solutions correspond to successful balancing of the pendulum.Mathematically, the reflection symmetry introduces non-generic symmetry-breakingbifurcations, including pitchfork bifurcations, which require special treatment withinDDE-BIFTOOL. The continuation of equilibria for DDEs is identical to that for ODEs: equilibria 𝑥 eq of 𝑥 ′ ( 𝑡 ) = 𝑓 ( 𝑥 ( 𝑡 ) , 𝑥 ( 𝑡 − 𝜏 ) , … , 𝑥 ( 𝑡 − 𝜏 𝑑 ) , 𝑝 ) are given by the nonlinear equation 𝐺 eq ( 𝑦 eq ) = 𝑓 ( 𝑥 eq , … , 𝑥 eq , 𝑝 ) = 0 , which is of the general form (6) where the dimension 𝑁 = 𝑛 is the dimension of thephysical space. One typically chooses the unknowns 𝑦 eq = ( 𝑥 eq , 𝑝 𝑖 ) for continuation ofequilibria in one of the system parameters 𝑝 𝑖 .The pendulum problem (8) has the upright position 𝑥 = (0 , as its reflection-invarianttrivial equilibrium solution. fig. 1(b) shows the line (dashed) for fixed 𝑏 = 1 . and varying 𝑎 in the ( 𝑎, 𝑏 ) -plane for 𝜏 = √ , which corresponds to the single-parameter family oftrivial equilibria (0 , , 𝑎 ) . 8 .3 Linear stability analysis of equilibria The natural next step is the computation of the linear stability for each equilibrium ( 𝑥 eq , 𝑝 ) along the computed branch, which is determined by the stability of the DDE, linearizedin the equilibrium (including 𝜏 = 0 into our list of discrete delays), 𝑀 𝑥 ′ ( 𝑡 ) = 𝑑 ∑ 𝑗 =0 𝐴 𝑗 𝑥 ( 𝑡 − 𝜏 𝑗 ) , where 𝐴 𝑗 = 𝜕𝑓𝜕𝑥 𝑗 ( 𝑥 eq , … , 𝑥 eq , 𝑝 ) ∈ ℝ 𝑛 × 𝑛 . (10)The stability of the origin for the linear DDE (10) is determined by the real parts of thespectrum of its infinitesimal generator , defined by [ 𝑥 ]( 𝜃 ) = 𝑥 ′ ( 𝜃 ) with domain 𝐷 ( ) = { 𝑥 ∈ 𝐶 ∶ 𝑀 𝑥 ′ (0) = 𝑑 ∑ 𝑖 =0 𝐴 𝑖 𝑥 (− 𝜏 𝑖 ) } .The eigenvalue problem 𝜆𝑥 = 𝑥 in the infinite-dimensional space 𝐷 ( ) is equivalentto the 𝑛 -dimensional eigenvalue problem for the characteristic matrix Δ( 𝜆 ) ∈ ℝ 𝑛 × 𝑛 , Δ( 𝜆 ) 𝑣 = 0 , where Δ( 𝜆 ) = 𝜆𝑀 − 𝑑 ∑ 𝑖 =0 𝐴 𝑖 exp(− 𝜆𝜏 𝑖 ) (11)(Kaashoek and Verduyn Lunel, 1992). The problem of finding the right-most eigenvalues 𝜆 of (or Δ ) is fundamentally different from an ODE stability problem, since typically has infinitely many eigenvalues. However, for the most common case where 𝑀 is theidentity, the infinitesimal generator has at most finitely many eigenvalues to the rightof any vertical line in the complex plane. For its linear stability analysis DDE-BIFTOOLdoes not solve det Δ( 𝜆 ) = 0 , but instead discretizes the eigenvalue problem 𝜆𝑥 = 𝑥 toobtain an eigenvalue problem for a pair of large matrices. The approach by Breda et al.(2005) is to discretize the ODE boundary-value problem 𝑥 ′ ( 𝜃 ) = 𝜆𝑥 ( 𝜃 ) , 𝑀 𝑥 ′ (0) = 𝑑 ∑ 𝑖 =0 𝐴 𝑖 𝑥 (− 𝜏 𝑖 ) ,on the interval [− 𝜏 max , with an 𝑚 th-order pseudospectral approximation (e.g., Cheby-shev polynomials) for 𝑥 ∶ [− 𝜏 max , → ℂ 𝑛 to obtain a matrix eigenvalue problem ofdimension ( 𝑚 + 1) 𝑛 . Engelborghs and Roose (2002) discretize (10) by using a linearmultistep ODE solver with a small time step ℎ = 𝜏 max ∕ 𝑚 and express the condition thatthe approximate solution after a single time step, 𝑥 ℎ , satisfies 𝑥 ℎ = 𝜇𝑥 on the uniformgrid with step ℎ on [− 𝜏 max , . This is also an eigenvalue problem for large matrices in 𝜇 ,from which the eigenvalues 𝜆 are obtained by the relation 𝜆 = (log 𝜇 )∕ ℎ . Both optionsare available in DDE-BIFTOOL, which automatically refines the discretization if thedesired right-most eigenvalues are not accurate up to a specified tolerance.For the pendulum model’s upright equilibrium 𝑥 = (0 , the characteristic matrix Δ( 𝜆 ) has the form Δ( 𝜆 ) = [ 𝜆 −1 𝑎 e − 𝜆𝜏 − 1 𝜆 + 𝑏 e − 𝜆𝜏 ] . 𝑎 from for 𝑏 = 1 . and 𝜏 = √ , one detects a change of linearstability at the points 𝑃 and 𝐻 in fig. 1(b). At the point 𝑃 (where 𝑎 = 1 ) a real eigenvaluecrosses the imaginary axis from the right to the left half of the complex plane. At thispoint 𝑃 , the trivial equilibrium gains linear stability (under increasing 𝑎 ) in a subcriticalpitchfork bifurcation (owing to the reflection symmetry), and a family of non-symmetricequilibria of the form ( 𝑥 , 𝑥 , 𝑎 ) = ( 𝑥 , , sin 𝑥 ∕( 𝑥 cos 𝑥 )) branches off. At the point 𝐻 (at 𝑎 ≈ 3 . ), the linear DDE has a pair of complex eigenvalues ±i 𝜔 ≈ ±0 . 𝜋 (crossingfrom left to right for increasing 𝑎 ) such that the upright equilibrium destabilizes in aHopf bifurcation. Consequently, for 𝑏 = 1 . and 𝜏 = √ the upright position of thependulum is linearly stable when 𝑎 ∈ (1 , . (that is, between the points 𝑃 and 𝐻 infig. 1). Once a bifurcation of an equilibrium (which will typically be of codimeonsion one)has been detected one may either branch off to follow another branch (of equilibria orperiodic orbits), or continue the bifurcation itself in additional parameters — as a curvein two chosen system parameters when it is indeed of codimension one. DDE-BIFTOOLsupports continuation of the generic codimension-one bifurcations, the saddle-node (orfold) and the Hopf bifurcation (Engelborghs and Roose, 1999). The continuation isimplemented for nonlinear systems that extend the nonlinear equation 𝐺 eq ( 𝑦 eq ) = 0 withthe eigenvalue problem (11) for the critical eigenvalue. For the fold of equilibria theeigenvalue 𝜆 is such that the unknowns are 𝑦 feq = ( 𝑥 feq , 𝑣 feq , 𝑝 𝑖 , 𝑝 𝑗 ) ∈ ℝ 𝑛 +2 (with twofree parameters 𝑝 𝑖 and 𝑝 𝑗 ). For the Hopf bifurcation the frequency 𝜔 𝐻 is unknown and theeigenvector 𝑣 𝐻 is complex such that we have the unknowns 𝑦 𝐻 = ( 𝑥 𝐻 , 𝑣 𝐻 , 𝜔 𝐻 , 𝑝 𝑖 , 𝑝 𝑗 ) ∈ ℝ 𝑛 +3 (counting the complex vector 𝑣 𝐻 as two real vectors of length 𝑛 ). The nonlinearequations for the fold and Hopf bifurcations are 𝐺 feq ( 𝑦 feq ) = ⎡⎢⎢⎣ 𝑓 ( 𝑥 feq , … , 𝑥 𝑓 , 𝑝 )Δ( 𝑥 feq , 𝑝 ; 0) 𝑣 feq 𝑣 𝑇 feq 𝑣 feq − 1 ⎤⎥⎥⎦ and 𝐺 𝐻 ( 𝑦 𝐻 ) = ⎡⎢⎢⎣ 𝑓 ( 𝑥 𝐻 , … , 𝑥 𝐻 , 𝑝 )Δ( 𝑥 𝐻 , 𝑝 ; i 𝜔 𝐻 ) 𝑣 𝐻 ̄𝑣 𝑇 ref 𝑣 𝐻 − 1 ⎤⎥⎥⎦ ,which are in ℝ 𝑛 +1 and ℝ 𝑛 +2 , respectively; here the base points ( 𝑥 feq , 𝑝 ) and ( 𝑥 𝐻 , 𝑝 ) areincluded as arguments of the characteristic matrix Δ (since the matrices 𝐴 𝑖 in (11) dependon them). The scale of the complex eigenvector 𝑣 𝐻 ∈ ℂ 𝑛 is fixed with the help of areference vector 𝑣 ref ∈ ℂ 𝑛 and one complex equation (i.e., two real equations).For the balancing pendulum model, the point 𝐻 is on a curve of Hopf bifurcations,which can be continued in two parameters. The resulting curve in the ( 𝑎, 𝑏 ) -plane isshown in green in fig. 1(b). The eigenvalue zero bifurcation at the point 𝑃 is not a fold, buta pitchfork bifurcation due to the equilibrium’s invariance under the reflection symmetry.Thus, the equation 𝐺 feq ( 𝑦 feq ) = 0 is singular. An experimental feature of DDE-BIFTOOLpermits the user to add constraints enforcing the symmetry of the equilibrium to make theextended nonlinear system regular (for example, 𝑥 feq , = 0 for the inverted pendulum);see also Section 3.10. For system (8) the pitchfork bifurcation is at 𝑎 = 1 for any 𝑏 and 𝜏 ,and this curve 𝑃 is shown in red in fig. 1(b).10long these two curves of codimension-one bifurcations one encounters degeneratepoints. At 𝜔 𝐻 = 0 the Hopf bifurcation meets the pitchfork bifurcation in a special pointwhere the linearization has a zero eigenvalue of algebraic multiplicity ; this point liesat 𝑏 = 𝜏 and 𝑎 = 1 in fig. 1(b). For small 𝜏 the Hopf bifurcation curve emerges to theright of the line of pitchfork bifurcations before it bends back to cross this line again at 𝑏 ≈ 4 . ; at this point there is a pitchfork-Hopf bifurcation, where the linearization has aneigenvalue and an imaginary eigenvalue pair ±i 𝜔 ℎ .The Hopf bifurcation curve and the pitchfork bifurcation line bound the region infig. 1(b) where the upright position 𝑥 = (0 , is linearly stable; this is indicated withthe label , which refers to the number of unstable eigenvalues. Bosschaert et al. (2020)implemented normal form analysis for generic equilibrium bifurcations of codimensionone or two. Their analysis produces normal form coefficients permitting us to branch offtoward secondary codimension-one branches from codimension-two bifurcation points.The expressions rely on genericity conditions that are violated at the pitchfork bifurcation(due to the reflection symmetry), such that the routines provided by Bosschaert et al.(2020) cannot be applied to the double-zero and pitchfork-Hopf interaction points wefound for the inverted pendulum DDE (8). However, the criticality of the Hopf bifurcationalong the branch of Hopf bifurcation can be determined by computing the Lyapunovcoefficient 𝓁 with the routines from (Bosschaert et al., 2020). The coefficient is negativeeverywhere, meaning that the Hopf bifurcation is supercritical. From this informationthe shown numbers of unstable eigenvalues of the upright equilibrium 𝑥 = (0 , in theother regions of fig. 1(b) can be determined. The region of stability shown in fig. 1(b) shrinks to the single point ( 𝑎 ∗ , 𝑏 ∗ ) = (1 , √ atthe critical delay 𝜏 ∗ = √ , at which the upright equilibirum 𝑥 = (0 , has a linearizationwith a zero eigenvalue of multiplicity . Hence, the critical delay 𝜏 ∗ is the largest possibledelay for which the delayed PD control can stabilize the upright position. In (Sieberand Krauskopf, 2004a,b) we performed an analysis of the bifurcation structure in theneighborhood of this singularity. Its results are shown in fig. 2 for parameters on anellipsoid around the singular point ( 𝑎 ∗ , 𝑏 ∗ , 𝜏 ∗ ) = (1 , √ , √ of codimension three. Theellipsoid is parametrized in the new coordinates ( 𝛼, 𝛽, 𝛾 ) , which are related to the originalcoordinates via 𝑎 = 𝑎 ∗ + 𝑟 𝛼 , 𝑏 = 𝑏 ∗ + 𝑟 𝜏 ∗ 𝛽 , 𝜏 = 𝑏 + 𝑟 𝜏 ∗ 𝛾 , (12) 𝛼 = sin( 𝜑𝜋 ∕2) , 𝛽 = cos( 𝜑𝜋 ∕2) cos(2 𝜋𝜓 ) , 𝛾 = cos( 𝜑𝜋 ∕2) sin(2 𝜋𝜓 ) , where the radius-type scaling parameter 𝑟 determines the size of the ellipsoid. The unitsphere in the rescaled ( 𝛼, 𝛽, 𝛾 ) -space is then parametrized by the two polar coordinateangles ( 𝜑, 𝜓 ) ∈ [0 ,
1] × [0 , , yielding the representation in fig. 2.The viewpoint of the sphere in ( 𝛼, 𝛽, 𝛾 ) -space is chosen in fig. 2 with a focus on theregion of stability of the upright pendulum, labeled (a) and bounded by the curves 𝑃 ofpitchfork bifurcation and 𝐻 of Hopf bifurcation, as well as on nearby regions with morecomplicated dynamics of the controlled inverted pendulum. Note from panel (a) of fig. 211 c)(a) (b)(a) (b)(c) ? P H P po ( ∗ ) ? P DT ∞ Figure 2:
Bifurcation diagram of the controlled inverted pendulum DDE on an ellipse around thetriple-zero point, with regions of equilibrium stabilization (a), small oscillations (b) and runawayoscillations (c) of the pendulum. The bounded chaotic dynamics in the area labelled ( ∗ ) is shownseparately in fig. 3(b). Curves of codimension-one bifurcations on the ellipse include those ofpitchfork bifurcation 𝑃 , Hopf bifurcation 𝐻 , pitchfork bifurcation of periodic orbits 𝑃 po , perioddoubling 𝑃 𝐷 , and connecting orbits 𝑇 ∞ . Here, 𝑟 = 1∕2 in the definition (12) of the ellipse. that successful stabilization of the upright position involves convergence of the position 𝜃 as well as of the velocity ̇𝛿 of the cart. As we discuss next, finding additional behaviorof the system beyond stabilization requires the continuation of periodic orbits and theirbifurcations. A periodic orbit 𝑥 ( 𝑡 ) with 𝑥 ( 𝑡 ) = 𝑥 ( 𝑡 − 𝑇 ) for some fixed period 𝑇 > and all times 𝑡 is given as the solution of a periodic DDE boundary-value problem (BVP). Hence,obtaining the nonlinear problem 𝐺 po ( 𝑦 po ) = 0 for the computation and continuation ofperiodic orbits requires a discretization of a periodic BVP. As is common, we rescale theperiod to the interval [0 , and use the variable 𝑠 = 𝑡 ∕ 𝑇 for the rescaled time (and recallthe convention that 𝜏 = 0 ) to obtainDDE: 𝑀 𝑥 po′ ( 𝑠 ) = 𝑇 𝑓 ( 𝑥 po ( 𝑠 − 𝜏 𝑇 ) [0 , , … , 𝑥 po ( 𝑠 − 𝜏 𝑑 𝑇 ) [0 , , 𝑝 ) , (13)BC: 𝑥 po (0) = 𝑥 po (1) , (14)PC: ∫ 𝑥 ′ref ( 𝑠 ) 𝑇 𝑥 po ( 𝑠 )d 𝑠 (15)for the unknowns 𝑥 po ( ⋅ ) ∈ 𝐶 ([0 , ℝ 𝑛 ) , 𝑇 ∈ (0 , ∞) and a single system parameter 𝑝 𝑖 . In (13) we use the notation 𝑥 po ( 𝑠 ) [0 , for the ‘wrapped’ evaluation of 𝑥 po ( 𝑠 ) , that12s, ( 𝑠 ) [0 , = ( 𝑠 − f loor( 𝑠 )) (and, hence, 𝑥 po ( 𝑠 ) [0 , = 𝑥 po ( 𝑠 − f loor( 𝑠 )) ), where f loor( 𝑠 ) is the largest integer less than or equal to 𝑠 . Thus, ( 𝑠 ) [0 , is always in [0 , and the boundary condition (BC) (14) simply expresses the periodicity of the solution. In theoriginal time, the periodic orbit is then 𝑡 ↦ 𝑥 po ( 𝑡 ∕ 𝑇 ) [0 , . A phase condition (PC) isrequired to select a unique and isolated solution of the overall BVP that represents theperiodic orbit. This is the case because, if 𝑥 po ( ⋅ ) [0 , is a solution of (13)–(14), thenso is 𝑥 po ( ⋅ + 𝑐 ) [0 , for any 𝑐 ∈ [0 , . We use here the integral phase condition (15),which fixes the free phase of the present solution 𝑥 po by minimizing the integral distance ∫ ( 𝑥 ref ( 𝑠 ) − 𝑥 po ( 𝑠 )) 𝑇 ( 𝑥 ref ( 𝑠 ) − 𝑥 po ( 𝑠 ))d 𝑠 to a periodic reference solution 𝑥 ref ; see Doedel(2007) for more details on phase conditions. We may write the infinite-dimensionalnonlinear problem (13)–(15) in the form (6) as 𝐺 po ( 𝑦 po ) = 0 , where 𝑦 po = ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝 𝑖 ) ∈ 𝐶 ([0 , ℝ 𝑛 ) × ℝ × ℝ . (16)Engelborghs et al. (2000a) and Engelborghs and Doedel (2002) constructed a fixed-degree piecewise polynomial collocation discretization for (13)–(15). The discretizationstores the approximate solution 𝑥 po ( ⋅ ) on a mesh 𝑠 e of 𝑁 𝑇 subintervals in polynomialpieces of degree 𝜅 , in the form of a vector 𝑥 po , d ∈ ℝ 𝑛 ( 𝜅𝑁 𝑇 +1) . The discretized residual atan arbitrary point 𝑠 ∈ [0 , has the form 𝐺 DE , d ( 𝑥 po , d , 𝑇 , 𝑝 ; 𝑠 ) = 𝑀 𝐸 (1) ( 𝑠 ) 𝑥 po , d − 𝑇 𝑓 ( 𝐸 (0) ( 𝑠 − 𝜏 𝑇 ) [0 , 𝑥 po , d , … , 𝐸 (0) ( 𝑠 − 𝜏 𝑑 𝑇 ) [0 , 𝑥 po , d , 𝑝 ) .The matrices 𝐸 ( 𝓁 ) ( 𝑠 ) are 𝑛 × ( 𝑛 ( 𝜅𝑁 𝑇 + 1)) interpolation (for 𝓁 = 0 ) and differentiation(for 𝓁 > ) matrices, such that 𝑥 ( 𝑠 ) = 𝐸 (0) ( 𝑠 ) 𝑥 po , d and 𝑥 ( 𝓁 ) ( 𝑠 ) = 𝐸 ( 𝓁 ) ( 𝑠 ) 𝑥 po , d for allpiecewise polynomials 𝑥 defined on the mesh 𝑠 e . The overall discretized periodic DDEBVP (13)–(15) has the form 𝐺 po , d ( 𝑦 po , d ) = ⎡⎢⎢⎢⎣ ( 𝐺 DE , d ( 𝑥 po , d , 𝑇 , 𝑝 ; 𝑠 c ,𝑗 ) ) 𝜅𝑁 𝑇 𝑗 =1 [ 𝐸 (0) (0) − 𝐸 (0) (1) ] 𝑥 po , d ∫ ( 𝐸 (1) ( 𝑠 ) 𝑥 ref , d ) 𝑇 𝐸 (0) ( 𝑠 ) 𝑥 po , d d 𝑠 ⎤⎥⎥⎥⎦ , (17)where ( 𝑠 c ,𝑗 ) 𝜅𝑁 𝑇 𝑗 =1 is a suitably chosen collocation mesh. The variable 𝑦 po , d = ( 𝑥 po , d , 𝑇 , 𝑝 𝑖 ) has dimension 𝑁 = 𝑛 ( 𝜅𝑁 𝑇 + 1) + 2 and 𝐺 po ,𝑑 ( 𝑦 po , d ) has 𝑛𝜅𝑁 𝑇 + 𝑛 + 1 components. Thematrix 𝐸 (0) ( ⋅ ) provides a natural embedding such that the function 𝑥 po , e ( ⋅ ) ∶ 𝑡 ↦ 𝑥 po , e ( 𝑡 ) = 𝐸 (0) ( 𝑡 ∕ 𝑇 ) [0 , 𝑥 po , d (18)is a piecewise polynomial approximation of the periodic orbit. The integral in the finalcomponent of 𝐺 po , d (the phase condition) can be computed exactly as 𝑥 𝑇 ref , d 𝑊 𝑐 𝑥 po , d withprecomputed quadrature weights 𝑊 c ; here the discretized reference solution 𝑥 ref , d isgenerally chosen as the discretized periodic orbit computed at the last continuation step.While the definitions for 𝐺 po , d look identical for constant and state-dependent delays, the 𝜏 𝑗 in the definition of 𝐺 DE , d are functions of 𝑠 and ( 𝑥 po , d , 𝑇 , 𝑝 ) if they depend on the state.13onvergence of the discretization for constant delays has only been proven recently byAndó and Breda (2020). Convergence proofs in earlier papers treated the variables 𝑇 and 𝜏 𝑗 as given constants, because considering (for example) the period 𝑇 as an unknowncreates analytical difficulties similar to those when considering state-dependent delays.The above description fits the current implementation in DDE-BIFTOOL, generalizingthe original construction by Engelborghs et al. (2000a) and Engelborghs and Doedel(2002) by separating the construction of the matrices 𝐸 ( 𝓁 ) ( ⋅ ) from the construction of thenonlinear problem. Since different discretization methods enter only through differentinterpolation matrices 𝐸 ( 𝓁 ) , different choices of discretization can be made dependingon the matrix 𝑀 . The mesh adaptation is based on the same error estimate and errorequidistribution as those in AUTO (Doedel et al., 1999).For the inverted pendulum model (8), the periodic orbits that emerge from the Hopfbifurcations can be computed as solutions of 𝐺 po , d ( 𝑦 po , d ) = 0 with 𝐺 po , d ( 𝑦 po , d ) given by(17), for example, with the angle 𝜙 as the continuation parameter 𝑝 𝑖 . These periodic orbitshave the spatio-temporal symmetry 𝑥 po ( 𝑡 + 𝑇 ∕2) = − 𝑥 po ( 𝑡 ) , and a typical time profileis shown in fig. 2(b). The symmetric periodic orbits lose their stability in a symmetrybreaking pitchfork bifurcation of periodic orbits, labelled 𝑃 po in fig. 2; see Sections 3.7and 3.8 for the linear stability analysis of periodic orbits. Moreover, a branch of symmetricperiodic orbits terminates when the period 𝑇 (which is part of the variable 𝑦 po and itsdiscretization 𝑦 po ,𝑑 along the solution branch) goes to infinity; at this point the periodicorbits approach a symmetric heteroclinic connection between the two non-symmetricsaddle equilibria, the locus of which is labeled 𝑇 ∞ in fig. 2. Together the curves 𝑃 po and 𝑇 bound the region of stable symmetric periodic orbits, which is labeled (b) in fig. 2. Thistype of periodic orbit can be interpreted as partial stabilization of the inverted pendulum:while the upright pendulum equilibrium is no longer stable, nevertheless, the periodicmotion of the pendulum as well as the velocity of the cart are bounded and initially (closeto the curve 𝐻 ) quite small. Linearizing around a solution ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝 ) of (13)–(15) yields a linear DDE with time-periodic coefficient matrices 𝑠 ↦ 𝐴 𝑗 ( 𝑠 ) and (if the delays are state-dependent) timeperiodic delays 𝑠 ↦ 𝜏 𝑗 ( 𝑠 ) of period . The eigenvalue problem involves a Floquetmultiplier 𝜇 ∈ ℂ and an eigenfunction 𝑥 ev ∶ [− 𝜏 max ∕ 𝑇 , ↦ ℂ 𝑛 , satisfying 𝑀 𝑥 ev′ ( 𝑠 ) − 𝑑 ∑ 𝑗 =0 𝐴 𝑗 ( 𝑠 ) 𝑥 ev ( 𝑠 − 𝜏 𝑗 ( 𝑠 ) 𝑇 ) if 𝑠 ∈ (0 , , (19) 𝜇𝑥 ev ( 𝑠 ) = 𝑥 ev ( 𝑠 + 1) if 𝑠 ∈ [− 𝜏 max ∕ 𝑇 , (20)(note the absence of wrapping). Here, for state-dependent delays with 𝜏 = 0 , 𝜏 𝑗 = 𝜏 𝑓𝑗 ( 𝑥 , … , 𝑥 𝑗 −1 , 𝑝 ) , 𝑋 = 𝐼 𝑛 , 𝑥 𝑗 = 𝑥 po ( 𝑠 − 𝜏 𝑗 ∕ 𝑇 ) [0 , , and 𝑥 ′ 𝑗 = 𝑥 po′ ( 𝑠 − 𝜏 𝑗 ∕ 𝑇 ) [0 , for14 = 0 , … , 𝑑 we have 𝐴 𝑗 ( 𝑠 ) = 𝑇 𝜕𝑓𝜕𝑥 𝑗 ( 𝑥 , … , 𝑥 𝑑 , 𝑝 ) 𝑋 𝑗 and (21) 𝑋 𝑗 = 𝐼 𝑛 − 𝑥 ′ 𝑗 𝑇 𝑗 −1 ∑ 𝓁 =0 𝜕𝜏 𝑓𝑗 𝜕𝑥 𝓁 ( 𝑥 , … , 𝑥 𝑗 −1 , 𝑝 ) 𝑋 𝓁 .In practice, the arguments of 𝜕𝑓 ∕ 𝜕𝑥 𝑗 have to be evaluated for the approximate periodicorbit defined by (18), that is, 𝑥 𝑗 = 𝑥 po , e ( 𝑠 − 𝜏 𝑗 ∕ 𝑇 ) [0 , and 𝑥 ′ 𝑗 = 𝑥 po , e ( 𝑠 − 𝜏 𝑗 ∕ 𝑇 ) [0 , .Szalai et al. (2006) and Sieber and Szalai (2011) showed that the eigenvalue problem(19)–(20) is equivalent to a finite-dimensional eigenvalue problem Δ( 𝜇 ) 𝑣 = 0 , wherethe dimension of the characteristic matrix Δ( 𝜇 ) may be larger than 𝑛 but is bounded forbounded 𝑇 and 𝜏 max . Yanchuk et al. (2019) constructed a characteristic matrix for thecase of a single delay 𝜏 (so, 𝑑 = 1 ), which can be large but keeps 𝑇 − 𝜏 bounded (whichis a common scenario for pulse-type periodic solutions in systems with a single largedelay).Extending the mesh 𝑠 e from the interval [0 , to the interval [− 𝜏 max ∕ 𝑇 , , a piecewisepolynomial 𝑥 of degree 𝜅 on [− 𝜏 max ∕ 𝑇 , is now stored in a vector 𝑥 ev , d of size 𝑛 ( 𝜅 ( 𝑁 𝑇 + 𝑁 𝜏 max ) + 1) . With the help of the interpolation and differentiation matrices 𝐸 ( 𝓁 ) ( ⋅ ) on theextended mesh 𝑠 e , we define the 𝑛 × ( 𝑛 ( 𝜅 ( 𝑁 𝑇 + 𝑁 𝜏 max ) + 1)) coefficient matrix for thediscretized residual of the linear DDE (19) at an arbitrary time 𝑠 ∈ (0 , 𝐴 DE , d ( 𝑠 ) = 𝑀 𝐸 (1) ( 𝑠 ) − 𝑑 ∑ 𝑗 =0 𝐴 𝑗 ( 𝑠 ) 𝐸 (0) ( 𝑠 − 𝜏 𝑗 ( 𝑠 ) 𝑇 ) .Then the discretized eigenvalue problem is a ( 𝑁 + 𝑁 ) -dimensional generalized matrixeigenvalue problem with ( 𝑁 , 𝑁 ) = ( 𝑛 ( 𝜅𝑁 𝜏 max + 1) , 𝑛𝜅𝑁 𝑇 ) of the form (Borgioli et al.,2020) 𝐴𝑥 ev , d = 𝜇𝐵𝑥 ev , d , with 𝐴 = [ ( 𝐴 DE , d ( 𝑠 c ,𝑗 )) 𝜅𝑁 𝑇 𝑗 =1 𝑁 × 𝑁 𝐼 𝑁 × 𝑁 ] and 𝐵 = [ 𝑁 ×( 𝑁 + 𝑁 ) 𝐼 𝑁 × 𝑁 𝑁 × 𝑁 ] .In the definition for 𝐴 we may use the same collocation mesh 𝑠 c as for the computationof the periodic orbit 𝑥 . The interpolation matrix 𝐸 (0) provides a natural embedding suchthat 𝑥 ev , e ∶ 𝑡 ↦ 𝐸 (0) ( 𝑡 ∕ 𝑇 ) 𝑥 ev , d is an approximate Floquet eigenfunction corresponding toa Floquet multiplier 𝜇 , where 𝑥 ev , d is the eigenvector of size ( 𝑁 + 𝑁 ) from the discretematrix eigenvalue problem. The dimension of the matrix eigenvalue problem can bereduced to an explicit 𝑁 × 𝑁 eigenvalue problem by using the first 𝑁 equations (where 𝜇 does not show up) to eliminate the components (( 𝑥 ev , d ) 𝑗 ) 𝑁 + 𝑁 𝑗 = 𝑁 +1 . This corresponds tosolving the discretized monodromy problem for the linear DDE (19), and works wheneverthe initial-value problem is well-posed. Yanchuk et al. (2019) observed that an adaptivemesh 𝑠 e that has been adapted for a good approximation of the periodic orbit 𝑥 po ( 𝑡 ) maygive poor approximations 𝑥 ev , e for the Floquet eigenfunctions 𝑥 ev ( 𝑡 ) even for multipliers15ith | 𝜇 | ≈ 1 . The eigenfunction 𝑥 ev ( 𝑡 ) may have rapid oscillations where 𝑥 po ( 𝑡 ) isapproximately constant. This is common when the delay 𝜏 and the period 𝑇 are relativelylarge and 𝑥 po is pulse-like (Yanchuk et al., 2019). For generic local bifurcations of periodic orbits DDE-BIFTOOL implements fully ex-tended defining systems (Govaerts, 2000) by appending the variational problem. Theextended system is formulated in the infinite-dimensional space in a manner that it hasagain the form (13)–(15) of a periodic DDE BVP with additional free parameters andintegral constraints. For the continuation of folds of periodic orbits, we consider theinfinite-dimensional nonlinear problem 𝐺 po ( 𝑦 po ) = 0 where 𝑦 po = ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝 𝑖 , 𝑝 𝑗 ) (thus,adding one free parameter), and append its variational problem, such that we have 𝐺 fpo ( 𝑦 fpo ) = 𝐺 fpo ( 𝑦 fpo , 𝑦 v ) = 𝐺 fpo ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝, 𝑥 v ( ⋅ ) , 𝑇 𝑣 )= ⎡⎢⎢⎣ 𝐺 po ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝 ) 𝜕 𝑥 po 𝐺 po ( 𝑦 po ) 𝑥 v + 𝜕 𝑇 𝐺 po ( 𝑦 po ) 𝑇 v ∫ 𝑥 v ( 𝑠 ) 𝑇 𝑥 v ( 𝑠 )d 𝑠 + 𝑇 v2 − 1 ⎤⎥⎥⎦ .Here, the variational variables 𝑦 v = ( 𝑥 v ( ⋅ ) , 𝑇 v ) ∈ 𝐶 ([0 , ℝ 𝑛 ) × ℝ have the sameformat as ( 𝑥 po ( ⋅ ) , 𝑇 ) . Thus, 𝐺 fpo ( 𝑦 fpo ) has the same format as (13)–(15): it consists ofa periodic BVP of dimension 𝑛 (for ( 𝑥 po ( ⋅ ) , 𝑥 v ( ⋅ )) ) with scalar integral conditionssuch that DDE-BIFTOOL uses the discretization 𝐺 po , d on this extended problem. Thediscretized problem has the overall dimension 𝑁 = 2 𝑛 ( 𝜅𝑁 𝑇 + 1) + 4 .For the pitchfork bifurcation of periodic orbits found in the bifurcation diagram of thecontrolled inverted pendulum this system is singular. However, an experimental feature ofDDE-BIFTOOL lets the user append conditions that enforce a spatio-temporal symmetryof 𝑥 po , such that the system becomes regular (for example, for the reflection symmetry ∫ ( 𝑥 po ) ( 𝑠 )d 𝑠 = 0 ); the bifurcation curve 𝑃 po in fig. 2 was computed in this way. Noticefrom fig. 2(c) that the bifurcating asymmetric periodic orbit in the corresponding regionrepresents oscillations around a value of 𝜃 that is not zero; in particular, the velocity ̇𝛿 ofthe cart decreases (or increases). Hence, this type of periodic orbit no longer satifies thecondition for generalized stability of the controlled inverted pendulum.For the continuation of period doubling and torus bifurcations we append the equationfor the critical eigenfunction in periodic form: if 𝑥 ev ( 𝑠 + 1) = 𝜇𝑥 ev ( 𝑠 ) and | 𝜇 | = 1 , wemay write 𝜇 = exp(i 𝜋𝜔 v ) , where 𝜋𝜔 v is the rotation number, and introduce 𝑥 v ( 𝑠 ) =exp(−i 𝜋𝜔 v 𝑠 ) 𝑥 ev ( 𝑠 ) . Then the defining system for torus and period doubling bifurcationsappends to (17) the equations 𝑥 v′ ( 𝑠 ) = −i 𝜋𝜔 v 𝑥 v ( 𝑠 ) + 𝑑 ∑ 𝑗 =0 𝐴 𝑗 ( 𝑠 ) 𝑥 v ( 𝑠 − 𝜏 𝑗 ( 𝑠 )∕ 𝑇 ) exp(−i 𝜋𝜔 v 𝜏 𝑗 ( 𝑠 )∕ 𝑇 ) , (22) ∫ 𝑥 v , ref ( 𝑠 ) 𝑇 𝑥 v ( 𝑠 )d 𝑠 − 1 , (23)16here the 𝐴 𝑗 are as defined by (21) in the section on linear stability of periodic orbits.Equation (23) has two real components (or one complex component). It fixes the lengthand phase of the complex Floquet vector 𝑥 v , depending on a reference function 𝑥 v , ref .Thus, the defining system 𝐺 tor ( 𝑦 tor ) = 0 for the torus bifurcation consists of 𝐺 po ( 𝑦 po ) = 0 defined by (17) together with (22)–(23) for the (extended) input vector variables 𝑦 tor =( 𝑦 po , 𝑥 v ( ⋅ ) , 𝜔 v ) = ( 𝑥 po ( ⋅ ) , 𝑇 , 𝑝 𝑖 , 𝑝 𝑗 , 𝑥 v ( ⋅ ) , 𝜔 v ) . Again, the solutions are both representedby discretization vectors 𝑥 po , d and 𝑥 v , d such that the overall dimension of 𝑦 tor is 𝑁 =3 𝑛 ( 𝜅𝑁 𝑇 + 1) + 5 . For the period doubling bifurcation the unknown rotation number 𝜋𝜔 v equals 𝜋 (such that 𝜔 v = 1 ).Apart from these generic codimension-one local bifurcation of periodic orbits, DDE-BIFTOOL is able to detect and continue connecting orbits between equilibria by means ofa defining system 𝐺 con ( 𝑦 con ) implemented by Samaey et al. (2002). As for the variable 𝑦 po for periodic orbits, the formulation is that of a DDE BVP and the variable 𝑦 con contains(discretized) orbit segments and an unknown integration time (which is, however, nolonger a period). Moreover, 𝑦 con also contains variables associated with the locations ofthe equilibria and their linearizations, because the setup uses projection boundary condi-tions to approximate the connecting orbit by an orbit segment with a finite integrationtime (Beyn, 1990). The DDE BVP 𝐺 con ( 𝑦 con ) = 0 can be used to detect connecting orbitsand continue them in a parameter plane when they are of codimension one; however, thepresent implementation is restricted to the constant-delay case. An alternative methodfor computing connecting orbits that arise as limits of periodic orbits, such as homoclinicorbits, is to continue a branch of periodic orbits until the period 𝑇 is very high, indicat-ing that the parameter is close to the locus of connecting orbits. The continuation ofperiodic orbits with a fixed, sufficiently high period in two system parameters is then anapproximation of the sought locus of connections.This latter approach was used for finding the curve 𝑇 ∞ for the controlled invertedpendulum in fig. 2 as the locus of symmetric periodic orbits of large period. Moreover,we found that for 𝜙 approaching the curve 𝑇 ∞ the branch of symmetric periodic orbits hasan infinite sequence of pitchfork and fold bifurcations when the heteroclinic connection isapproached; this is the case because the dominant eigenvalues of the linearizations at thenon-symmetric saddle equilibria satisfy the (saddle-quantity) condition for a complicated(or chaotic) Shilnikov bifurcation (Kuznetsov, 2013). The non-symmetric periodic orbitsbranching off at these pitchfork bifurcations encounter period doublings for larger 𝜙 andthen also approach a connecting orbit, namely a non-symmetric homoclinic connectionto a single non-symmetric equilibrium. The curve 𝑃 𝐷 in fig. 2, bounding the dark redregion, is the first of these period doubling bifurcations and it has been computed withthe defining system 𝐺 tor ( 𝑦 tor ) = 0 as defined above. The repeated period doublings of the non-symmetric oscillations in fig. 2 suggest that thereis a region where oscillations can be chaotic. We expect this to occur inside the wedgesformed by the first period doubling curve
𝑃 𝐷 , and a region with chaotic oscillations isindeed readily identified inside the largest wedge by numerical simulations. As is shownin row (a) of fig. 3, initially, these chaotic oscillations of (8) are non-symmetric and they17 a1)(b1) (a2)(b2) (a3)(b3)
Figure 3:
Chaotic runaway motion (a) and chaotic bounded motion (b) of the controlled invertedpendulum DDE; here, panels (a1) and (b1) show the trajectory and (a2) and (b2) the Poincare mapin projection onto the ( 𝑥 , 𝑥 ) -plane, while (a3) and (b3) show the time evolution of the velocity ̇𝛿 of the cart. From [Sieber and Krauskopf (2004b)] © 2004 Elsevier; reproduced with permission. do not correspond to successful balancing, because they feature run-away accelerationof the cart balancing the pendulum; see fig. 3(a3). However, further inside the chaoticregion, inside the subregion labelled ( ∗ ) in fig. 2, the two symmetrically related non-symmetric chaotic attractors collide near a homoclinic tangency of the symmetric saddleperiodic orbit. As a result, symmetric chaotic oscillations are possible, as is shown inrow (b) of fig. 3. For these symmetric chaotic oscillations of position 𝑥 and velocity 𝑥 of the pendulum, the cart performs a chaotic walk around its zero velocity; see fig. 3(b3).Therefore, this type of symmetric chaotic attractor represents a quite extreme form ofgeneralized stabilization of the inverted pendulum. The homoclinic tangency of thesaddle periodic orbit can be approximated by computing the unstable manifold of thesymmetric saddle orbit and checking if it returns to the symmetric saddle orbit in away that is tangent to the stable linear subspace of the symmetric saddle orbit. A briefdescription of computations of unstable manifolds of periodic orbits is given as part ofthe case study in Section 5; see Sieber and Krauskopf (2004b) for further details of theoverall dynamics of the DDE model (8) of the PD controlled inverted pendulum. We finish the description of the capabilities of DDE-BIFTOOL by mentioning brieflysome features that are still experimental.
Problems with discrete symmetry
Similar to the reflection symmetry of the pendulummodel (9), systems with discrete symmetries (for example where 𝑅𝑓 ( 𝑥 , … , 𝑥 𝑑 , 𝑝 ) = 𝑓 ( 𝑅𝑥 , … , 𝑅𝑥 𝑑 , 𝑝 ) 𝑅 of the identity matrix) have additional degeneracies when symmetrybreaking occurs. The user can add additional constraints enforcing the symmetry, andrequest that the defining system automatically appends dummy variables to create regularcontinuation problems for symmetry-breaking bifurcations. This feature has been usedto find the curve 𝑃 po of pitchfork bifurcation of periodic orbits of (8) on the ellipse inparameter space shown in fig. 2. Problems with rotational symmetry
A common feature of problems with delays inoptics is that they have rotational symmetry, in the simplest case with a single free rotation.That is, the right-hand side 𝑓 satisfies exp( 𝐴𝑡 ) 𝑓 ( 𝑥 , … , 𝑥 𝑑 , 𝑝 ) = 𝑓 (exp( 𝐴𝑡 ) 𝑥 , … , exp( 𝐴𝑡 ) 𝑥 𝑑 , 𝑝 ) for a fixed anti-symmetric matrix 𝐴 ∈ ℝ 𝑛 × 𝑛 (that is, 𝐴 𝑇 = − 𝐴 ), and arbitrary 𝑡 ∈ ℝ and 𝑥 , … , 𝑥 𝑑 ∈ ℝ 𝑛 . In this case one is interested in tracking (1) any rotating wave, which isa periodic solution of the form 𝑥 ( 𝑡 ) = exp( 𝐴𝜔𝑡 ) 𝑥 eq , as a relative equilibrium, that is, as 𝑥 eq ; and (2) any modulated rotating wave, which is a quasi-periodic solution of the form 𝑥 ( 𝑡 ) = exp( 𝐴𝜔𝑡 ) 𝑥 po ( 𝑡 ) where 𝑥 po has period 𝑇 , as a relative periodic orbit, that is, as 𝑥 po .Typical examples are the DDE models for semiconductor lasers with optical feedbackfrom mirrors or with delayed optical coupling, which are invariant under rotation of thecomplex electric field (Pieroux et al., 2001; Haegeman et al., 2002; Krauskopf et al., 2000;Krauskopf, 2005). In fact, the wish to perform the bifurcation analysis of this type of lasersystem provided considerable motivation for the early development of DDE-BIFTOOL.In order to select unique solutions in the presence of this rotational symmetry, the phase ofthe electric field was pushed into an additional parameter, which was then set to a specific,fixed value to implement a phase condition. Generalizing this approach, DDE-BIFTOOLprovides a wrapper around the defining systems in Sections 3.2, 3.4, 3.6 and 3.8 thatintroduces the mean rotation frequency 𝜔 as an additional free parameter and adds anadditional phase condition. Interface with COCO
DDE-BIFTOOL can only continue curves of invariant opbjectsand their bifurcations, that is, it computes one-dimensional solution manifolds of therespective defining system. The package COCO, on the other hand, has implementedalgorithms for multi-parameter continuation by means of growing atlases of solutionmanifolds of arbitrary dimension, which is based on the original algorithm by Henderson(2002); see Dankowicz and Schilder (2013) for more details. An interface is availablethat feeds the defining systems implemented in DDE-BIFTOOL into COCO’s moregeneral continuation algorithm, and provides constructors to start such a multi-parametercontinuation.
The capabilities we just reviewed, which allow DDE-BIFTOOL to perform the varioustasks required for bifurcation analysis, can be applied to standard DDEs with constant orstate-dependent delays. The latter requires the user to specify the state-dependence in the19orm (2); how this is done in practice is demonstrated in Section 4 for a conceptual DDEmodel of the ENSO system and in Section 5 for a scalar DDE with two state-dependentdelays.Moreover, DDE-BIFTOOL can be used for the study of DDEs beyond the standard form,whose formulation requires that the matix 𝑀 in (1) is not the identity matrix. Permittingthe matrix 𝑀 to be singular drastically expands the class of DDE one can consider. Inparticular, the general form (1) used by DDE-BIFTOOL includes the following types ofDDEs. Neutral DDEs feature delayed derivatives, and this case can be formulated in the frame-work of (1) by a suitable choice of 𝑀 . As an example, including an acceleration depen-dence into the feedback term 𝐹 in (7) changes the governing equations for the deviation 𝜃 ( 𝑡 ) of the pendulum angle from the upright position to 𝜃 ′′ ( 𝑡 ) = sin 𝜃 ( 𝑡 ) − cos 𝜃 ( 𝑡 )[ 𝑎𝜃 ( 𝑡 − 𝜏 ) + 𝑏𝜃 ′ ( 𝑡 − 𝜏 ) + 𝑐𝜃 ′′ ( 𝑡 − 𝜏 )] ,where 𝑐 is an additional control gain (Sieber and Krauskopf, 2005; Insperger et al., 2013).This can be formulated within (1) for 𝑥 ( 𝑡 ) = ( 𝜃 ( 𝑡 ) , 𝜃 ′ ( 𝑡 ) , 𝜃 ′′ ( 𝑡 )) ∈ ℝ and 𝑝 = ( 𝑎, 𝑏, 𝑐, 𝜏 ) ∈ ℝ by setting 𝑀 = diag(1 , , and 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = ⎡⎢⎢⎣ 𝑥 𝑥 sin 𝑥 − cos 𝑥 [ 𝑝 𝑥 + 𝑝 𝑥 + 𝑝 𝑥 ] − 𝑥 ⎤⎥⎥⎦ . In a number of applications one encounters algebraic constraints on the variables of aDDE, leading to a system of DAEs with delays. For example, one may define a state-dependent delay implicitly, as is done in the position control problem with echo locationmeasurements discussed by Walther (2002) 𝑦 ′ ( 𝑡 ) = 𝑘 [ 𝑦 ref − 𝑐 𝑠 ( 𝑡 − 𝜏 ) ] , where 𝑐𝑠 ( 𝑡 ) = 𝑦 ( 𝑡 ) + 𝑦 ( 𝑡 − 𝑠 ( 𝑡 )) .Here 𝑦 ( 𝑡 ) is the position to be kept at target 𝑦 ref by feedback control. The other state, 𝑠 ( 𝑡 ) ,is the travel time of an echo location signal sent from position 𝑦 ( 𝑡 − 𝑠 ( 𝑡 )) with speed 𝑐 to 𝑦 ref and reflected back to 𝑦 ( 𝑡 ) to estimate the current position offset as 𝑐𝑠 ( 𝑡 )∕2 . Theconstant delay 𝜏 is a reaction delay in the application of the feedback control with targetposition 𝑦 ref (similar to the delay in the pendulum feedback in (7), where the referenceposition is the upright angle ). Here 𝑥 ( 𝑡 ) = ( 𝑦 ( 𝑡 ) , 𝑠 ( 𝑡 )) ∈ ℝ and 𝑝 = ( 𝑦 ref , 𝑘, 𝑐, 𝜏 ) ∈ 𝑅 .This can be formulated by setting 𝑀 = diag(1 , , 𝑑 = 2 and 𝜏 𝑓 ( 𝑥 , 𝑝 ) = 𝑝 , 𝜏 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = 𝑥 , and 𝑓 ( 𝑥 , 𝑥 , 𝑥 , 𝑝 ) = ⎡⎢⎢⎢⎣ 𝑝 ( 𝑝 − 𝑝 𝑥 ) 𝑝 𝑥 − 𝑥 − 𝑥 ⎤⎥⎥⎥⎦ .20 .11.3 Forward-backward/mixed-type equations Even when all delays are positive in (1), the possibility of adding algebraic equationspermits one to introduce both negative and positive delays. These types of equations donot describe well-posed initial-value problems but may occur when modeling travelingwaves or periodic wave trains on a space-discrete lattice (Abell et al., 2005). For example,a wave in a discrete linear diffusion equation traveling with speed 𝜏 satisfies 𝑢 ′ ( 𝑡 ) =Δ[ 𝑢 ( 𝑡 + 𝜏 ) + 𝑢 ( 𝑡 − 𝜏 ) − 2 𝑢 ( 𝑡 )] . This could be formulated by setting 𝑀 = diag(1 , , 𝑝 = (Δ , 𝜏 ) ∈ ℝ , and 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = [ 𝑝 [ 𝑥 + 𝑥 − 2 𝑥 ] 𝑥 − 𝑥 ] ,such that 𝑥 ( 𝑡 ) ∈ ℝ and 𝑢 ( 𝑡 ) = 𝑥 ( 𝑡 ) . While the types of problems above can be formulated for input in DDE-BIFTOOL byusing a singular matrix 𝑀 in (1), we stress that any subsequent computations of invariantobjects and their bifurcation must be considered as being of an experimental nature.Namely, the accuracy of the results obtained by the different numerical computationswe described for standard DDEs is not always guaranteed. See Barton et al. (2006) foran analysis of convergence properties for neutral DDEs and note that forward-backwardproblems and delayed DAEs with index higher than are yet untested. Feedback loops are crucial ingredients in the dynamics of climate systems, where theyarise due to the interactions between various subsystems, including distinct bodies ofwater, the atmosphere, land and ice masses; see, for example, (Bar-Eli and Field, 1998;Dijkstra, 2008, 2013; Kaper and Engler, 2013; Keane et al., 2017; Simonnet et al., 2009)as entry points to the literature. Such feedback loops are subject to inherent time delays,mainly as a result of the time it takes to transport mass or energy across the globeand/or throughout the atmosphere, or due to delayed reactions of subsystems to changingconditions. Whenever the time delays of feedback loops in climate systems are largecompared to the forcing time scales under consideration, explicit modeling of the delaymakes sense in conceptual models. As a specific example of immediate human and widermathematical interest, we consider the El Niño phenomenon — a large increase of thesea surface temperature in the eastern equatorial Pacific Ocean that occurs about every3–7 years. This oceanic phenomenon is associated with an atmospheric component, theSouthern Oscillation, and they are jointly known as El Niño Southern Oscillation (ENSO)variability (Dijkstra, 2008; Graham and White, 1988; Kaper and Engler, 2013; Tzipermanet al., 1998; Zaliapin and Ghil, 2010). Large peaks in the sea surface temperature of theeastern Pacific Ocean near the coast off Peru represent El Niño events, the warm phaseof ENSO, while large drops represent the cool phase known as La Niña.21 tmosphere h + + R o s s b y w a v e s K e l v i n w a v e s Equator thermoclinedepth
Figure 4:
Schematic of the negative feedback loop of the thermocline ℎ at the eastern PacificOcean due to energy transport via Rossby and Kelvin wave from the central ocean-atmosphereinteraction zone. El Niño events have major consequences world-wide, yet they remain notoriouslyhard to predict even with sophisticated global climate models (Barnston et al., 2012).An important aspect of ENSO is that El Niño events tend to occur at the same timeof year, always around Christmas. This suggests locking to the seasonal cycle (with aperiod of 1 year), which represents the characteristic forcing time scale of the ENSOsystem. Feedback mechanisms in ENSO arise naturally from ocean-atmosphere couplingprocesses in the eastern and central equatorial Pacific Ocean, and they have delay timesof many months due to the time it takes waves to propagate across the Pacific Ocean.In light of the overall complexity of climate systems, conceptual models have much tooffer in terms of elucidating underlying mechanisms behind observed dynamics. Concep-tual DDE models for ENSO (and some other climate phenomena) have been developedby Bar-Eli and Field (1998); Dijkstra (2013); Falkena et al. (2019); Ghil et al. (2008);Kaper and Engler (2013); Tziperman et al. (1994, 1998); Zaliapin and Ghil (2010) toprovide insights into the interplay between delayed feedback loops and different typesof external forcing. Such DDE ENSO models constitute a significant model reduction,compared to the full description of atmospheric and oceanic dynamics and interaction,including their velocity and temperature fields. Since the feedback loops and their delaytimes are explicit parts of the DDE model, their roles for observed system behaviour canbe investigated readily. Generally, the delays that arise in such models are estimated,from quantities such as average wave speeds and distances, and taken to be constant.
We consider here an ENSO DDE model that follows the prominent delayed actionoscillator (DAO) paradigm that was first introduced by Suarez and Schopf (1988). Thereexist a number of models based on the DAO paradigm; for example, see (Suarez andSchopf, 1988; Battisti and Hirst, 1989; Tziperman et al., 1994, 1998). The DDE model22ntroduced by Ghil et al. (2008) is one of the simplest in that it focuses on the interactionbetween the negative delayed feedback and additive seasonal forcing. Its ingredientsare illustrated in fig. 4. The main quantity of interest for the DAO is the depth of thethermocline, which is the thin and distinct layer in the ocean that separates deeper coldwaters from shallower warm surface water. The thermocline is deeper in the Westand shallower in the East, and this is represented by the tilted bottom plane of fig. 4.The variable ℎ ( 𝑡 ) denotes the deviation of the thermocline depth from the long-termthermocline mean in the eastern equatorial Pacific, off the coast of Peru. A positive valueof ℎ corresponds to a larger layer of warm water and, hence, an increased sea-surfacetemperature (SST) in the eastern Pacific, while negative ℎ means a decreased SST. Inother words, the variable ℎ can be seen as a proxy for SST. We are concerned here withthe main negative feedback loop: off the equator, a negative anomalous thermoclinedepth signal is carried to the western boundary of the ocean via so-called Rossby waves,which are reflected as Kelvin waves. The oceanic waves of the negative feedback are thegreen arrows in fig. 4, and they carry the shallow thermocline perturbation back to theeastern boundary of the ocean, which is a process that takes on the order of 8 months.There is also a positive feedback of ℎ with a shorter delay of only about a month. It isrepresented by the black arrow in fig. 4 and arises from the fact that a warm SST anomalyslows down the easterly trade winds, leading to westerly wind anomalies that deepen thethermocline; see, for example, (Dijkstra, 2013; Keane et al., 2017, 2019) for more details. As was done by Ghil et al. (2008), we now consider only the above-mentioned negativefeedback loop and the seasonal cycle, with the goal of demonstrating that their interplayis sufficient to produce rich dynamical behaviour that is relevant to ENSO. The effects ofincluding the positive feedback loop into this model are studied in detail in Keane et al.(2016). The model from (Ghil et al., 2008), which we refer to as the GZT model fromnow on, takes the form ℎ ′ ( 𝑡 ) = − 𝑏 tanh[ 𝜅ℎ ( 𝑡 − 𝜏 )] + 𝑐 cos(2 𝜋𝑡 ) . (24)Here, 𝜏 is the delay time of the negative feedback loop with amplification factor 𝑏 , whichis further characterized by the coupling parameter 𝜅 ; note that 𝜅 is the slope at 0 of the tanh -function and see (Münnich et al., 1991) for a justification for this simple type ofocean-atmosphere coupling. Throughout, we fix the parameters 𝑏 and 𝜅 to the values 𝑏 = 1 and 𝜅 = 11 that were used and justified in previous investigations of the ENSOphenomenon (Ghil et al., 2008; Zaliapin and Ghil, 2010). In (24) the periodic forcing ofstrength 𝑐 enters as an additive term. Alternatively, one may shift the period-1 responseto the origin such that the seasonal forcing is parametric, as was considered in a simpleDDE ENSO model by Tziperman et al. (1998); Krauskopf and Sieber (2014). We remarkthat simple conceptual ENSO DDEs such as the GZT model are of wider interest becausethey are rather prototypical: DDE models of much the same structure can also be foundin control theory and machining; see, for example, (Just et al., 2007; Milton et al., 2009;Purewal et al., 2014; Stépán, 1989). 23 key feature of system (24) is the periodic forcing term with its explicit dependenceon time 𝑡 ; hence, this DDE is non-autonomous. Since DDE-BIFTOOL is designed forautonomous DDEs, we transform (24) into autonomous form by introducing an artificialstable oscillation that generates the periodic forcing. For any periodically forced DDE(or ODE) this can be achieved with the Hopf normal form for a stable periodic orbit ofradius 1; it can be written in complex form as ̇𝑧 ( 𝑡 ) = (1 + 𝜔𝑖 ) 𝑧 ( 𝑡 ) − 𝑧 ( 𝑡 ) | 𝑧 ( 𝑡 ) | , (25)and we set 𝜔 = 2 𝜋 to have the required forcing period of 1 (year). This two-dimensionalsystem then drives (24) in its rewritten form ̇ℎ ( 𝑡 ) = − 𝑏 tanh [ 𝜅ℎ ( 𝑡 − 𝜏 )] + 𝑐 Re( 𝑧 ( 𝑡 )) . (26)The equivalent autonomous system (26) with (25) is readily implemented in DDE-BIFTOOL with physical state 𝑥 ( 𝑡 ) = ( ℎ ( 𝑡 ) , Re( 𝑧 ( 𝑡 )) , Im( 𝑧 ( 𝑡 ))) ∈ ℝ , parameter vector 𝑝 = ( 𝑝 , 𝑝 , 𝑝 , 𝑝 ) = ( 𝑏, 𝑐, 𝜅, 𝜏 ) ∈ ℝ (such that 𝜏 = 𝑝 ), and right-hand side 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = ⎡⎢⎢⎣ − 𝑝 tanh[ 𝑝 𝑥 ] + 𝑝 𝑥 𝑥 − 2 𝜋𝑥 − 𝑥 (( 𝑥 ) + ( 𝑥 ) ) 𝑥 + 2 𝜋𝑥 − 𝑥 (( 𝑥 ) + ( 𝑥 ) ) ⎤⎥⎥⎦ . (27)fig. 5 presents five stable solutions of (24) as obtained by numerical integration withthe Euler method from initial conditions ℎ ≡ and/or ℎ ≡ after transients have settleddown. Shown are the respective time series of ℎ in panels (a1)–(e1), which are intuitivein the context of ENSO system since the variable ℎ is a proxy for the SST: maxima andminima of ℎ represent El Niño and La Niña events, respectively. Panels (a2)–(e2) offig. 5 are projections of the corresponding attractors onto the ( ℎ ( 𝑡 − 𝜏 ) , ℎ ( 𝑡 )) -plane.For zero seasonal forcing 𝑐 = 0 there is an attracting periodic solution; see row (a) offig. 5. Its zigzag-like shape and period of 𝑇 = 4 𝜏 years is due to the fact that the slope 𝜅 is quite large at 𝜅 = 11 , so that the tanh -function is rather close to a discontinuousswitching function. In the context of ENSO, this stable periodic orbit corresponds to anEl Niño event exaclty every 4.8 years as driven by the delay time of 𝜏 = 1 . years. Onthe other hand, when the periodic forcing is large compared to the negative feedback onefinds a periodic solution that is quite close to sinusoidal with a period 𝑇 = 1 year; anexample is shown in row (b) of fig. 5. The observed dynamics is clearly dominated bythe seasonal forcing, meaning that the SST varies exactly with the seasonal cycle.The interesting case is that of an interplay between negative feedback and the seasonalcycle when 𝑏 and 𝑐 are of the same order. In this regime one may find dynamics oninvariant tori, which may be locked or quasiperiodic. Row (c) of fig. 5 shows a stablelocked periodic solution; in fact, it coexists with the seasonally driven periodic solutionin row (b). Another example of multi-stability are the stable solutions shown in rows (d)and (e). The projection onto the ( ℎ ( 𝑡 − 𝜏 ) , ℎ ( 𝑡 )) -plane in panel (d2) clearly shows thatthere is an attracing torus with dynamics that is quasiperiodic (or periodic with a very highperiod); the attractor in row (e), on the other hand, is clearly periodic and likely a lockedsolution on a different torus. Notice the difference in amplitude between the respective24 ✺ ✶✵✲✶✵✶ ✲✶ ✵ ✶ (a1) (a2) h ( t ) (cid:0) ✁ ✂(cid:0)✄✂(cid:0)✂ ✄✂ (cid:0) ✂ (b1) (b2) h ( t ) ☎ ✆ ✝☎✞✟☎✟ ✞✟ ☎ ✟ (c1) (c2) h ( t ) ✠ ✡ ☛✠☞✠✌✡✠✠✌✡ ☞✠✌✡ ✠ ✠✌✡ (d1) (d2) h ( t ) ✍ ✎ ✏✍✑✏✍✏ ✑✏ ✍ ✏ (e1) (e2) h ( t ) t h ( t − τ ) Figure 5:
Stable solutions of (24), shown as time series in panels (a1)–(e1) and as projectionsonto the ( ℎ ( 𝑡 − 𝜏 ) , ℎ ( 𝑡 )) -plane in panels (a2)–(e2); throughout 𝑏 = 1 , 𝜅 = 11 and 𝜏 = 1 . , 𝑐 = 0 for (a), 𝜏 = 1 . , 𝑐 = 3 for (b) and (c), and 𝜏 = 0 . , 𝑐 = 3 for (d) and (e). From [Keane et al.(2015)] © 2015 Society for Industrial and Applied Mathematics; reproduced with permission. coexisting stable solutions in rows (b) and (c) and in rows (d) and (e), respectively. Aninterpretation of the locked periodic solutions in row (c) and (e) would be a build-up ofthe maxima of the SST from year to year until a global maximum, interpreted as an ElNiño event, is reached and the SST decreases until a global minimum, interpreted as LaNiña event, is reached and the process repeats.fig. 6 shows two maximum maps in the ( 𝑐, 𝜏 ) -plane of (24), together with curves ofsaddle-node bifurcations of periodic orbits SL, period-doubling bifurcations PD, andtorus bifurcations T. Two-parameter maximum maps, which plot for each point of a gridin parameter space the maximum of a sufficiently long time series after transients havesettled down, have been considered as a convenient way of obtaining an overview of theoverall dynamics (Ghil et al., 2008; Keane et al., 2015). We show two maximum mapsin panels (a) and (b) of fig. 6 in a greyscale where, for each row of fixed delay 𝜏 , theparameter 𝑐 is swept up or down in small steps as is indicated by the arrows; here the25 ax( h ( t )) τ c c SLPD T SLPD T
Figure 6:
Maximum maps and bifurcation set of (24) in the ( 𝑐, 𝜏 ) -plane with curves of saddle-node bifurcations of periodic orbits (SL), period-doubling (PD) and torus bifurcations (T), andlabelled lower resonance tongues. From [Keane et al. (2015)] © 2015 Society for Industrial andApplied Mathematics; reproduced with permission. respective previous solution is used as the initial history for the next value of 𝑐 . In thisway, hysteresis loops in 𝑐 are detected and regions of multistability are identied as regionsin the ( 𝑐, 𝜏 ) -plane where the two maximum maps do not agree.The bifurcation curves in fig. 6 explain features of the two maximum maps. In particular,the curves SL of saddle-node bifurcations of periodic orbits delineate the elongated shapes.In fact, they bound resonance tongues that emerge from the line 𝑐 = 0 of zero forcing andfrom the curve T of torus bifurcations, namely at points of 𝑝 ∶ 𝑞 resonance, some of whichare labelled. In the regions bounded by respective curves SL one finds stable frequencylocked solutions of the fixed frequency ratio 𝑝 ∶ 𝑞 ; compare with fig. 5(c) and (e). Thecurve T lies near the (roughly diagonal) boundary along which one finds sudden jumpsof the maxima. Notice that this boundary is different for increasing 𝑐 in panel (a) versusdecreasing 𝑐 in panel (b) of fig. 6, showing that the curve T is associated with regionsof multistability. As is discussed by Keane and Krauskopf (2018), this involves foldingresonance tongues and the break-up of invariant tori in what are known as Chencinerbubbles. Overall, fig. 6 confirms that for sufficiently large 𝑐 solutions are dominatedby the seasonal forcing, while there is an interplay between the forcing and the delayedfeedback for lower values of 𝑐 , specifically, to the left of the torus bifurcation curve T,where one finds dynamics on invariant tori.We focus here on the existence of chaotic dynamics caused by this interplay, becauseit has been suggested that irregular locked motion of (24) captures important aspects ofENSO (Ghil et al., 2008; Zaliapin and Ghil, 2010). Regions where such dynamics mayoccur are those in fig. 6 that are bounded by curves PD of period-doubling bifurcations,which are found inside some of the shown resonance tongues. This is a known featurethat occurs when resonance tongues overlap and the corresponding tori lose their normalhyperbolicity and break up; see, for example, (Broer et al., 1998; Kuznetsov, 2013). Toidentify where chaotic dynamics can be found, fig. 7 shows the bifurcation curves inthe ( 𝑐, 𝜏 ) -plane overlaid on a map of the maximal Lyapunov exponent, as computed foran upsweep of 𝑐 with the algorithm for DDEs from (Farmer, 1982). It shows that apositive maximal Lyapunov exponent indicating chaotic dynamics is generally associated26 < 000.10.20.30.4 ¯ τ m a x i m a l L y a pun o v e x p o n e n t PPi
QQk (cid:27)
SLPD T
Figure 7:
Bifurcation set and maximal Lyapunov exponent of solutions in the ( 𝑐, 𝜏 ) -plane of (24)[corresponding to that of (32) with 𝜂 𝑐 = 0 and 𝜂 𝑒 = 0 ]. From [Keane et al. (2019)] © 2019 TheRoyal Society; reproduced with permission. with period-doubling cascades and can be found only in small regions of the parameterplane. Note that some of these regions of positive maximal Lyapunov exponents can befound where no curves of period-doublings are shown; indeed, there exist infinitely manyhigher-order resonance tongues in between those we have shown, and they are expectedto overlap. As demonstrated for a related DDE model by Keane et al. (2016), regions ofchaotic dynamics may also be entered via intermittent transitions that are characterisedby the sudden appearance of chaos at a saddle-node bifurcation (Pomeau and Manneville,1980). The conclusion to be drawn from fig. 7 is that chaotic dynamics of (24) can befound, but only for quite specific and small ranges of 𝑐 and 𝜏 . As we will show next, statedependence of the feedback loop changes this picture considerably. It is important to recognise that taking a constant value for any delay in a DDE model isa modelling assumption that must be justified. The assumption of delays being constantis well justified in certain applications, such as machining (Insperger and Stépán, 2000)and laser dynamics (Kane and Shore, 2005). On the other hand, delays in many applica-tions, and certainly in climate modelling, are definitely not constant. While the delays inconceptual DDE climate models have generally been taken to be constant (Keane et al.,2017), there are many reasons to suspect that this is not actually be the case. Generally,the delay times will depend on the state of the system itself, which leads to DDE modelswith state-dependent delays. The main questions that need to be addressed from a moregeneral modelling perspective are: 27 a s t e r l i e s ocean surfaceinteraction zone upwelling t h e r m o c l i n e m e a n t h e r m o c l i n e w a v e p r o p a g a t i o n ocean adjustment[ τ c ] [ τ e ] [ τ w ] [ ] Figure 8:
Ocean adjustment and upwelling as sources of state-dependence in the ocean-atmosphereinteraction in the central and eastern equatorial Pacific Ocean. From [Keane et al. (2019)] © 2019The Royal Society; reproduced with permission. (1) When does state dependence arise from physical processes and what mathematicalforms does it take?(2) Does state dependence of delays have a significant effect for the observed dynamicsof the respective DDE model?Specifically for the GZT model (24), a non-constant delay in the negative feedback looparises from the physics of the coupling of the ocean surface with the thermocline below.fig. 8 illustrates the heuristic argument for considering two terms with state dependencein the overall negative delay loop of ENSO, which are not described by and go beyondthe original DAO mechanism; more details can be found in Keane et al. (2019). Thehorizontal direction represents longitude along the equator between the basin boundariesof the Pacific Ocean and the vertical direction represents depth below the ocean surfacewith the atmosphere above. The thermocline is sketched as a deviation from its mean,which is about 50 metres deep in the East and 150 metres deep in the central equatorialPacific Ocean. The black arrows represent the four components of the negative feedbackloop. A positive perturbation in the thermocline depth ℎ ( 𝑡 ) in the eastern equatorialPacific increases the SST after an upwelling process with associated delay time 𝜏 𝑒 . Theeasterly winds forming the atmospheric component that transports such a perturbation tothe interaction zone in the central Pacific Ocean are considered fast and are modelled asinstanteaneous (as is the case in all DAO models). The interaction zone is reasonablylocalised, and it is simplified to a point in mathematical derivations of DAO models (Caneet al., 1990; Jin, 1997). There is then a delay 𝜏 𝑐 due to the coupling process known asocean adjustment of the SST influencing the thermocline, which depends on the currentthermocline depth ℎ ( 𝑡 ) . As in the GZT model without state dependence, Rossby wavesthen carry the signal to the western basin boundary and are reflected as Kelvin waves,which carry the signal back to the East with an associated delay time 𝜏 𝑤 , which weassume here is constant. The total delay time associated with the negative feedback loopis therefore 𝜏 = 𝜏 𝑒 + 𝜏 𝑐 + 𝜏 𝑤 . (28)28ere 𝜏 𝑒 and 𝜏 𝑐 are state dependent, that is, depend on the thermocline depth ℎ . Todetermine their functional form we summarize briefly the modelling exercise in Keaneet al. (2019), where more details can be found. It is convenient to define the constant partof the delay 𝜏 (with respect to the mean thermocline depth) as ̄𝜏 = ̄𝜏 𝑒 + ̄𝜏 𝑐 + 𝜏 𝑤 . (29)Here ̄𝜏 𝑒 is the constant time it takes the signal to travel from the mean thermoclinedepth to the surface, and ̄𝜏 𝑐 is the constant time of the ocean adjustment at the centralPacific associated with the mean thermocline depth. From a correlation analysis ofobservational SST and thermocline depth data (Zelle et al., 2004) one concludes that thetwo constant delay times ̄𝜏 𝑒 and ̄𝜏 𝑐 for the long-term average of the thermocline are 2weeks and 4 months, respectively, which gives the values ̄𝜏 𝑒 = 2∕52 and ̄𝜏 𝑐 = 4∕12 (inyears) that we use from now on. Moreover, based on oceanic wave speeds calculatedfrom TOPEX/POSEIDON satellite data in (Boulanger and Menkes, 1995; Chelton andSchlax, 1996), realistic values of 𝜏 𝑤 lie between 5.2 and 7.2 months, that is, in the range [0 . , . when scaled to years. Hence, one obtains the estimated range [0 . , . forthe constant part ̄𝜏 of the overall delay.The upwelling delay can be modelled by 𝜏 𝑒 = ̄𝜏 𝑒 + 𝜂 𝑒 ℎ ( 𝑡 − ̄𝜏 ) , (30)where 𝜂 𝑒 is the inverse of the upwelling speed. Note that the state-dependent term isitself subject to a delay because the thermocline depth signal that ultimately returns tothe eastern equatorial Pacific at time 𝑡 began its journey at the thermocline one feedbackcycle ago; in (30) this implicitly defined state dependence is resolved by consideringthe first-order approximation given by the constant part ̄𝜏 . Maximum deviations in thethermocline depth in the eastern equatorial Pacific Ocean are about 50 metres (Harrisonand Vecchi, 2001) and it follows, with time measures in years, that the nominal value ofthe inverse upwelling speed is 𝜂 𝑒 ≈ 2∕52 ≈ 0 . .The dependence of the delay time 𝜏 𝑐 due to mass transport between ocean surface andthe thermocline can be modelled by 𝜏 𝑐 = ̄𝜏 𝑐 + 𝜂 𝑐 ℎ ( 𝑡 ) , (31)where 𝜂 𝑐 is the ocean-adjustment speed. Since the maximum deviations in thermoclinedepth in the central equatorial Pacific Ocean of 150 metres corresponds to about onethird of its mean depth, we obtain similarly the nominal value 𝜂 𝑐 ≈ (4∕3)∕12 ≈ 0 . (inunits of years per meter). The resulting state-dependent GZT ENSO DDE model we consider in what follows isgiven by (24) with the overall state-dependent delay 𝜏 ( ℎ ) = ̄𝜏 + 𝜂 𝑒 ℎ ( 𝑡 − ̄𝜏 ) + 𝜂 𝑐 ℎ ( 𝑡 ) , (32)29hich has the additional parameters 𝜂 𝑒 and 𝜂 𝑐 that allow us to ‘switch on’ the two types ofstate dependence. Clearly, for 𝜂 𝑒 = 𝜂 𝑐 = 0 this model reduces to the constant-delay GZTDDE. When implementing the state-dependent delay in DDE-BIFTOOL, the expressionfor the right-hand side is very similar to that given in (27), but the parameter vector ischanged to 𝑝 = ( 𝑝 , … , 𝑝 ) = ( 𝑏, 𝑐, 𝜅, ̄𝜏, 𝜂 𝑐 , 𝜂 𝑒 ) ∈ ℝ , the number of delays is specifiedas 𝑑 = 2 , and the delays are given as functions: 𝜏 𝑓 ( 𝑥 , 𝑝 ) = 𝑝 , 𝜏 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = 𝑝 + 𝑝 𝑥 + 𝑝 𝑥 , 𝑓 ( 𝑥 , 𝑥 , 𝑥 , 𝑝 ) = ⎡⎢⎢⎣ − 𝑝 tanh[ 𝑝 𝑥 ] + 𝑝 𝑥 𝑥 − 2 𝜋𝑥 − 𝑥 (( 𝑥 ) + ( 𝑥 ) ) 𝑥 + 2 𝜋𝑥 − 𝑥 (( 𝑥 ) + ( 𝑥 ) ) ⎤⎥⎥⎦ .Note that the delayed argument appearing in 𝑓 is now 𝑥 , instead of 𝑥 as was the casein (27). The question is what effects the two types of state dependence have on theobserved dynamics as represented by the bifurcation set in the ( 𝑐, ̄𝜏 ) -plane. This wasconsidered by Keane et al. (2019) for ranges of 𝜂 𝑒 and 𝜂 𝑐 up to 𝜂 𝑒 = 0 . and 𝜂 𝑐 = 0 . ,that is, twice their nominal values. It turns out that, within the ranges of the parametersconsidered, state dependence of 𝜏 𝑒 alone has a negligible effect on the bifurcation set.State dependence of 𝜏 𝑐 , on the other hand, has a significant impact on the bifurcationset in the ( 𝑐, ̄𝜏 ) -plane, featuring considerably increased and more overlapping resonanceregions. Surprisingly, for < 𝜏 𝑐 , state dependence of 𝜏 𝑒 does have a definite influence onthe bifurcation set, namely that of increasing the observed complexity even further.As an example of the effect of both types of state dependence, fig. 9 shows the bifur-cation set and maximal Lyapunov exponent of solutions in the ( 𝑐, ̄𝜏 ) -plane of (24) with(32) for the case where the upwelling 𝜂 𝑒 and the ocean adjustment 𝜂 𝑐 are at the maximumof their considered ranges at 𝜂 𝑒 = 0 . and 𝜂 𝑐 = 0 . , respectively. As for 𝜂 𝑒 = 𝜂 𝑐 = 0 infig. 7, shown in fig. 9 are curves SL, PD and T of saddle-node of periodic orbits, period-doubling and torus bifurcations, respectively. They were computed with DDE-BIFTOOL,thus, demonstrating that such computations can be performed readily also for DDEs thatfeature state dependence. The maximal Lyapunov exponent of solutions was computedagain for an upsweep of 𝑐 with the algorithm for DDEs from (Farmer, 1982). Note that inthe dark grey region for low values of ̄𝜏 the delay becomes negative during the integration,so that a sufficiently long time series to determine the Lyapunov exponent cannot befound. Similarly, some curves of period-doubling bifurcations stop in this region becausethe delay becomes negative along the respective periodic orbit during the continuation.Comparing fig. 9 with fig. 7 clearly drives home the point that state dependence has alarge effect on the bifurcation set and, hence, on the observable dynamics of the GZTmodel (24). In fig. 9 the bifurcation set now extends substantially further into the regionof large forcing 𝑐 : the curve T of torus bifurcation has moved, as have resonance regionsassociated with it. In particular, there is now a cluster of overlapping resonance regionsnear ̄𝜏 = 1 , that is, in the physically relevant range of the ( 𝑐, ̄𝜏 ) -plane. Notice thatthis cluster is associated with large positive Lyapunov exponents. More generally, statedependence results in considerably more and larger regions where chaotic dynamics canbe found. Interestingly, there is also a large region of chaotic dynamics for very low30 ¯ τ c • PPi
QQk (cid:27)
SL PD T
Figure 9:
Bifurcation set and maximal Lyapunov exponent of solutions in the ( 𝑐, ̄𝜏 ) -plane of (24)with (32) for 𝜂 𝑒 = 0 . and 𝜂 𝑐 = 0 . ; the blue dot indicates the parameter point for the timeseries in Figure 10. From [Keane et al. (2019)] © 2019 The Royal Society; reproduced withpermission. values of ̄𝜏 , near the boundary where the delay becomes negative.Clearly, introducing a physically motivated state-dependent delay time changes theoverall observed dynamics of the GZT model. We finish by demonstrating that this modelmodification leads to dynamics that represents realistic aspects of the ENSO system, moreso than those found in the absence of state dependence. Namely, in the constant delay case,irregular (chaotic) behaviour could only be found for small pockets of the ( 𝑐, ̄𝜏 ) -planeand, as such, could not be considered a prominent feature of the model behaviour. In thepresence of state dependence due to upwelling and ocean adjustment, on the other hand,this type of behaviour is more prominent, especially in the physically relevant clusterof period-doublings near ̄𝜏 = 1 . The blue dot in this cluster indicates the parameterpoint ( 𝑐, ̄𝜏 ) = (7 . , . , and fig. 10 shows the corresponding times series and powerspectra in direct comparison with those for the measured Nino3 index. This data is thespatially averaged SST over 5 ◦ N–5 ◦ S and 150 ◦ W–90 ◦ W as derived from the OptimumInterpolation SST V2 data by the National Oceanic and Atmospheric Administrationin Boulder, Colorado. The Nino3 time series, which was linearly detrended, is shownin fig. 10(a1). Prominent in the time series data is the strong annual forcing, which isrepresented by the large peak at 1 year in the power spectrum in panel (a2), which wascalculated by using the Welch method with windows of length 15 years and overlappingacross 12 years. Moreover, the Nino3 time series shows characteristic larger maxima, thatis, El Niño events, about every 4 to 7 years, which give rise to the distinct but broad peakin the power spectrum that is centred near the frequency of about . years. Indeed,there is clearly a high degree of variability in the timing of larger maxima and, as we31
982 1987 1992 1997 2002 2007 2012 2017-202 0 1-3-2-10 N i n o3 i nd e x l og [ a . u .] (a1) (a2) h l og [ a . u .] (b1) (b2) t [years] freq [years − ] Figure 10:
Time series of the Nino3 index (a1) obtained from the observational data set NOAAOptimum Interpolation SST V2 (Jan. 1982 – Dec. 2017) by linear detrending and the correspond-ing power spectrum (a2); and times series (b1) and corresponding power spectrum (b2) for theattractor of (24) with (32) for 𝜂 𝑒 = 0 . and 𝜂 𝑐 = 0 . at the parameter point ( 𝑐, ̄𝜏 ) = (7 . , . indicated in Figure 9. Data for row (a) is provided by the Physical Sciences Division, EarthSystem Research Laboratory, NOAA, Boulder, Colorado. From [Keane et al. (2019)] © 2019The Royal Society; reproduced with permission. checked, they tend to be seasonally locked. As row (b) of fig. 10 shows, time series andpower spectrum of the thermocline deviation ℎ of the GZT model (24) with delays givenby (32) and 𝜂 𝑒 = 0 . and 𝜂 𝑐 = 0 . at ( 𝑐, ̄𝜏 ) = (7 . , . also possess these importantcharacteristics of ENSO. The solution from which the time series is derived evolves on achaotic attractor that lies at the intersection of several resonance tongues. The times seriesin panel (b1) clearly lacks certain aspects of the data in panel (a1), and it is not obvioushow exactly the thermocline deviation ℎ as described by the rather simple conceptualGZT model translates to an observable such as Nino3. Nevertheless, the time series of ℎ features irregularity in the form of relatively large peaks that occur every 2–7 years, with asimilar broad peak centered near the frequency of about . years ) as well as seasonallocking, which is very robust with respect to the choice of parameters. The number anddistribution of large maxima, on the other hand, depends on where the parameters pointis chosen to lie in the regions of overlapping 𝑝 ∶ 𝑞 resonance tongues. Moreover, therelative strengths between the peaks in the power spectrum, representing both seasonaland El Niño time scales can be influenced by the choice of the seasonal forcing strength 𝑐 . Overall, we conclude that solutions with fundamental ENSO characteristics can befound in appropriate regions of parameter space of the state-dependent GZT model asconsidered here. 32 Resonance phenomena in a scalar DDE with twostate-dependent delays
The previous section demonstrated that state dependence of delays can have a seriousimpact on the observed dynamics of a given DDE. On the other hand, the GZT modelfor ENSO features complicated dynamics already when the delays are constant. As wewill discuss now, state dependence of delays alone can create complicated nonlineardynamics, even when the constant-delay DDE has only trivial, linear dynamics. Thissurprising result was obtained in Calleja et al. (2017) for the scalar DDE 𝑢 ′ ( 𝑡 ) = − 𝛾𝑢 ( 𝑡 ) − 𝜅 𝑢 ( 𝑡 − 𝑎 − 𝑐 𝑢 ( 𝑡 )) − 𝜅 𝑢 ( 𝑡 − 𝑎 − 𝑐 𝑢 ( 𝑡 )) . (33)Here, < 𝛾 is the linear decay rate and ≤ 𝜅 , 𝜅 are the strengths of the two negativefeedback loops with the constant delay times < 𝑐 , 𝑎 and linear state dependence ofstrengths ≤ 𝑐 , 𝑐 . For 𝜅 = 𝜅 = 0 , this system is simply a linear scalar equationwhose solutions decay exponentially to the origin with rate < 𝛾 . For < 𝜅 , 𝜅 , on theother hand, (33) is a DDE with two negative feedback loops. For 𝑐 = 𝑐 = 0 this DDE islinear with the two fixed delays 𝑎 and 𝑎 and all trajectories of (33) decay to the originor blow up to infinity, depending on the values of 𝛾 , 𝜅 and 𝜅 ; see Bellman and Cooke(1963); Hale (1977); Hale and Verduyn Lunel (1993). In other words, the dynamics ofthe system without state dependence in the delay terms is indeed trivial.The situation is very different with state dependence, that is, for < 𝑐 , 𝑐 , in whichcase (33) may show a wide range of behaviours. The two-delay state-dependent DDE(33) was introduced by Humphries et al. (2012). It is a generalisation of the single-delaystate-dependent DDE, corresponding to setting 𝜅 = 0 , which was first introduced ina singularly perturbed form as an example problem in (Mallet-Paret et al., 1994) andconsidered extensively in (Mallet-Paret and Nussbaum, 2011a). A singularly perturbedversion of the two-delay state-dependent DDE (33) was studied by Humphries et al.(2016) and Kozyreff and Erneux (2013). Specifically, solutions near the singular Hopfbifurcations were considered by Kozyreff and Erneux (2013), while large amplitudesingular solutions are constructed and studied by Humphries et al. (2016). We reporthere on the work by Calleja et al. (2017) and consider (33) for 𝑎 < 𝑎 without loss ofgenerality. It was shown by Humphries et al. (2012) that the state-dependent delays cannever become advanced when 𝜅 < 𝛾 , which we assume from now on. Hence, for anyLipschitz contiuous initial condition the initial value problem given by (33) has a uniquesolution. Moreover, Humphries et al. (2012) showed that state dependence of the delayterms changes the dynamics in an essential way. In particular, although it is only linear,the state dependence of the delays for < 𝑐 , 𝑐 is responsible for nonlinearity in thesystem, and the dynamics of the DDE (33) is no longer linear. Given that it has twofeedback loops, the system is, colloquially speaking, potentially at least as complicatedas two coupled damped nonlinear oscillators.This realization was the starting point of the extensive bifurcation analysis of the twodelay state-dependent DDE (33) by Calleja et al. (2017), where more details can be found.The first step is to bring (33) into the form required by DDE-BIFTOOL. To this end, one33 .8 2 2.2 2.4 2.6 (a1) κ κ HH H H u T T u (a2) κ HH H H u T T u (b) κ κ HH T u T H H H u CCCW
CCW ? (cid:8)(cid:8)(cid:8)* (cid:8)(cid:8)* (cid:8)(cid:8)* Figure 11:
Local torus bifurcation curves 𝑇 and 𝑇 𝑢 of (33) emerging from a Hopf-Hopf bifur-cation point HH at the intersection of curves 𝐻 and 𝐻 𝑢 , as computed from the normal form(a1) and by numerical continuation (a2). The bifurcation diagram in the ( 𝜅 , 𝜅 ) -plane (b) showsresonance tongues connecting different points of resonance on 𝑇 and 𝑇 𝑢 , which are actually partof a single closed curve of torus bifurcation. From [Calleja et al. (2017)] © 2017 Society forIndustrial and Applied Mathematics; reproduced with permission. has to specify the number of delays, 𝑑 = 2 , and define 𝜏 𝑓 ( 𝑥 , 𝑝 ) = 𝑎 + 𝑐 𝑥 , 𝜏 𝑓 ( 𝑥 , 𝑥 , 𝑝 ) = 𝑎 + 𝑐 𝑥 , 𝑓 ( 𝑥 , 𝑥 , 𝑥 , 𝑝 ) = − 𝛾𝑥 − 𝜅 𝑥 − 𝜅 𝑥 .Note that here the physical space is one-dimensional since 𝑥 ( 𝑡 ) ∈ ℝ , 𝑀 = 1 and theparameter vector is 𝑝 = ( 𝑝 , … , 𝑝 ) = ( 𝛾, 𝜅 , 𝜅 , 𝑎 , 𝑎 , 𝑐 , 𝑐 ) ∈ ℝ . We focus here on resonance phenomena associated with a point HH , where a Hopf-Hopf bifurcation occurs. As in (Humphries et al., 2012; Calleja et al., 2017), we fix theparameters of (33) to 𝛾 = 4 . , 𝑎 = 1 . , 𝑎 = 6 , and 𝑐 = 𝑐 = 1 . Thus, parameters34 and 𝜅 are bifurcation parameters, where we restrict to 𝜅 ∈ (0 , . so that indeed 𝜅 < 𝛾 . The bifurcation diagram of (33) in the ( 𝜅 , 𝜅 ) -plane is shown in fig. 11; hererow (a) focuses on the immediate vicinity of the point HH , while panel (b) shows therelevant bifurcation set associated with HH over a wider range of 𝜅 and 𝜅 .An important contribution of (Calleja et al., 2017) is the computation of the third-ordernormal form of the state-dependent DDE in the form of an ordinary differential equation(ODE) on the center manifold near Hopf-Hopf points; see, for example, (Guckenheimerand Holmes, 1983; Kuznetsov, 2013) for the respective ODE normal forms. This isachieved by expanding the state dependence to derive a DDE with terms up to a givenorder and with only constant delays. This constant-delay DDE can then be reduced to therequired four-dimensional ODE normal form with standard techniques; see (Bélair andCampbell, 1994; Guo and Wu, 2013; Wage, 2014). Computation with DDE-BIFTOOLshows that the point HH lies at ( 𝜅 , 𝜅 ) = (2 . , . where the Hopf bifurcationcurves 𝐻 and 𝐻 𝑢 intersect; it features a double pair of purely complex conjugate eigen-values with frequencies (imaginary parts) 𝜔 = 2 . and 𝜔 = 1 . . The normalform computation at this Hopf-Hopf point, details of which can be found in (Callejaet al., 2017), shows that HH is subcase III of what is referred to as the simple case in(Kuznetsov, 2013). This means that there are two curves of torus (or Neimark-Sacker)bifurcations emerging from the codimension-two Hopf-Hopf point. When the normalform coordinates are transformed back into the ( 𝜅 , 𝜅 ) -plane, one obtains the bifurcationdiagram shown in fig. 11(a1), featuring the curves 𝐻 and 𝐻 𝑢 and the torus bifurcationcurves 𝑇 and 𝑇 𝑢 . Note that all curves are straight lines, whose slopes are determined bythe respective normal form coefficient. The bifurcation diagram in fig. 11(a2) shows thesame bifurcation curves 𝐻 , 𝐻 𝑢 , 𝑇 and 𝑇 𝑢 but now computed for (33) by continuationwith DDE-BIFTOOL. Note that these curves are no longer straight lines. Comparisonwith panel (a1) shows that the nature, order and slopes of the respective bifurcation curvesis indeed as determined by the normal form computation, which strongly supports thecorrectness of the expansion method used to derive the Hopf-Hopf normal form of thefull state-dependent DDE (33). This approach has been extended by Sieber (2017) toall codimension-two bifurcations of steady-states that are defined by conditions on thelinearization. Hence, normal form calculations for these codimension-two bifurcations,which had been incorporated into the capabilities of the package DDE-BIFTOOL forconstant-delay DDEs by Wage (2014), are now also available for state-dependent DDEs;see (Sieber et al., 2015).fig. 11(b) shows that the local curves 𝑇 and 𝑇 𝑢 , when continued beyond a neighborhoodof the Hopf-Hopf point HH , actually form a single curve in the ( 𝜅 , 𝜅 ) -plane. Asexpected from theory, along the branches 𝑇 and 𝑇 𝑢 of torus bifurcations one finds pointsof 𝑝 ∶ 𝑞 resonance, which we include for 𝑞 ≤ . At each such point the Floquet multiplieris a rational multiple of 𝜋 and a resonance tongue emerges where the dynamics on thetorus is 𝑝 ∶ 𝑞 locked. For parameter points that do not lie in a resonance tongue the rotationnumber is an irrational multiple 𝛼 of 𝜋 and the dynamics of the torus is quasiperiodic.In either case, the bifurcating torus is normally hyperbolic, and hence smooth, near therespective torus bifurcation. Each resonance region is bounded locally near the point of 𝑝 ∶ 𝑞 resonance by a pair of saddle-node bifurcations of periodic orbit. Tori with fixedirrational rotation number 𝛼 , on the other hand, lie on smooth curves that connect to35 (c) u t − a ( θ ) u t − a ( θ ) θ (a) Σ −6 −4 −2 0−0.500.5 (b) u t ( θ ) θ Figure 12:
Normally hyperbolic quasiperiodic torus of (33) for 𝜅 = 4 . and 𝜅 = 3 . inprojection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space (a), represented by a single trajectory (light blue)together with the Poincaré trace (blue dots) on the (projected) section Σ (green); also shown arethe corresponding function segments represented by 𝑢 𝑡 ( 𝜃 ) (b) and by ( 𝑢 𝑡 − 𝑎 ( 𝜃 ) , 𝑢 𝑡 − 𝑎 ( 𝜃 )) . From[Calleja et al. (2017)] © 2017 Society for Industrial and Applied Mathematics; reproduced withpermission. the point on the torus bifurcation curve with the corresponding Floquet multiplier. Alsoshown in fig. 11(b) are the bounding curves of saddle-node bifurcations for the resonancepoints with 𝑞 ≤ . They have been found by identifying, by means of numericalintegration, stable locked periodic orbits near the respective branch of torus bifurcation;this is possible because the tori bifurcating from 𝑇 and 𝑇 𝑢 are actually attracting (whichis in agreement with the normal form calculation). The subsequent continuation of thesewith DDE-BIFTOOL in 𝜅 identifies the pair of saddle-node bifurcations, at two specificvalues of 𝜅 , that form the boundary of the resonance tongue. Once these two points ofsaddle-node bifurcations of periodic orbits have been found, they can be continued inboth 𝜅 and 𝜅 towards the curves 𝑇 and 𝑇 𝑢 to obtain the respective curves shown in the ( 𝜅 , 𝜅 ) -plane. Note that the gap between the two bounding curves of the 𝑝 ∶ 𝑞 resonancetongue becomes smaller for increasing 𝑞 . As we show now, invariant tori of DDEs can also be computed, including when the delayis state dependent; see Krauskopf and Green (2003); Green et al. (2003); Calleja et al.(2017). This is quite straigthforward for an attracting quasiperiodic torus, because it36s densely filled by any trajectory on it. Such a trajectory can be obtained readily bynumerical integration from an initial condition sufficiently near the torus, after transientshave settled down. Such a torus is a smooth two-dimensional submanifold that lives inthe infinite-dimensional phase space 𝐶 of the DDE. Therefore, the question is how torepresent it via a suitable low-dimensional projection. fig. 12 shows with the example ofthe smooth attracting quasiperiodic torus of (33) for 𝜅 = 4 . and 𝜅 = 3 . how this canbe achieved. Panel (a) shows the computed long trajectory on the torus in the natural andconvenient projection onto the three-dimensional ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space of (33).Also shown is the intersection set of the trajectory on the torus with the shown section Σ ,which forms a smooth invariant curve as is expected for a quasiperiodic torus.Note that the representation of the torus in fig. 12(a) looks very much like a two-dimensional smooth torus in a three-dimensional phase space of an ODE, with an as-sociated image of the dynamics of the Poincaré map in the two-dimensional section Σ .However, it is important to recognize that this image is a projection from the infinite-dimensional phase space 𝐶 . In particular, it is an interesting question how best to definea Poincaré map for a DDE. In general terms, given a section Σ of codimension one in thephase space 𝐶 that is transverse in some region of interest to the (semi)flow Φ 𝑡 generatedby the DDE, the (local) Poincaré map 𝑃 is defined as 𝑃 Σ ∶ Σ → Σ ,𝑞 ↦ Φ 𝑡 𝑞 ( 𝑞 ) , (34)where 𝑡 𝑞 > is the return time to Σ . The main issue from a practical perspective is how todefine the section Σ . When the DDE has a physical space ℝ 𝑛 of sufficient dimension (atleast three), then it is convenient to consider a codimension-one section Σ ⊂ ℝ 𝑛 ; requiringthat the headpoint 𝑞 (0) of the point 𝑞 lies in Σ induces a codimension-one section inthe infinite dimensional phase space 𝐶 , which we also refer to as Σ for simplicity; see(Krauskopf and Green, 2003). Moreover, there is a natural projection onto ℝ 𝑛 and justconsidering the headpoints in the section Σ ⊂ ℝ 𝑛 gives what we refer to as the finite-dimensional Poincaré trace of the dynamics. Unfortunately, this approach is not workablefor (33) because it is a scalar DDE and, moreover, state dependent. Instead, we make useof the fact that all periodic and quasi-periodic orbits cross repeatedly cross { 𝑢 = 0} ⊂ ℝ when the parameters are all positive as considered here; therefore, it is natural to use thiscondition for defining the Poincaré section as Σ = { 𝑞 ∈ 𝐶 ∶ 𝑞 (0) = 0} . (35)Clearly, the section Σ is infinite dimensional itself, and the local Poincaré map 𝑃 Σ on Σ is defined as the map that takes a downward transversal crossing of zero (where 𝑞 (0) = 0 with 𝑞 ′ (0) < ) to the next such crossing. The Poincaré trace in the 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane, which is the projection of the infinite-dimensional Σ from (35), is obtained fromthe projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space simply by also requiring that 𝑢 ( 𝑡 ) = 0 .This projected section, which we again refer to as Σ for convenience, is shown in fig. 12(a).The underlying projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space generalises an idea ofMackey and Glass (1977), who were the first to project solutions of DDEs into finitedimensions by plotting values of 𝑢 ( 𝑡 − 𝜏 ) against 𝑢 ( 𝑡 ) for a single-delay DDE. Defining the37 (c) u t − a ( θ ) u t − a ( θ ) θ (a1) (a2) Σ −6 −4 −2 0−0.500.5 (b) u t ( θ ) θ Figure 13:
Normally hyperbolic phase-locked torus of (33) for 𝜅 = 5 . and 𝜅 = 2 . in projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space (a1), represented by the stable periodic orbit(blue), the saddle periodic orbit (red), and its unstable manifold (grey curves), together with thePoincaré trace on the (projected) section Σ (green) (a2); also shown are the corresponding functionsegments represented by 𝑢 𝑡 ( 𝜃 ) (b) and by ( 𝑢 𝑡 − 𝑎 ( 𝜃 ) , 𝑢 𝑡 − 𝑎 ( 𝜃 )) . From [Calleja et al. (2017)] © 2017Society for Industrial and Applied Mathematics; reproduced with permission. section Σ for (33) by 𝑞 (0) = 0 has the advantage that the two delays are exactly 𝑎 and 𝑎 ,that is, constant. fig. 12(b) and (c) illustrates this by showing the function segments of allthe points that generate the Poincaré trace on Σ ; tey are represented by 𝑢 𝑡 ( 𝜃 ) in panel (b)and by ( 𝑢 𝑡 − 𝑎 ( 𝜃 ) , 𝑢 𝑡 − 𝑎 ( 𝜃 )) in panel (c), both as functions of the argument 𝜃 , which runsover the interval [−6 , since 𝑎 = 𝑎 = 6 . In particular, fig. 12(c) clearly shows the‘history tails’ over the time interval [−6 , associated with headpoints that form the tracein (the two-dimensional projection of) Σ (given by 𝜃 = 0 ); see also (Krauskopf and Green,2003).fig. 13 shows an example of an attracting smooth invariant torus with locked dynamics;namely this example for 𝜅 = 5 . and 𝜅 = 2 . is from the resonance tongue.Hence, there are a stable periodic orbit and a saddle periodic peridic orbit that bothform a torus knot. The presentation is as for the quasi-periodic torus in fig. 12.38anel (a1) of fig. 13 shows the shows the torus rendered as a surface in projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space, together with the (projection of the) section Σ . ThePoincaré trace in the ( 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane, that is, the intersection set of the toruswith Σ , is shown on its own in panel (a2). The associated function segments or ‘historytails’ are shown in fig. 13(c) and (d) as represented by 𝑢 𝑡 ( 𝜃 ) and by ( 𝑢 𝑡 − 𝑎 ( 𝜃 ) , 𝑢 𝑡 − 𝑎 ( 𝜃 )) ,respectively, for 𝜃 ∈ [−6 , . The locked dynamics on the torus as represented infig. 13(a) is again very reminiscent of what one would expect to find for a torus of athree-dimensional ODE: its two-dimensional Poincaré trace in panel (a2) clearly shows asingle smooth curve with four points of a stable period-four orbit and four points of anunstable period-four orbit; see fig. 13(a2). Notice that the invariant curve has a point ofself-intersection; this is due to projection and a reminder that we are dealing with a DDEwith an infinite-dimensional phase space.As opposed to the case of a quasiperiodic torus, a torus with locked dynamics cannot befound by numerical integration alone. Indeed, any initial condition will, after transientshave settled down, trace out only the attracting periodic orbit. The torus on which it liescan be computed as follows. Continuation of the stable periodic orbit in the parameter 𝜅 gives, after a fold or saddle-node bifurcation of periodic orbits, the coexisting saddleperiodic orbit for the intitial value of 𝜅 = 5 . . As theory predicts, this saddle peridicorbit has one unstable Floquet multiplier and, hence, one unstable eigenfunction, whichwe extracted from the DDE-BIFTOOL data; see also (Green et al., 2003). We then usedthe eigenfunction to define two initial functions in the local unstable manifold of theperiodic orbit, one on each side and sufficiently cloese to of the saddle periodic orbit.Trajectory segments that lie on the unstable manifold were then found with numericalintegration from initial functions along the unstable eigenfunction; a selection of them isshown in fig. 13(b) and (c). The torus was rendered as a surface in ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space in panel (a1) and as an invariant curve in the ( 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane in panel (a2)by ordering a suitable selection of trajectory segments from the Poincaré section back toitself. As was already mentioned in Section 5.1, any computed pair of saddle-node bifurcationcurves shown in fig. 11(b) emerges at a resonance point on the torus bifurcation curve 𝑇 𝑢 and connects to a point of resonance on the torus bifurcation curve 𝑇 (or vice versa).More generally, this is evidence for the observation that, near the Hopf-Hopf point HH ,any 𝑝 ∶ 𝑞 resonance point on the upper branch 𝑇 𝑢 is connected by a pair of saddle-nodebifurcation curves with a 𝑝 ∶ ( 𝑝 + 𝑞 ) resonance point on the lower branch 𝑇 . In the regionof locked dynamics bounded by such a pair one, hence, finds 𝑝 ∶ 𝑞 locked dynamics on asmooth invariant torus near 𝑇 𝑢 and 𝑝 ∶ ( 𝑝 + 𝑞 ) locked dynamics on a smooth invarianttorus near 𝑇 . Since a 𝑝 ∶ 𝑞 torus knot and a 𝑝 ∶ ( 𝑝 + 𝑞 ) torus knot cannot exist onone and the same smooth torus for topological reasons, the torus inside the respectiveresonance tongue cannot be smooth throughout in the transition inside the tongue fromthe 𝑝 ∶ 𝑞 to the 𝑝 ∶ ( 𝑝 + 𝑞 ) resonance point (or vice versa). On the other hand, withoutthe requirement that the periodic orbit lies on an invariant two-dimensional torus, it ispossible to transform a periodic orbit with 𝑞 loops into one with ( 𝑝 + 𝑞 ) loops; in other39 (a) κ κ T u T SL • • (b) κ u ( t − a ) (c) u ( t − a ) u ( t − a ) u ( t ) (d) u ( t − a ) u ( t − a ) Figure 14:
The resonance tongue in the ( 𝜅 , 𝜅 ) -plane of (33) that connects the resonanceon 𝑇 𝑢 with the resonance on 𝑇 (a); the one-parameter bifurcation diagram in 𝜅 for fixed 𝜅 = 3 . (b) showing the values of 𝑢 ( 𝑡 − 𝑎 ) of the Poincaré traces of stable (blue) and saddle(red) periodic orbits (red) inside the resonance tongue, and of other solutions on tori outside theresonance tongue; the phase-locked torus-like object (grey) for 𝜅 = 6 . with the stableand saddle periodic orbits in projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space (c); and its Poincarétrace in the ( 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane. From [Calleja et al. (2017)] © 2017 Society for Industrialand Applied Mathematics; reproduced with permission. words, there is indeed no topological obstruction for the bounding curves of saddle-nodebifurcations of periodic orbits to connect in the way we found in fig. 11(b).fig. 14 shows that it is possible to perform computations that show how an invariantlocked torus of a DDE loses its smoothness and bifurcates further. The key idea here isto compute the one-dimensional unstable manifold in the Poincaré trace that is associatedwith a saddle-periodic orbit with a single unstable Floquet multiplier. In other words, theapproach we used to compute the smooth phase-locked torus in fig. 13 also workswhen the torus is no longer smooth; the only requirement is that the saddle-periodic orbitstill exists and has a single unstable Floquet multiplier; see also (Krauskopf and Green,2003; Green et al., 2003; Calleja et al., 2017; Keane and Krauskopf, 2018). Panel (a) offig. 14 shows the resonance tongue in the ( 𝜅 , 𝜅 ) -plane that connects the resonanceon 𝑇 𝑢 with the resonance on 𝑇 . The line segment at 𝜅 = 3 . indicates the 𝜅 -rangeof the one-parameter bifurcation diagram shown in panel (b). Specifically, shown is thevalue of 𝑢 ( 𝑡 − 𝑎 ) when 𝑢 ( 𝑡 ) = 0 , that is, the first component of the ( 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane of the Poincaré trace. Inside the resonance tongue we find locking: there arefour branches of the stable periodic orbit and four branches of the saddle periodic orbit,which meet and disappear in saddle-node bifurcations; these branches were found withDDE-BIFTOOL by continuation in 𝜅 , and this computation also confirms that the saddle40eriodic orbit has a single unstable Floquet multiplier throughout the 𝜅 -range of theresonance tongue when 𝜅 = 3 . . Outside the resonance tongue we find quasi-periodicdynamics or dynamics of very high period; it was identified by numerical integration andis represented by many points in the Poincaré trace whose 𝑢 ( 𝑡 − 𝑎 ) -values effectively fillout intervals.Panels (c) and (d) of fig. 14 show the result of computing the unstable manifold ofthe saddle periodic orbit for 𝜅 = 6 . with the approach from Section 5.2. Whilethis is hard to see in the projection onto ( 𝑢 ( 𝑡 ) , 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -space in fig. 14(c),the Poincaré trace in the ( 𝑢 ( 𝑡 − 𝑎 ) , 𝑢 ( 𝑡 − 𝑎 )) -plane in panel (d) clearly shows that theone-dimensional unstable manifold of the saddle periodic orbit now spirals around thestable periodic orbit. This means that the stable periodic orbit has two dominant Floquetmultipliers that are complex conjugate, which is confirmed by the computation of theFloquet multipliers during the continuation of the periodic orbits with DDE-BIFTOOL.The attracting periodic orbit on the torus developing a pair of complex conjugate Floquetmultiplier is a mechanism for the loss of normal hyperbolicity of an invariant torus thatis known from ODE theory (Aronson et al., 1982). As fig. 14(c) and (d) shows, there isstill a continuous two-dimensional torus, formed by the closure of this unstable manifold,but this torus is indeed no longer smooth.Further bifurcations may occur that change the nature of the invariant set in the Poincarétrace, including homoclinic and heteroclinic tangencies of unstable manifold of saddleperiodic orbits. This happens, for example, when the resonance tongue is crossedagain at 𝜅 = 3 . for larger values of 𝜅 ; see Calleja et al. (2017) for the details. Otherexamples where such global bifurcations in DDEs have been identified via unstablemanifold computations are the transition to chaos in a laser with phase-conjugate feedbackin (Krauskopf and Green, 2003; Green et al., 2003) and the break-up of a torus in theGZT model of Section 4.2 due to the transition through a bifurcation structure known asa Chenciner bubble in (Keane and Krauskopf, 2018). The case studies of the GZT ENSO model and of the DDE with two state-dependent delayswe presented show that core tasks of numerical bifurcation analysis can be performedfor DDEs with finitely many discrete delays, even when the delays are state dependent.More specifically, the routines that are implemented within the package DDE-BIFTOOLinclude the detection and continuation of equilibria, periodic orbits and their bifurcationsof codimension one, of codimension-one connecting orbits between equilibria, as well asthe computation of normal forms of bifurcations of equilibria up to including codimensiontwo. This suite of tools puts the present capabilities practically at the same level that isavailable for ODEs. In other words, the numerical bifurcation analysis of DDEs, whetherthey arise in applications or in a theory context, is now perfectly feasible.One area where the capabilities for DDEs still lag behind that for ODEs is the compu-tation of normal forms for bifurcation of periodic orbits. The approach to normal formanalysis designed by Dhooge et al. (2003) for MATCONT and extended by Bosschaertet al. (2020) to equilibria of DDEs with constant delays is, in principle, applicable also41o periodic orbits and their bifurcations. However, there remain some technical diffcul-ties. For example, in case delays are bifurcation parameters or are state dependent, thecomputation of normal form coefficients may involve the computation of high-order timederivatives of the piecewise polynomials representing the solution.We considered here chiefly DDEs in standard form with a finite number of discretedelays. For this class the discussed tools for the numerical bifurcation analysis are onvery firm ground theoretically when the delays are constant. On the other hand, somepresent capabilities of the numerical methods and the software assume properties ofthe underlying DDE that have not yet been proven rigorously when the delays are statedependent. For example, convergence of the collocation schemes used for the representingperiodic orbits has been proved for standard DDE with constant delays, but remains anopen question when the delays are part of the unknowns or state-dependent. Case studiessuch as the ones presented here clearly suggest that collocation ‘works well’ also insuch wider circumstances; moreover, the techniques introduced by Andó and Breda(2020) look promising as a tool for proving this. Similarly and as we also demonstrated,associated normal form calculations appear to be working perfectly fine when delays arestate dependent and are in agreement with the results of numerical bifurcation analysis.Yet the proof that the suggested expansion of the state dependence gives the correctnormal form is still outstanding — owing to the fact that regularity results for local centermanifolds in DDEs with state-dependent delays are strictly speaking still open. In spiteof these technical difficulties, we would argue that the tools we presented can be usedwith considerable confidence also for DDEs with discrete state-dependent delays.The methods as implemented in DDE-BIFTOOL actually permit the bifurcation analy-sis of systems from a far larger class of problems, including neutral DDEs (with constantor state-dependent delays), differential algebraic equations with delay, possibly of higherindex, and advanced-delayed systems. We explained briefly how these types of systemscan be defined within the framework of the software, so that the different tasks of bifur-cation theory can be performed also for such DDEs that are not in standard form (with anon-identity matrix multiplying the left-hand side). However, rigorous regularity results(such as the existence of smooth local center manifolds or branches of periodic orbits)and numerical convergence statements are not available yet for many of these problems.Therefore, when attempting a numerical bifurcation analysis in this wider context it ispresently the responsibility of the user to experiment and test for convergence a-posteriori.Finally, we hope that this review may encourage the use of numerical tools frombifurcation theory in the study of systems with delays in different application contexts.In particular, we would like to stress again that these tools are available and reliable notonly when the delays are constant, but also for the case that delays are state dependent.Hence, there is no need to approximate state-dependent delays with constant delays. Thismessage is important from a practical perspective because state dependence may beresponsible for layers of additional dynamics. Indeed, as we have demonstrated, in somesituations it may even generate all of the nontrivial dynamics. Case studies of specificDDEs beyond the standard form would also be very interesting and are encouraged. Whilethis is more challenging in terms of ensuring that the results stand up to scrutiny, suchinvestigations have a role in guiding the further development of theory and methods — inmuch the same way as case studies of standard DDEs with constant and state-dependent42elays have helped us get to where we are now.
Acknowledgements
We are very grateful to Dimitri Breda for organizing the school
Controlling DelayedDynamics: Advances in Theory, Methods and Applications at the International Centrefor Mechanical Sciences (CISM) in November 2019 in Udine, and we thank CISM forproviding location, management and support. A big thank you also to the other presenters,Tamás Insperger, Wim Michiels, Silviu Niculescu, Sjoerd Verduyn Lunel and of courseDimitri himself, and to all participants for making the school such an enjoyable andengaging event. As this volume shows, the school achieved its goal of providing anup-to-date snapshot of the field. A lot of credit goes again to Dimitri for pulling thisvolume together in the very challenging time of a global pandemic.The work reviewed in this chapter is the result of collaborations with colleagues andfriends over quite a number of years. In particular, we would like to acknowledge thecontributions of Renato Calleja, Henk Dijstra, Tony Humphries, Andrew Keane andClaire Postlethwaite to the two case studies presented here. Moreover, we thank Elsevier,the Society for Industrial and Applied Mathematics, and The Royal Society for permissionto reproduce figures and accompanying text passages from our previously published work.The research of BK was supported by Royal Society Te Ap¯arangi Marsden Fund grant
References
K. A. Abell, C. E. Elmer, A. R. Humphries, and E. S. Van Vleck. Computation ofmixed type functional differential boundary value problems.
SIAM Journal on AppliedDynamical Systems , 4(3):755–781, 2005.A. Andó and D. Breda. Convergence analysis of collocation methods for computingperiodic solutions of retarded functional differential equations.
SIAM Journal onNumerical Analysis , 2020. to appear, arXiv: 2008.07604, arxiv.org/abs/2008.07604 .D. G. Aronson, M. A. Chory, G. R. Hall, and R. P. McGehee. Bifurcations from aninvariant circle for two-parameter families of maps of the plane: a computer-assistedstudy.
Communications in Mathematical Physics , 83(3):303–354, 1982.K. Bar-Eli and R. J. Field. Earth-average temperature: a time delay approach.
Journal ofGeophysical Research: Atmospheres , 103(D20):25949–25956, 1998.A. G. Barnston, M. K. Tippett, M. L. L’Heureux, S. Li, and D. G. DeWitt. Skill of real-time seasonal ENSO model predictions during 2002–11: is our capability increasing?
Bulletin of the American Meteorological Society , 93(5):631–651, 2012.43. A. W. Barton, B. Krauskopf, and R. E. Wilson. Collocation schemes for periodicsolutions of neutral delay differential equations.
Journal of Difference Equations andApplications , 12(11):1087–1101, 2006.D. S. Battisti and A. C. Hirst. Interannual variability in a tropical atmosphere-oceanmodel: influence of the basic state, ocean geometry and nonlinearity.
Journal of theAtmospheric Sciences , 46(12):1687–1712, 1989.J. Bélair and S. A. Campbell. Stability and bifurcations of equilibria in a multiple-delayeddifferential equation.
SIAM Journal on Applied Mathematics , 54(5):1402–1424, 1994.R. E. Bellman and K. L. Cooke.
Differential-Difference Equations . Academic Press,1963.W. Beyn. The numerical computation of connecting orbits in dynamical systems.
IMAJournal of Numerical Analysis , 3:379–405, 1990.F. Borgioli, D. Hajdu, T. Insperger, G. Stépán, and W. Michiels. Pseudospectral methodfor assessing stability robustness for linear time-periodic delayed dynamical systems.
International Journal for Numerical Methods in Engineering , 121(16):3505–3528,2020.M. M. Bosschaert, S. G. Janssens, and Yu. A. Kuznetsov. Switching to nonhyperboliccycles from codimension two bifurcations of equilibria of delay differential equations.
SIAM Journal on Applied Dynamical Systems , 19(1):252–303, 2020.J.-P. Boulanger and C. Menkes. Propagation and reflection of long equatorial waves inthe Pacific Ocean during the 1992–1993 El Niño.
Journal of Geophysical Research:Oceans , 100(C12):25041–25059, 1995.D. Breda, S. Maset, and R. Vermiglio. Pseudospectral differencing methods for charac-teristic roots of delay differential equations.
SIAM Journal on Scientific Computing ,27(2):482–495, 2005.H. Broer, C. Simó, and J. C. Tatjer. Towards global models near homoclinic tangenciesof dissipative diffeomorphisms.
Nonlinearity , 11(3):667, 1998.R. C. Calleja, A. R. Humphries, and B. Krauskopf. Resonance phenomena in a scalardelay differential equation with two state-dependent delays.
SIAM Journal on AppliedDynamical Systems , 16(3):1474–1513, 2017.M. A. Cane, M. Münnich, and S. F. Zebiak. A study of self-excited oscillations of thetropical ocean-atmosphere system. Part I: Linear analysis.
Journal of the AtmosphericSciences , 47(13):1562–1577, 1990.D. B. Chelton and M. G. Schlax. Global observations of oceanic Rossby waves.
Science ,272(5259):234–238, 1996. 44. Dankowicz and F. Schilder.
Recipes for Continuation . Computer Science and En-gineering. SIAM, 2013. COCO download: http://sourceforge.net/projects/cocotools .A. Dhooge, W. J. F. Govaerts, and Yu. A. Kuznetsov. MatCont: A Matlab package fornumerical bifurcation analysis of ODEs.
ACM Transactions on Mathematical Software ,29(2):141–164, 2003. https://sourceforge.net/projects/matcont/ .O. Diekmann, S. van Gils, S. M. Verduyn Lunel, and H.-O. Walther.
Delay Equations ,volume 110 of
Applied Mathematical Sciences . Springer-Verlag, 1995.H. A. Dijkstra.
Dynamical Oceanography . Springer Science & Business Media, 2008.H. A. Dijkstra.
Nonlinear Climate Dynamics . Cambridge University Press, 2013.E. J. Doedel. Lecture notes on numerical analysis of nonlinear equations. In B. Krauskopf,H. M. Osinga, and J. Galán-Vioque, editors,
Numerical Continuation Methods forDynamical Systems: Path Following and Boundary Value Problems , pages 1–49.Springer-Verlag, Dordrecht, 2007.E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, andX. Wang.
AUTO97: Continuation and bifurcation software for ordinary differentialequations . URL http://indy.cs.concordia.ca/auto/ , 1999.J. W. Eaton, D. Bateman, S. Hauberg, and R. Wehbring.
GNU Octave version 4.2.1manual: a high-level interactive language for numerical computations , 2017. URL .K. Engelborghs and E. J. Doedel. Stability of piecewise polynomial collocation forcomputing periodic solutions of delay differential equations.
Numerische Mathematik ,91(4):627–648, 2002.K. Engelborghs and D. Roose. Numerical computation of stability and detection ofHopf bifurcations of steady state solutions of delay differential equations.
Advances inComputational Mathematics , 10(3–4):271–289, 1999.K. Engelborghs and D. Roose. On stability of lms methods and characteristic roots ofdelay differential equations.
SIAM Journal on Numerical Analysis , 40(2):629–650,2002.K. Engelborghs, T. Luzyanina, K. J. in ´t Hout, and D. Roose. Collocation methods forthe computation of periodic solutions of delay differential equations.
SIAM Journal onScientific Computing , 22:1593–1609, 2000a.K. Engelborghs, T. Luzyanina, and G. Samaey. DDE-BIFTOOL: a Matlab packagefor bifurcation analysis of delay differential equations. Report TW 305, KatholiekeUniversiteit Leuven, 2000b. 45. Engelborghs, T. Luzyanina, and G. Samaey. DDE-BIFTOOL v.2.00: a Matlab packagefor bifurcation analysis of delay differential equations. Report TW 330, KatholiekeUniversiteit Leuven, 2001.S. K. J. Falkena, C. Quinn, J. Sieber, J. Frank, and H. A. Dijkstra. Derivation of delayequation climate models using the Mori-Zwanzig formalism.
Proceedings of the RoyalSociety A , 475(2227), 2019.J. D. Farmer. Chaotic attractors of an infinite-dimensional dynamical system.
Physica D ,4(3):366–393, 1982.M. Ghil, I. Zaliapin, and S. Thompson. A delay differential model of ENSO variabil-ity: parametric instability and the distribution of extremes.
Nonlinear Processes inGeophysics , 15(3):417–433, 2008.W. J. F. Govaerts.
Numerical Methods for Bifurcations of Dynamical Equilibria . Miscel-laneous Titles in Applied Mathematics Series. SIAM, 2000.N. E. Graham and W. B. White. The El Niño cycle: a natural oscillator of the Pacificocean-atmosphere system.
Science , 240:1293–âĂŞ1302, 1988.K. Green, B. Krauskopf, and K. Engelborghs. One-dimensional unstable eigenfunctionand manifold computations in delay differential equations.
Journal of ComputationalPhysics , 197(1):86–98, 2003.J. Guckenheimer and P. Holmes.
Nonlinear Oscillations, Dynamical Systems, andBifurcations of Vector Fields . Springer Verlag, 1983.S. Guo and J. Wu.
Bifurcation Theory of Functional Differential Equations . AppliedMathematical Sciences. Springer-Verlag, 2013.B. Haegeman, K. Engelborghs, D. Roose, D. Pieroux, and T. Erneux. Stability and ruptureof bifurcation bridges in semiconductor lasers subject to optical feedback.
PhysicalReview E , 66:046216, 2002.J. K. Hale.
Theory of Functional Differential Equations . Applied Mathematical Sciences.Springer-Verlag, 1977.J. K. Hale and S. M. Verduyn Lunel.
Introduction to Functional Differential Equations .Applied Mathematical Sciences. Springer-Verlag, 1993.D. E. Harrison and G. A. Vecchi. El Niño and La Niña equatorial Pacific thermoclinedepth and sea surface temperature anomalies, 1986–98.
Geophysical Research Letters ,28(6):1051–1054, 2001.F. Hartung, T. Krisztin, H.-O. Walther, and J. Wu. Functional differential equationswith state-dependent delays: theory and applications. In P. Drábek, A. Cañada, andA. Fonda, editors,
Handbook of Differential Equations: Ordinary Differential Equa-tions , volume 3, chapter 5, pages 435–545. North-Holland, 2006.46. E. Henderson. Multiple parameter continuation: Computing implicitly defined k-manifolds.
International Journal of Bifurcation and Chaos , 12(03):451–476, 2002.A. R. Humphries, O. A. DeMasi, F. M. G. Magpantay, and F. Upham. Dynamics of a delaydifferential equation with multiple state-dependent delays.
Discrete and ContinuousDynamical Systems , 32(8):2701–2727, 2012.A. R. Humphries, D. A. Bernucci, R. Calleja, N. Homayounfar, and M. Snarski. Periodicsolutions of a singularly perturbed delay differential equation with two state-dependentdelays.
Journal of Dynamics and Differential Equations , 8:1215–1263, 2016.T. Insperger and G. Stépán. Stability of the milling process.
Periodica PolytechnicaMechanical Engineering , 44(1):47–57, 2000.T. Insperger, J. Milton, and G. Stépán. Acceleration feedback improves balancing againstreflex delay.
Journal of the Royal Society Interface , 10(79):20120763, 2013.F.-F. Jin. An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-downcoupled model.
Journal of the Atmospheric Sciences , 54(7):830–847, 1997.W. Just, B. Fiedler, M. Georgi, V. Flunkert, P. Hövel, and E. Schöll. Beyond the oddnumber limitation: A bifurcation analysis of time-delayed feedback control.
PhysicalReview E , 76:026210, 2007.M. A. Kaashoek and S. M. Verduyn Lunel. Characteristic matrices and spectral propertiesof evolutionary systems.
Transactions of the American Mathematicsl Society , 334(2):479–517, 1992.D. M. Kane and K. A. Shore.
Unlocking Dynamical Diversity: Optical Feedback Effectson Semiconductor Lasers . John Wiley & Sons, 2005.H. Kaper and H. Engler.
Mathematics and Climate . SIAM, 2013.A. Keane and B. Krauskopf. Chenciner bubbles and torus break-up in a periodicallyforced delay differential equation.
Nonlinearity , 31(6):R165, 2018.A. Keane, B. Krauskopf, and C. M. Postlethwaite. Delayed feedback versus seasonalforcing: Resonance phenomena in an El Niño Southern Oscillation model.
SIAMJournal on Applied Dynamical Systems , 14(3):1229–1257, 2015.A. Keane, B. Krauskopf, and C. M. Postlethwaite. Investigating irregular behavior in amodel for the El Niño Southern Oscillation with positive and negative delayed feedback.
SIAM Journal on Applied Dynamical Systems , 15(3):1656–1689, 2016.A. Keane, B. Krauskopf, and C. M. Postlethwaite. Climate models with delay differentialequations.
Chaos , 27(11):114309, 2017.A. Keane, B. Krauskopf, and H. A. Dijkstra. Delayed feedback versus seasonal forcing:Resonance phenomena in an El Niño Southern Oscillation model.
PhilosophicalTransactions of the Royal Society A , 377(2153):2153, 2019.47. Kozyreff and T. Erneux. Singular Hopf bifurcation in a differential equation withlarge state-dependent delay.
Procedings of the Royal Society A , 470(2162), 2013.B. Krauskopf. Bifurcation analysis of lasers with delay. In D. M. Kane and K. A. Shore,editors,
Unlocking Dynamical Diversity: Optical Feedback Effects on SemiconductorLasers , chapter 5, pages 147–183. John Wiley and Sons, Hoboken, NJ, USA, 2005.B. Krauskopf and K. Green. Computing unstable manifolds of periodic orbits in delaydifferential equations.
Journal of Computational Physics , 186(1):230–249, 2003.B. Krauskopf and J. Sieber. Bifurcation analysis of delay-induced resonances of theEl-Niño Southern Oscillation.
Procedings of the Royal Society A , 470(2169):20140348,2014.B. Krauskopf, G. H. M. Van Tartwijk, and G. R. Gray. Symmetry properties of laserssubject to optical feedback.
Optics Communications , 177(1–6):347–353, 2000.T. Krisztin. A local unstable manifold for differential equations with state-dependentdelay.
Discrete and Continuous Dynamical Systems , 9(4):993–1028, 2003.T. Krisztin. Smooth center manifolds for differential equations with state-dependentdelay. In
AIMS Conference Poitiers , 2006.Yu. A. Kuznetsov.
Elements of Applied Bifurcation Theory . Springer Science & BusinessMedia, 2013.M. C. Mackey and L. Glass. Oscillation and chaos in physiological control systems.
Science , 197(4300):287–289, 1977.J. Mallet-Paret and R. D. Nussbaum. Superstability and rigorous asymptotics in singu-larly perturbed state-dependent delay-differetnial equations.
Journal of DifferentialEquations , 250:4037–4084, 2011a.J. Mallet-Paret and R. D. Nussbaum. Stability of periodic solutions of state-dependentdelay-differential equations.
Journal of Differential Equations , 250(11):4085–4103,2011b.J. Mallet-Paret, R. D. Nussbaum, and P. Paraskevopoulos. Periodic solutions for functionaldifferential equations with multiple state-dependent time lags.
Topological Methods inNonlinear Analysis , 3:101–162, 1994.MATLAB. version 9.5.0.944444 (R2018b) . The MathWorks Inc., Natick, Massachusetts,2018.J. Milton, J. L. Townsend, M. A. King, and T. Ohira. Balancing with positive feedback:the case for discontinuous control.
Philosophical Transactions of the Royal Society A ,367(1891), 2009. 48. Münnich, M. A. Cane, and S. E. Zebiak. A study of self-excited oscillations of thetropical ocean-atmosphere system. Part II: Nonlinear cases.
Journal of the AtmosphericSciences , 48(10):1238–1248, 1991.D. Pieroux, T. Erneux, B. Haegeman, K. Engelborghs, and D. Roose. Bridges of periodicsolutions and tori in semiconductor lasers subject to delay.
Physics Letters , 87:193901,2001.Y. Pomeau and P. Manneville. Intermittent transition to turbulence in dissipative dynami-cal systems.
Communications in Mathematical Physics , 74(2):189–197, 1980.A. S. Purewal, C. M. Postlethwaite, and B. Krauskopf. A global bifurcation analysis ofthe subcritical Hopf normal form subject to Pyragas time-delayed feedback control.
SIAM Journal on Applied Dynamical Systems , 13(4):1879–1915, 2014.D. Roose and R. Szalai. Continuation and bifurcation analysis of delay differentialequations. In B. Krauskopf, H. M. Osinga, and J. Galán-Vioque, editors,
NumericalContinuation Methods for Dynamical Systems: Path Following and Boundary ValueProblems , pages 359–399. Springer-Verlag, Dordrecht, 2007.G. Samaey, K. Engelborghs, and D. Roose. Numerical computation of connecting orbitsin delay differential equations.
Numerical Algorithms , 30:335–352, 2002.J. Sieber. Finding periodic orbits in state-dependent delay differential equations as roots ofalgebraic equations.
Discrete & Continuous Dynamical Systems-A , 32(8):2607–2651,2012. updated with corrections on arXiv:1010.2391.J. Sieber. Local bifurcations in differential equations with state-dependent delay.
Chaos ,27(11):114326, 2017.J. Sieber and B. Krauskopf. Bifurcation analysis of an inverted pendulum with delayedfeedback control near a triple-zero eigenvalue singularity.
Nonlinearity , 17(1):85–103,2004a.J. Sieber and B. Krauskopf. Complex balancing motions of an inverted pendulum subjectto delayed feedback control.
Physica D , 197(3–4):332–345, 2004b.J. Sieber and B. Krauskopf. Extending the permissible control loop latency for thecontrolled inverted pendulum.
Dynamical Systems , 20(1):189–199, 2005.J. Sieber and R. Szalai. Characteristic matrices for linear periodic delay differentialequations.
SIAM Journal on Applied Dynamical Systems , 10(1):129–147, 2011.J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, and D. Roose. DDE-BIFTOOLmanual — bifurcation analysis of delay differential equations. arXiv preprintarXiv:1406.7144 , 2015. download: https://sourceforge.net/projects/ddebiftool . 49. Simonnet, H. A. Dijkstra, and M. Ghil. Bifurcation analysis of ocean, atmosphere,and climate models. In
Handbook of Numerical Analysis , volume 14, pages 187–229.Elsevier, 2009.A. L. Skubachevskii and H.-O. Walther. On the Floquet multipliers of periodic solutionsto nonlinear functional differential equations.
Journal of Dynamics and DifferentialEquations , 18:257–355, 2006.G. Stépán.
Retarded Dynamical Systems: Stability and Characteristic Functions . Long-man Scientific and Technical, 1989.E. Stumpf. The existence and 𝐶 -smoothness of local center-unstable manifolds for dif-ferential equations with state-dependent delay. Rostocker Mathematisches Kolloquium ,66:3–44, 2011.M. J. Suarez and P. S. Schopf. A delayed action oscillator for ENSO.
Journal of theAtmospheric Sciences , 45(21):3283–3287, 1988.R. Szalai, G. Stépán, and S.J. Hogan. Continuation of bifurcations in periodic delay differ-ential equations using characteristic matrices.
SIAM Journal on Scientific Computing ,28(4):1301–1317, 2006.E. Tziperman, L. Stone, M. A. Cane, and H. Jarosh. El Niño chaos: overlapping ofresonances between the seasonal cycle and the Pacific ocean-atmosphere oscillator.
Science , 264:72–74, 1994.E. Tziperman, M. A. Cane, S. E. Zebiak, Y. Xue, and B. Blumenthal. Locking of ElNiño’s peak time to the end of the calendar year in the delayed oscillator picture ofENSO.
Journal of Climate , 11(9):2191–2199, 1998.B. Wage. Normal form computations for delay differential equations in DDE-BIFTOOL.
Masters thesis, Universiteit Utrecht , 2014.H.-O. Walther. Stable periodic motion of a system with state-dependent delay.
Differentialand Integral Equations , 15:923–944, 2002.H.-O. Walther. The solution manifold and 𝐶 -smoothness for differential equations withstate-dependent delay. Journal of Differential Equations , 195(1):46–65, 2003.S. Yanchuk, S. Ruschel, J. Sieber, and M. Wolfrum. Temporal dissipative solitons intime-delay feedback systems.
Physical Review Letters , 123(5):053901, 2019.I. Zaliapin and M. Ghil. A delay differential model of ENSO variability — Part 2:Phase locking, multiple solutions and dynamics of extrema.
Nonlinear Processes inGeophysics , 17(2):123–135, 2010.H. Zelle, G. Appeldoorn, G. Burgers, and G. J. van Oldenborgh. The relationship betweensea surface temperature and thermocline depth in the eastern equatorial Pacific.