Manual
User Manual:
Open the PDF directly: View PDF .
Page Count: 44
Download | |
Open PDF In Browser | View PDF |
DDE-BIFTOOL v. 3.1.1 Manual — Bifurcation analysis of delay differential equations J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, D. Roose April 5, 2017 Keywords nonlinear dynamics, delay-differential equations, stability analysis, periodic solutions, collocation methods, numerical bifurcation analysis, state-dependent delay. 1 Citation, license, and obtaining the package 1 2 Version history 2.1 Changes from 3.1 to 3.1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Changes from 3.0 to 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Changes from 2.03 to 3.0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 2 3 3 Capabilities and related reading and software 3 4 Structure of DDE-BIFTOOL 5 5 Delay differential equations 5.1 Equations with constant delays . . . . . 5.1.1 Steady states . . . . . . . . . . . . 5.1.2 Periodic orbits . . . . . . . . . . . 5.1.3 Connecting orbits . . . . . . . . . 5.2 Equations with state-dependent delays 6 . . . . . . . . . . . . . . . . . . . . 6 6 6 7 7 7 System definition 6.1 Change from DDE-BIFTOOL v. 2.03 to v. 3.0 . . . . . . . . . . . . . . . . . . . . 6.2 Equations with constant delays . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Right-hand side — sys_rhs . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.2 Delays — sys_tau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.3 Jacobians of right-hand side — sys_deri (optional, but recommended) 6.3 Equations with state-dependent delays . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Right-hand side — sys_rhs . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Delays — sys_tau and sys_ntau . . . . . . . . . . . . . . . . . . . . . . . 6.3.3 Jacobian of right-hand side — sys_deri (optional, but recommended) . 6.3.4 Jacobians of delays — sys_dtau (optional, but recommended) . . . . . . 6.4 Extra conditions — sys_cond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Collecting user functions into a structure — call set_funcs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 9 9 9 10 10 11 12 12 13 13 14 14 . . . . . i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Data structures 7.1 Problem definition (functions) structure 7.2 Point structures . . . . . . . . . . . . . . 7.3 Stability structures . . . . . . . . . . . . 7.4 Method parameters . . . . . . . . . . . . 7.5 Branch structures . . . . . . . . . . . . . 7.6 Scalar measure structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 15 16 18 18 19 21 8 Point manipulation 21 9 Branch manipulation 24 10 Numerical methods 10.1 Determining systems . . . . . . . . 10.2 Extra conditions . . . . . . . . . . . 10.3 Continuation . . . . . . . . . . . . . 10.4 Roots of the characteristic equation 10.5 Floquet multipliers . . . . . . . . . . . . . . 27 27 29 30 31 33 11 Concluding comments 11.1 Existing extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 34 A Jacobians of tutorial examples neuron and sd_demo 38 B Octave compatibility considerations 42 . . . . . . . . . . . . . . . . . . . . ii . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. Citation, license, and obtaining the package DDE-BIFTOOL was started by Koen Engelborghs as part of his PhD at the Computer Science Department of the K.U.Leuven under supervision of Prof. Dirk Roose. Scientific publications for which the package DDE-BIFTOOL has been used shall mention usage of the package DDE-BIFTOOL, and shall cite the following publications to ensure proper attribution and reproducibility: • K. Engelborghs, T. Luzyanina, and D. Roose. Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL, ACM Trans. Math. Softw. 28 (1), pp. 1-21, 2002. • (this manual) J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, D. Roose . DDE-BIFTOOL v. 3.1.1 Manual — Bifurcation analysis of delay differential equations, http://arxiv.org/abs/ 1406.7144. The implementation of normal forms for equilibria is based on • M. M. Bosschaert, B. Wage and Y. Kuznetsov: Description of the extension ddebiftool nmfm http://ddebiftool.sourceforge.net/nmfm_extension_description.pdf, 2015. • B. Wage: Normal form computations for Delay Differential Equations in DDE-BIFTOOL. Master Thesis, Utrecht University (NL), supervised by Y.A. Kuznetsov (http://dspace.library.uu. nl/handle/1874/296912, 2014. • M. M. Bosschaert: Switching from codimension 2 bifurcations of equilibria in delay differential equations. Master Thesis, Utrecht University (NL), supervised by Y.A. Kuznetsov, http: //dspace.library.uu.nl/handle/1874/334792, 2016. All versions of this manual from v. 3.0 onward are available at http://arxiv.org/abs/1406.7144. Citation License The following terms cover the use of the software package DDE-BIFTOOL: BSD 2-Clause license Copyright (c) 2017, K.U. Leuven, Department of Computer Science, K. Engelborghs, T. Luzyanina, G. Samaey. D. Roose, K. Verheyden, J. Sieber, B. Wage, D. Pieroux All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 1 Download Upon acceptance of the above terms, one can obtain the package DDE-BIFTOOL (version 3.1.1) from https://sourceforge.net/projects/ddebiftool. Versions up to 3.0 and resources on theoretical background continue to be available (under a different license) on http://twr.cs.kuleuven.be/research/software/delay/ddebiftool.shtml. 2. Version history 2.1. Changes from 3.1 to 3.1.1 • Small bug fix: when changing focus on plot window during br_contn, the online plotting no longer follows focus. This keeps the online bifurcation diagram in the same window. An optional named argument 'plotaxis' has been added to br_contn to explicitly set the plot axes. • Added support for rotations (phase oscillators). Periodicity is enforced only up to multiples of 2π. Demo ../demos/phase_oscillator/html/phase_oscillator.html shows how one can track rotations. Demo was contributed by Azamat Yeldesbay. • Demos showing detection and computation of Bodganov-Takens bifurcation in ../demos/Holling-Tanner/html/HollingTanner_demo.html and cusp in ../demos/cusp/html/cusp_demo.html. Contributed by M. M. Boschaert and Y. Kuznetsov. • A description of the mathematical formulas behind the normal form computations is now in nmfm_extension_description.pdf, by M.Bossschaert, B. Wage, Y. Kuznetsov. 2.2. Changes from 3.0 to 3.1 D. Roose has permitted to change the license to a Sourceforge-compliant BSD License. Thus, code and newest releases from version 3.1 onward are now available from https://sourceforge.net/projects/ddebiftool. Older versions will continue to be available from http://twr.cs.kuleuven.be/research/software/delay/ddebiftool. shtml. Change of License and move to Sourceforge The new functionality is only applicable for equations with constant delay. Normal form coefficients can be computed through the extension ddebiftool extra nmfm. This extension is included in the standard DDE-BIFTOOL archive, but the additional functions are kept in a separate folder. The following bifurcations are currently supported. New feature: Normal form computation for bifurcations of equilibria • Hopf bifurcation (coefficient L1 determining criticality), • generalized Hopf (Bautin) bifurcation (of codimension two, typically encountered along Hopf curves) • Zero-Hopf interaction (Gavrilov-Guckenheimer bifurcation, of codimension two, typically encountered along Hopf curves) 2 • Hopf-Hopf interaction (of codimension two, typically encountered along Hopf curves) The extension comes with a demo nmfm demo. The demos neuron, minimal demo and Mackey-Glass illustrate the new functionality, too. Background theory is given in [46]. 2.3. Changes from 2.03 to 3.0 New features • Continuation of local periodic orbit bifurcations for systems with constant or state-dependent delay is now supported through the extension ddebiftool extra psol. This extension is included in the standard DDE-BIFTOOL archive, but the additional functions are kept in a separate folder. • User-defined functions specifying the right-hand side and delays (such as sys_rhs and sys_tau) can have arbitrary names. These user functions (with arbitrary names) get collected in a structure funcs, which gets then passed on to the DDE-BIFTOOL routines. This interface is similar to other functions acting on Matlab functions such as fzero or ode45. It enables users to add extensions such as ddebiftool extra psol without changing the core routines. • State-dependent delays can now have arbitrary levels of nesting (for example, periodic orbits of ẋ (t) = µ − x (t − x (t − x (t − x (t)))) and their bifurcations can be tracked). • Vectorization Continuation of periodic orbits and their bifurcations benefits (moderately) from vectorization of the user-defined functions. • Utilities Recurring tasks (such as branching off at bifurcations, defining initial pieces of branches, or extracting the number of unstable eigenvalues) can now be performed more conveniently with some auxiliary functions provided in a separate folder ddebiftool utilities. See ../demos/neuron/html/demo1_simple.html for a demonstration. • Bugs fixed Some bugs and problems have been fixed in the implementation of the heuristics applied to choose the stepsize in the computation of eigenvalues of equilibria [45]. • Continuation of relative equilibiria and relative periodic orbits and their local bifurcations for systems with constant delay and rotational symmetry (saddle-node bifurcation, Hopf bifurcation, period-doubling, and torus bifurcation) is now supported through the extension ddebiftool extra rotsym. Extensions come with demos and separate documentation. Versions from 3.0 onward have a different the user interface for many DDE-BIFTOOL functions than versions up to 2.03. They add one additional input argument funcs (this new argument comes always first). Since DDE-BIFTOOL v. 3.0 changed the user interface, scripts written for DDE-BIFTOOL v. 2.0x or earlier will not work with versions later than 3.0. For this reason version 2.03 will continue to be available. Users of both versions should ensure that only one version is in the Matlab path at any time to avoid naming conflicts. Change of user interface 3. Capabilities and related reading and software DDE-BIFTOOL consists of a set of routines running in Matlab1 [29] or octave2 , both widely used environments for scientific computing. The aim of the package is to provide a tool for numerical 1 http://www.mathworks.com 2 http://www.gnu.org/software/octave 3 bifurcation analysis of steady state solutions and periodic solutions of differential equations with constant delays (called DDEs) or state-dependent delays (here called sd-DDEs). It also allows users to compute homoclinic and heteroclinic orbits in DDEs (with constant delays). DDE-BIFTOOL can perform the following computations: • continuation of steady state solutions (typically in a single parameter); • approximation of the rightmost, stability-determining roots of the characteristic equation which can further be corrected using a Newton iteration; • continuation of steady state folds and Hopf bifurcations (typically in two system parameters); • continuation of periodic orbits using orthogonal collocation with adaptive mesh selection (starting from a previously computed Hopf point or an initial guess of a periodic solution profile); • approximation of the largest stability-determining Floquet multipliers of periodic orbits; • branching onto the secondary branch of periodic solutions at a period doubling bifurcation or a branch point; • continuation of folds, period doublings and torus bifurcations (typically in two system parameters) using the extension ddebiftool extra psol; • computation of normal form coefficients for Hopf bifurcations and codimension-two bifurcations along Hopf bifurcation curves (typically in two system parameters) using the extension ddebiftool extra nmfm; • continuation of connecting orbits (using the appropriate number of parameters). All computations can be performed for problems with an arbitrary number of discrete delays. These delays can be either parameters or functions of the state. (The only exception are computations of connecting orbits, which support only problems with delays as parameters at the moment.) A practical difference to AUTO or MatCont is that the package does not detect bifurcations automatically because the computation of eigenvalues or Floquet multipliers may require more computational effort than the computation of the equilibria or periodic orbits (for example, if the system dimension is small but one delay is large). Instead the evolution of the eigenvalues can be computed along solution branches in a separate step if required. This allows the user to detect and identify bifurcations. Capabilities This manual documents version 3.1.1. Earlier versions of the manual for earlier versions of DDE-BIFTOOL continue to be available at the web addresses About this manual — related reading versions ≥ 3.0 versions ≤ 2.03 http://arxiv.org/abs/1406.7144 http://twr.cs.kuleuven.be/research/software/delay/ddebiftool.shtml. For readers who intend to analyse only systems with constant delays, the parts of the manual related to systems with state-dependent delays can be skipped (sections 5.2, 6.3). In the rest of this manual we assume the reader is familiar with the notion of a delay differential equation and with the basic concepts of bifurcation analysis for ordinary differential equations. The theory on delay differential equations and a large number of examples are described in several books. Most notably the early [4, 13, 14, 25, 32] and the more recent [2, 30, 26, 10, 31]. Several excellent books contain introductions to dynamical systems and bifurcation theory of ordinary differential equations, see, e.g., [1, 6, 23, 33, 40]. The tutorial demos demo1 and sd demo, providing a step-by-step walk-through for the typical working mode with DDE-BIFTOOL are included as separate html files, published directly from the comments in the demo code. See ../demos/index.html for links to all demos, many of which are extensively commented. Turorial demos 4 USER Layer 3 : BRANCH MANIPULATION Layer 2 : POINT MANIPULATION Layer 1 : NUMERICAL METHODS Layer 0 : SYSTEM DEFINITION Figure 1: The structure of DDE-BIFTOOL. Arrows indicate the calling (−) or writing (·−) of routines in a certain layer. Related software A large number of packages exist for numerical continuation and bifurcation analysis of systems of ordinary differential equations. Currently maintained packages are AUTO url: http://sourceforge.net/projects/auto-07p MatCont url: http://sourceforge.net/projects/matcont/ Coco url: http://sourceforge.net/projects/cocotools using FORTRAN or C [12, 11], for Matlab [9, 22], and for Matlab [8]. For delay differential equations the package knut url: http://gitorious.org/knut using C++ (formerly PDDECONT) is available as a stand-alone package (written in C++, but with a user interface requiring no programming). This package was developed in parallel with DDE-BIFTOOL but independently by R. Szalai [44, 38]. For simulation (time integration) of delay differential equations the reader is, e.g., referred to the packages ARCHI, DKLAG6, XPPAUT, DDVERK, RADAR and dde23, see [37, 7, 21, 20, 41, 24]. Of these, only XPPAUT has a graphical interface (and allows limited stability analysis of steady state solutions of DDEs along the lines of [35]). TRACE-DDE is a Matlab tool (with graphical interface) for linear stability analysis of linear constant-coefficient DDEs [5]. 4. Structure of DDE-BIFTOOL The structure of the package is depicted in figure 1. It consists of four layers. Layer 0 contains the system definition and consists of routines which allow to evaluate the right hand side f and its derivatives, state-dependent delays and their derivatives and to set or get the parameters and the constant delays. It should be provided by the user and is explained in more detail in section 6. All user-provided functions are collected in a single structure (called funcs in this manual), and are passed on by the user as arguments to layer-3 or layer-2 functions. Note that this is a change in user interface between version 2.03 and version 3.0! Layer 1 forms the numerical core of the package and is (normally) not directly accessed by the user. The numerical methods used are explained briefly in section 10, more details can be found in the papers [35, 18, 17, 16, 19, 34, 39] and in [15]. Its functionality is hidden by and used through layers 2 and 3. Layer 2 contains routines to manipulate individual points. Names of routines in this layer start with ”p ”. A point has one of the following five types. It can be a steady state point (abbreviated 5 'stst'), steady state Hopf (abbreviated 'hopf') or fold (abbreviated 'fold') bifurcation point, a periodic solution point (abbreviated 'psol') or a connecting orbit point (abbreviated 'hcli'). Furthermore, a point can contain additional information concerning its stability. Routines are provided to compute individual points, to compute and plot their stability and to convert points from one type to another. Layer 3 contains routines to manipulate branches. Names of routines in this layer start with ”br ”. A branch is structure containing an array of (at least two) points, three sets of method parameters and specifications concerning the free parameters. The 'point' field of a branch contains an array of points of the same type ordered along the branch. The 'method' field contains parameters of the computation of individual points, the continuation strategy and the computation of stability. The 'parameter' field contains specification of the free parameters (which are allowed to vary along the branch), parameter bounds and maximal step sizes. Routines are provided to extend a given branch (that is, to compute extra points using continuation), to (re)compute stability along the branch and to visualize the branch and/or its stability. Layers 2 and 3 require specific data structures, explained in section 7, to represent points, stability information, branches, to pass method parameters and to specify plotting information. Usage of these layers is demonstrated through a step-by-step analysis of the demo systems neuron, sd demo and hom demo (see ../demos/index.html). Descriptions of input/output parameters and functionality of all routines in layers 2 and 3 are given in sections 8 respectively 9. 5. Delay differential equations This section introduces the mathematical notation that we refer to in this manual to describe the problems solved by DDE-BIFTOOL. 5.1. Equations with constant delays Consider the system of delay differential equations with constant delays (DDEs), d x (t) = f ( x (t), x (t − τ1 ), . . . , x (t − τm ), η ), dt (1) where x (t) ∈ Rn , f : Rn(m+1) × R p → Rn is a nonlinear smooth function depending on a number of parameters η ∈ R p , and delays τi > 0, i = 1, . . . , m. Call τ the maximal delay, τ = max τi . i =1,...,m The linearization of (1) around a solution x ∗ (t) is the variational equation, given by, m d y(t) = A0 (t)y(t) + ∑ Ai (t)y(t − τi ), dt i =1 (2) where, using f ≡ f ( x0 , x1 , . . . , x m , η ), Ai ( t ) = ∂f ∗ ( x (t), x ∗ (t − τ1 ), . . . , x ∗ (t − τm ), η ), i = 0, . . . , m. ∂xi 5.1.1. Steady states If x ∗ (t) corresponds to a steady state solution, x ∗ (t) ≡ x ∗ ∈ Rn , with f ( x ∗ , x ∗ , . . . , x ∗ , η ) = 0, 6 (3) then the matrices Ai (t) are constant, Ai (t) ≡ Ai , and the corresponding variational equation (2) leads to a characteristic equation. Define the n × n-dimensional matrix ∆ as m ∆(λ) = λI − A0 − ∑ Ai e−λτi . (4) i =1 Then the characteristic equation reads, det(∆(λ)) = 0. (5) Equation (5) has an infinite number of roots λ ∈ C which determine the stability of the steady state solution x ∗ . The steady state solution is (asymptotically) stable provided all roots of the characteristic equation (5) have negative real part; it is unstable if there exists a root with positive real part. It is known that the number of roots in any right half plane Re(λ) > γ, γ ∈ R is finite, hence, the stability is always determined by a finite number of roots. Bifurcations occur whenever roots move through the imaginary axis as one or more parameters are changed. Generically a fold bifurcation (or turning point) occurs when the root is real (that is, equal to zero) and a Hopf bifurcation occurs when a pair of complex conjugate roots crosses the imaginary axis. 5.1.2. Periodic orbits A periodic solution x ∗ (t) is a solution which repeats itself after a finite time, that is, x ∗ (t + T ) = x ∗ (t), for all t. Here T > 0 is the period. The stability around the periodic solution is determined by the time integration operator S( T, 0) which integrates the variational equation (2) around x ∗ (t) from time t = 0 over the period. This operator is called the monodromy operator and its (infinite number of) eigenvalues, which are independent of the starting moment t = 0, are called the Floquet multipliers. Furthermore, if S( T, 0)k is compact for k > τ/T. Thus, there are at most finitely many Floquet multipliers outside of any ball around the origin of the complex plane. For autonomous systems there is always a trivial Floquet multiplier at unity, corresponding to a perturbation along the time derivative of the periodic solution. The periodic solution is exponentially stable provided all multipliers (except the trivial one) have modulus smaller than unity, it is exponenially unstable if there exists a multiplier with modulus larger than unity. 5.1.3. Connecting orbits We call a solution x ∗ (t) of (1) at η = η ∗ a connecting orbit if the limits lim x ∗ (t) = x − , lim x ∗ (t) = x + , t→−∞ t→+∞ (6) exist. For continuous f , x − and x + are steady state solutions. If x − = x + , the orbit is called homoclinic, otherwise it is heteroclinic. 5.2. Equations with state-dependent delays Consider the system of delay differential equations with state-dependent delays (sd-DDEs), d x ( t ) = f ( x , x , . . . , x , η ), m 0 1 dt x j = x (t − τj ( x0 , . . . , x j−1 , η ) (τ0 = 0, j = 1, . . . , m), 7 (7) where x (t) ∈ Rn , and f :Rn(m+1) × R p → Rn τj :Rn j × R p → [0, ∞) are smooth functions depending of their arguments. The right-hand side f depends on m + 1 states x j = x (t − τj ) ∈ Rn (j = 0, . . . , m) and p parameters η ∈ R p . The jth delay function τj depends on all previously defined j − 1 states x (t − τi ) ∈ Rn (i = 0, . . . , j − 1) and p parameters η ∈ R p . This definition permits the user to formulate sd-DDEs with arbitrary levels of nesting in their function arguments. The linearization around a solution ( x ∗ (t), η ∗ ) of (7) (the variational equation) with respect to x is given by (see [27], we are using the notation x0∗ = x ∗ (t), τj∗ (t) = τj ( x0∗ , . . . , x ∗j−1 , η ∗ ) and x ∗j = x ∗ (t − τj∗ (t)) for j ≥ 1) d y(t) = dt m ∑ A j (t)Yk j =0 Y0 = y(t) (8) Yj = y(t − τj∗ (t)) − ( x ∗ )0 (t − τj (t)) j −1 ∑ Bk,j (t)Yk , (j = 1 . . . m) k =0 0 where ( x ∗ ) (t) = dx ∗ (t)/dt, and ∂f ∗ ∗ ∗ ( x , x , . . . , xm , η ) ∈ Rn×n , (j = 0, . . . , m), ∂x j 0 1 ∂τj Bj,k (t) = k ( x0∗ , x1∗ , . . . , x ∗j−1 , η ∗ ) ∈ R1×n , (k = 0, . . . , j − 1, j = 1, . . . , m). ∂x A j (t) = (9) ∗ ≡ x ∗ ∈ Rn , If ( x ∗ (t), τ̃ ∗ (t)) corresponds to a steady state solution, then x ∗ (t) = x0∗ = . . . = xm ∗ ∗ ∗ ∗ and τj (t) ≡ τj ( x , . . . , x , η ) for all j ≥ 1, with f (x∗ , x∗ , . . . , x∗ , η ∗ ) = 0 then the matrices Ai (t) are constant, Ai (t) ≡ Ai , and the vectors Bi,j (t) consist of zero elements only. In this case, the corresponding variational equation (8) is a constant delay differential equation and it leads to the characteristic equation (5), i.e. a characteristic equation with constant delays. Hence the stability analysis of a steady state solution of (7) is similar to the stability analysis of (1). Note that the right-hand side f , when considered as a functional mapping a history segment into Rn is not locally Lipschitz continuous. This creates technical difficulties when considering an sd-DDE of type (7) as an infinite-dimensional system, because the solution does not depend smoothly on the initial condition (see Hartung et al [27] for a detailed review). However, periodic boundary-value problems for (7) can be reduced to finite-dimensional systems of algebraic equations that are as smooth as the coefficient functions f and τj [43]. This implies that all periodic orbits and their bifurcations and stability as computed by DDE-BIFTOOL behave as expected. In particular, branching off at Hopf bifurcations and period doubling works in the same way as for constant delays (the proof for the Hopf bifurcation in sd-DDEs is also given in [43]). Moreover, Mallet-Paret and Nussbaum [36] proved that the stability of the linear variational equation (8) indeed reflects the local stability of the solution ( x ∗ (t), τ̃ ∗ (t)) of (7). For details on the relevant theory and numerical bifurcation analysis of differential equations with state-dependent delay see [34, 27] and the references therein. 8 6. System definition 6.1. Change from DDE-BIFTOOL v. 2.03 to v. 3.0 Note that in DDE-BIFTOOL 3.1.1 all user-provided functions can be arbitrary function handles, collected into a structure using the function set_funcs, explained in section 6.5. The only typical mandatory functions for the user to provide are the right-hand side ('sys_rhs') and the function returning the delay indices ('sys_tau'). The names of the user functions can be arbitrary, and user functions can be anonymous. This is a change from previous versions. See the tutorials in ../demos/index.html for examples of usage, and the function description in section 6.5 for details. 6.2. Equations with constant delays As an illustrative example we will use the following system of delay differential equations, taken from [42], x˙1 (t) = −κx1 (t) + β tanh( x1 (t − τs )) + a12 tanh( x2 (t − τ2 )) (10) x˙2 (t) = −κx2 (t) + β tanh( x2 (t − τs )) + a21 tanh( x1 (t − τ1 )). This system models two coupled neurons with time delayed connections. It has two components (x1 and x2 ), three delays (τ1 , τ2 and τs ), and four parameters (κ, β, a12 and a21 ). The demo neuron (see ../demos/index.html) walks through the bifurcation analysis of system (10) step by step to demonstrate the working pattern for DDE-BIFTOOL. To define a system, the user should provide the following Matlab functions, given in the following paragraphs for system (10). 6.2.1. Right-hand side — sys_rhs The right-hand side is a function of two arguments. For our example (10), this would have the form (giving the right-hand side the name neuron_sys_rhs) neuron_sys_rhs = @ ( xx , par )[... - par (1)* xx (1 ,1)+ par (2)* tanh ( xx (1 ,4))+ par (3)* tanh ( xx (2 ,3));... - par (1)* xx (2 ,1)+ par (2)* tanh ( xx (2 ,4))+ par (4)* tanh ( xx (1 ,2))]; % par =[\ kappa ,\ beta , a_ {12} , a_ {21} ,\ tau_1 ,\ tau_2 , \ tau_s ] Listing 1: Definition for right-hand side of (10) as a variable. Meaning of the arguments of the right-hand side function: • xx ∈ Rn×(m+1) contains the state variable(s) at the present and in the past, • par ∈ R1× p contains the parameters, par = η. The delays τi (i = 1 . . . , m) are considered to be part of the parameters (τi = η j(i) , i = 1, . . . , m). This is natural since the stability of steady solutions and the position and stability of periodic solutions depend on the values of the delays. Furthermore delays can occur both as a ‘physical’ parameter and as delay, as in ẋ = τx (t − τ ). From these inputs the right hand side f is evaluated at time t. Notice that the parameters have a specific order in par indicated in the comment line. An alternative (vectorized) form would be neuron_sys_rhs = @ ( xx , p )[... -p (1)* xx (1 ,1 ,:)+ p (2)* tanh ( xx (1 ,4 ,:))+ p (3)* tanh ( xx (2 ,3 ,:));.... -p (1)* xx (2 ,1 ,:)+ p (2)* tanh ( xx (2 ,4 ,:))+ p (4)* tanh ( xx (1 ,2 ,:))]; Listing 2: Alternative definition of the right-hand side of (10), vectorized for speed-up of periodic orbit computations. 9 Note the additional colon in argument xx and compare to Listing 1. The form shown in Listing 2 can be called in many points along a mesh simultaneously, speeding up the computations during analysis of periodic orbits. 6.2.2. Delays — sys_tau For constant delays another function is required which returns the position of the delays in the parameter list. For our example, this is neuron_tau = @ ()[5 6 7]; This function has no arguments for constant delays, and returns a row vector of indices into the parameter vector. 6.2.3. Jacobians of right-hand side — sys_deri (optional, but recommended) Several derivatives of the right hand side function f need to be evaluated during bifurcation analysis. By default, DDE-BIFTOOL uses a finite-difference approximation, implemented in df deriv.m. For speed-up or in case of convergence difficulties the user may provide the Jacobians of the right-hand side analytically as a separate function. Its header is of the format function J = sys_deri ( xx , par , nx , np , v ) Arguments: • xx ∈ Rn×(m+1) contains the state variable(s) at the present and in the past (as for the right-hand side); • par ∈ R1× p contains the parameters, par = η (as for the right-hand side); • nx (empty, one integer or two integers) index (indices) of xx with respect to which the righthand side is to be differentiated • np (empty or integer) whether right-hand side is to be differentiated with respect to parameters • v (empty or Cn ) for mixed derivatives with respect to xx, only the product of the mixed derivative with v is needed. The result J is a matrix of partial derivatives of f which depends on the type of derivative requested via nx and np multiplied with v (when nonempty), see table 1. J is defined as follows. Initialize J with f . If nx is nonempty take the derivative of J with respect to those arguments listed in nx’s entries. Each entry of nx is a number between 0 and m based on f ≡ f ( x0 , x1 , . . . , x m , η ). E.g., if nx has only one element take the derivative with respect to xnx(1) . If it has two elements, take, of the result, the derivative with respect to xnx(2) and so on. Similarly, if np is nonempty take, of the resulting J, the derivative with respect to ηnp(i) where i ranges over all the elements of np, 1 ≤ i ≤ p. Finally, if v is not an empty vector multiply the result with v. The latter is used to prevent J from being a tensor if two derivatives with respect to state variables are taken (when nx contains two elements). Not all possible combinations of these derivatives have to be provided. In the current version, nx has at most two elements and np at most one. The possibilities are further restricted as listed in table 1. In the last row of table 1 the elements of J are given by, ! n ∂ ∂ ∂ fi J i,j = Anx(1) v = ∑ nx(1) vk , ∂xnx(2) i,j ∂xnx(2) ∂x j with Al as defined in (3). 10 k =1 k length(nx) length(np) v J 1 0 empty ∂f = Anx(1) ∈ Rn×n nx(1) ∂x 0 1 empty ∂f ∈ Rn ×1 ∂ηnp(1) 1 1 empty 2 0 ∈ Cn ×1 ∂2 f ∈ Rn × n ∂xnx(1) ∂ηnp(1) ∂ Anx(1) v ∈ Cn×n ∂xnx(2) Table 1: Results of the function sys deri depending on its input parameters nx, np and v using f ≡ f ( x 0 , x 1 , . . . , x m , η ). The resulting routine is quite long, even for the small system (10); see Listing 6 in Appendix A for a printout of the function body. Furthermore, implementing so many derivatives is an activity prone to a number of typing mistakes. Hence a default routine df_deriv is available which implements finite difference formulas to approximate the requested derivatives (using several calls to the right-hand side. It is, however, recommended to provide at least the first order derivatives with respect to the state variables using analytical formulas. These derivatives occur in the determining systems for fold and Hopf bifurcations and for connecting orbits, and in the computation of characteristic roots and Floquet multipliers. All other derivatives are only necessary in the Jacobians of the respective Newton procedures and thus influence only the convergence speed. 6.3. Equations with state-dependent delays DDE-BIFTOOL also permits the delays to depend on parameters and the state. If at least one delay is state-dependent then the format and semantics of the function specifying the delays, sys_tau, is different from the format used for constant delays in section 6.2.2 (it now provides the values of the delays). Note that for a system with only constant delays we recommend the use of the system definitions as described in section 6.2 to reduce the computational effort. As an illustrative example we will use the following system of delay differential equations, d x (t) dt 1 d x2 ( t ) dt d x3 ( t ) dt d x (t) dt 4 d x5 ( t ) dt 1 (1 − p2 x1 (t) x1 (t − τ3 ) x3 (t − τ3 ) + p3 x1 (t − τ1 ) x2 (t − τ2 )) , p1 + x2 ( t ) p4 x1 ( t ) = + p5 tanh( x2 (t − τ5 )) − 1, p1 + x2 ( t ) = = p6 ( x2 (t) − x3 (t)) − p7 ( x1 (t − τ6 ) − x2 (t − τ4 ))e− p8 τ5 , = x1 (t − τ4 )e− p1 τ5 − 0.1, = 3( x1 (t − τ2 ) − x5 (t)) − p9 , 11 (11) function f = sd_rhs ( xx , par ) f (1 ,1)=(1/( par (1)+ xx (2 ,1)))*(1 - par (2)* xx (1 ,1)* xx (1 ,4)* xx (3 ,4)+... par (3)* xx (1 ,2)* xx (2 ,3)); f (2 ,1)= par (4)* xx (1 ,1)/( par (1)+ xx (2 ,1))+ par (5)* tanh ( xx (2 ,6)) -1; f (3 ,1)= par (6)*( xx (2 ,1) - xx (3 ,1)) - par (7)*( xx (1 ,7) -... xx (2 ,5))* exp ( - par (8)* xx (4 ,1)); f (4 ,1)= xx (1 ,5)* exp ( - par (1)* xx (4 ,1)) -0.1; f (5 ,1)=3*( xx (1 ,3) - xx (5 ,1)) - par (9); end Listing 3: Listing of right-hand side 'sys_rhs' function for (11) function for (11), here called sd_rhs. where τ1 , τ2 are constant delays, τ3 = 2 + p5 τ1 x2 (t) x2 (t − τ1 ), 1 τ4 = 1 − , 1 + x1 (t) x2 (t − τ2 ) τ5 = x4 (t), τ6 = x5 (t). This system has five components ( x1 , . . . , x5 ), six delays (τ1 , . . . , τ6 ) and eleven parameters ( p1 , . . . , p11 ), where p10 = τ1 and p11 = τ2 . A step-by-step tutorial for analysis of sd-DDEs is given in demo sd demo (see ../demos/index.html) of this system using (11). To define a system with state-dependent delays, the user should provide the following Matlab functions, given in the following sections for system (11). 6.3.1. Right-hand side — sys_rhs The definition and functionality of this routine is equivalent to the one described in section 6.2.1. Notice that the argument xx contains the state variable(s) at the present and in the past, xx = [ x (t) x (t − τ1 ) . . . . . . x (t − τm )]. Possible constant delays (τ1 and τ2 in example (11)) are also considered to be part of the parameters. See Listing 3 for the right-hand side to be provided for example (11). 6.3.2. Delays — sys_tau and sys_ntau The format and semantics of the routines specifying the delays differ from the one described in section 6.2.2. The user has to provide two functions: function ntau = sys_ntau () function tau = sys_tau ( ind , xx , par ) The function sys_ntau has no arguments and returns the number of (constant and state-dependent) delays. For the example (11), this could be the anonymous function @()6;. The function sys_tau has the three arguments: • ind (integer ≥ 1) indicates, which delay is to be returned; • xx (n × ind-matrix) is the state: xx(:,1) = x (t), xx(:,k) = x (t − τk−1 ) for k = 2 . . . ind; • par (row vector) is the vector of system parameters 12 function tau = sd_tau ( ind , xx , par ) if ind ==1 tau = par (10); elseif ind ==2 tau = par (11); elseif ind ==3 tau =2+ par (5)* par (10)* xx (2 ,1)* xx (2 ,2); elseif ind ==4 tau =1 -1/(1+ xx (2 ,3)* xx (1 ,1)); elseif ind ==5 tau = xx (4 ,1); elseif ind ==6 tau = xx (5 ,1); end end Listing 4: Listing of 'sys_tau' function for (11), here called sd_tau. The output is the value of τind , the ind’th delay. DDE-BIFTOOL calls sys_tau m times where m =sys_ntau(). In the first call xx is the n × 1 vector, x (t) such that τ1 may depend on x (t) and par. After the first call to sys_tau, DDE-BIFTOOL computes x (t − τ1 ). In the second call to sys_tau xx is a n × 2 matrix, consisting of [ x (t), x (t − τ1 )] such that the delay τ2 may depend on x (t) , x (t − τ1 ) and par, etc. In this way, the user can define state-dependent point delays with arbitrary levels of nesting. The delay function for the example (11) is given in Listing 4. Note: order of delays The order of the delays requested in sys_tau corresponds to the order in which they appear in xx as passed to the functions sys_rhs and sys_deri. Note: difference to section 6.2.2 When calling sys_tau for a constant delay, the value of the delay is returned. This is in contrast with the definition of sys_tau in section 6.2.2, where the position in the parameter list is returned. 6.3.3. Jacobian of right-hand side — sys_deri (optional, but recommended) The definition and functionality of this routine is equivalent to the one described in section 6.2.3. We do not present here the routine since it is quite long, see the Matlab code sd deri.m in the demo example sd demo. If the user does not provide a function for the Jacobians the finite-difference approximation (df deriv.m) will be used by default. However, as for constant delays, it is recommended to provide at least the first order derivatives with respect to the state variables using analytical formulas. 6.3.4. Jacobians of delays — sys_dtau (optional, but recommended) The routine of the format function dtau = sd_dtau ( ind , xx , par , nx , np ) supplies derivatives of all delays with respect to the state and parameters. Its functionality is similar to the function sys_deri. Inputs: • ind (integer ≥ 1) the number of the delay, 13 • par (n × ind matrix) is the state: xx(:,1) = x (t), xx(:,k) = x (t − τk−1 ) for k = 2 . . . ind; • nx (empty, one integer or two integers) index (indices) of xx with respect to which the righthand side is to be differentiated; • np (empty or integer) whether right-hand side is to be differentiated with respect to parameters. The result dtau is a scalar, vector or matrix of partial derivatives of the delay with number ind, which depends on the type of derivative requested via nx and np, see table 2. The resulting routine is quite long, even for the small system (11); see Listing 7 in Appendix A for a printout of the function body. If the user does not provide a 'sys_dtau' function then the default routine df_derit will be used, which implements finite difference formulas to approximate the requested derivatives (using several calls to sys_tau), analogously to df_deriv. As in the case of sys_deri, it is recommended to provide at least the first order derivatives with respect to the state variables using analytical formulas. length(nx) length(np) dtau 1 0 ∂τind ∈ Rn ∂x nx(1) 0 1 ∂τind ∈R ∂ ηnp(1) 1 1 2 0 ∂2 τind ∈ Rn ∂xnx(1) ∂ηnp(1) ∂τind ∂ ∈ Rn × n ∂xnx(2) ∂xnx(1) Table 2: Results of the function sys_dtau depending on its input parameters nx and np (ind= 1, . . . , m). 6.4. Extra conditions — sys_cond A system routine sys_cond with a header of the type function [ res , p ]= sys_cond ( point ) can be used to add extra conditions during corrections and continuation, see section 10.2 for an explanation of arguments and outputs. 6.5. Collecting user functions into a structure — call set_funcs The user-provided functions are passed on as an additional argument to all routines of DDEBIFTOOL (similar to standard Matlab routines such as ode45). This was changed in DDE-BIFTOOL 3.0 from previous versions. The additional argument is a structure funcs containing all the handles to all user-provided functions. In order to create this structure the user is recommended to call the function set_funcs at the beginning of the script performing the bifurcation analysis: function funcs = set_funcs (...) Its argument format is in the form of name-value pairs (in arbitrary order, similar to options at the end of a call to plot). For the example (10) of a neuron, discussed in section 6.2 and in demo neuron (see ../demos/index.html), the call to set_funcs could look as follows: 14 funcs = set_funcs ( ' sys_rhs ' , neuron_sys_rhs , ' sys_tau ' ,@ ()[5 ,6 ,7] ,... ' sys_deri ' , @neuron_sys_deri ); Note that neuron_sys_rhs is a variable (a function handle pointing to an anonymous function defined as in section 6.2.1), and neuron sys deri.m is the filename in which the function providing the system derivatives are defined (see section 6.2.3). The delay function 'sys_tau' is directly specified as an anonymous function in the call to set_funcs (not needing to be defined in a separate file or as a separate variable). If one does wish to not provide analytical derivatives, one may drop the 'sys_deri' pair (then a finite-difference approximation, implemented in df_deriv, is used): funcs = set_funcs ( ' sys_rhs ' , neuron_sys_rhs , ' sys_tau ' ,@ ()[5 ,6 ,7]); For the sd-DDE example (11), the call could look as follows: funcs = set_funcs ( ' sys_rhs ' , @sd_rhs , ' sys_tau ' , @sd_tau ,... ' sys_ntau ' ,@ ()6 , ' sys_deri ' , @sd_deri , ' sys_dtau ' , @sd_dtau ); Possible names to be used in the argument sequence equal the resulting field names in funcs (see Table 3 in section 7.1 later): • 'sys_rhs': handle of the user function providing the right-hand side, described in sections 6.2.1 and 6.3.1; • 'sys_tau': handle of the user function providing the indices of the delays in the parameter vector for DDEs with constant delays, described in sections 6.2.2, or providing the values of the delays for sd-DDEs as described in section 6.3.2; • 'sys_ntau' (relevant for sd-DDEs only): handle of the user function providing the number of delays in sd-DDEs as described in section 6.3.2; • 'sys_deri': handle of the user function providing the Jacobians of right-hand side, described in sections 6.2.3 and 6.3.3; • 'sys_dtau' (relevant for sd-DDEs only): handle of the user function providing the Jacobians of delays (for sd-DDEs), described in section 6.2.3; • 'x_vectorized' (logical, default false): if the functions in 'sys_rhs', 'sys_deri' (if provided), 'sys_tau' (for sd-DDEs) and 'sys_dtau' (for sd-DDEs if provided) can be called with 3d arrays in their xx argument. Vectorization will speed up computations for periodic orbits only. An example for a necessary modification of the right-hand side to permit vectorization is given for the neuron example in Listing 2. The output funcs is a structure containing all user-provided functions and defaults for the Jacobians if they are not provided. This output is passed on as first argument to all DDE-BIFTOOL routines during bifurcation analysis. 7. Data structures In this section we describe the data structures used to define the problem, and to present individual points, stability information, branches of points, method parameters and plotting information. 7.1. Problem definition (functions) structure The user-provided functions described in Section 6 get passed on to DDE-BIFTOOL’s routines collected in a single argument funcs, a structure containing at least the fields listed in Table 3. Only fields marked with [!] are mandatory (for sd-DDEs 'sys_ntau' is mandatory, too). The user does 15 field content 'sys_rhs' (mandatory) function handle 'sys_ntau' function handle default function in file sys rhs.m if found in current working folder @()0 'sys_tau' (mandatory) function handle function in file sys tau.m if found in current working folder 'sys_cond' function handle @(p)dummy_cond, a built-in routine that adds no conditions 'sys_deri' function handle df_deriv, coming with DDE-BIFTOOL and using finite-difference approximation 'sys_dtau' function handle df_derit, coming with DDE-BIFTOOL and using finite-difference approximation 'x_vectorized' logical false ('tp_del') logical n/a (automatically determined from the number of arguments expected by funcs.sys_tau) ('sys_deri_provided') logical n/a (set automatically) ('sys_dtau_provided') logical n/a (set automatically) Table 3: Problem definition structure containing (at least) the user-provided functions. Fields in brackets should not normally be set or manually changed by the user. See also section 6.5. usually not have to set the fields of this structure manually, but calls the routine set_funcs, which returns the structure funcs to be passed on to other functions. The usage of set_funcs and the meaning of the fields of its output are described in detail in section 6.5. See also the tutorial demos neuron and sd demo (see ../demos/index.html for examples of usage. 7.2. Point structures Table 4 describes the structures used to represent a single steady state, fold, Hopf, periodic and homoclinic/heteroclinic solution point. A steady state solution is represented by the parameter values η (which contain also the constant delay values, see section 6) and x ∗ . A fold bifurcation is represented by the parameter values η, its position x ∗ and a null-vector of the characteristic matrix ∆(0). A Hopf bifurcation is represented by the parameter values η, its position x ∗ , a frequency ω and a (complex) null-vector of the characteristic matrix ∆(iω ). A periodic solution is represented by the parameter values η, the period T and a time-scaled profile x ∗ (t/T ) on a mesh in [0,1]. The mesh is an ordered collection of interval points {0 = t0 < t1 < . . . < t L = 1} and representation points ti+ j , i = 0, . . . , L − 1, j = 1, . . . , d − 1 which need to be chosen in function of the interval points as d j t i + j = t i + ( t i +1 − t i ). d d The profile is a continuous piecewise polynomial on the mesh. More specifically, it is a polynomial of degree d on each subinterval [ti , ti+1 ], i = 0, . . . , L − 1. Each of these polynomials is uniquely represented by its Warning: this assumption is not checked but needs to be fulfilled for correct results! 16 field field content content field content 'kind' 'stst' 'parameter' R1× p 'x' Rn ×1 'stability' empty or struct 'kind' 'fold' 'parameter' R1 × p 'x' Rn ×1 'v' Rn ×1 'stability' empty or struct 'kind' 'hopf' 'parameter' R1 × p 'x' Rn ×1 'v' Cn ×1 'omega' R 'stability' empty or struct (a) Steady state (b) Steady state fold (c) Steady state Hopf field content field content 'kind' 'parameter' 'mesh' 'degree' 'profile' 'period' 'stability' 'psol' 'kind' 'parameter' 'mesh' 'degree' 'profile' 'period' 'x1' 'x2' 'lambda_v' 'lambda_w' 'v' 'w' 'alpha' 'epsilon' 'hcli' R1 × p 1 ×( Ld +1) or empty [0, 1] N0 n ×( R Ld+1) R0+ empty or struct (d) Periodic orbit R1 × p or empty N0 [0, 1]1×( Ld+1) Rn×( Ld+1) R0+ Rn Rn Cs1 Cs2 Cn × s1 Cn × s2 Cs1 R (e) Connecting orbit Table 4: Point structures: Field names and corresponding content for the point structures used to represent steady state solutions, fold and Hopf points, periodic solutions and connecting orbits. Here, n is the system dimension, p is the number of parameters, L is the number of intervals used to represent the periodic solution, d is the degree of the polynomial on each interval, s1 is the number of unstable modes of x − and s2 is the number of unstable modes of x + . values at the points {ti+ j } j=0,...,d . Hence the complete profile is represented by its value at all the mesh points, d x ∗ (ti+ j ), i = 0, . . . , L − 1, j = 0, . . . , d − 1; and x ∗ (t L ). d Because polynomials on adjacent intervals share the value at the common interval point, this representation is automatically continuous (it is, however, not continuously differentiable). (As indicated in table 4, the mesh may be empty, which indicates the use of an equidistant, fixed mesh.) A connecting orbit is represented by the parameter values η, the period T, a time-scaled profile x ∗ (t/T ) on a mesh in [0,1], the steady states x − and x + (fields 'x1' and 'x2' in the data structure), the unstable eigenvalues of these steady states, λ− and λ+ (fields 'lambda_v' and 'lambda_w' in the data structure), the unstable right eigenvectors of x − ('v'), the unstable left eigenvectors of x + ('w'), the direction in which the profile leaves the unstable manifold, determined by α, and the distance of the first point of the profile to x − , determined by e. For the mesh and profile, the same remarks as in the case of periodic solutions hold. The point structures are used as input to the point manipulation routines (layer 2) and are used 17 inside the branch structure (see further). The order of the fields in the point structures is important (because they are used as elements of an array inside the branch structure). No such restriction holds for the other structures (method, plot and branch) described in the rest of this section. 7.3. Stability structures Most of the point structures contain a field 'stability' storing eigenvalues or Floquet multipliers. (The exception is the 'hcli' structure for which stability does not really make sense.) During bifurcation analysis the computation of stability is typically performed as a separate step, after computation of the solution branches, because stability computation can easily be more expensive than the solution finding. If no stability has been computed yet, the field 'stability' is empty, otherwise, it contains the computed stability information in the form described in Table 5. field content field content 'h' 'l0' 'l1' 'n1' R Cn l Cn c ({−1} ∪ N0 )nc or empty 'mu' Cn m (a) Structure in field 'stability' for steady state, fold and Hopf points of Table 4 (b) Structure in field 'stability' for periodic orbit points of Table 4 Table 5: Stability structures for roots of the characteristic equation (in steady state, fold and Hopf structures) (left) and for Floquet multipliers (in the periodic solutions structure) (right). Here, nl is the number of approximated roots, nc is the number of corrected roots and nm is the number of Floquet multipliers. For steady state, fold and Hopf points, approximations to the rightmost roots of the characteristic equation are provided in field 'l0' in order of decreasing real part. The steplength that was used to obtain the approximations is provided in field 'h'. Corrected roots are provided in field 'l1' and the number of Newton iterations applied for each corrected root in a corresponding field 'n1'. If unconverged roots are discarded, 'n1' is empty and the roots in 'l1' are ordered with respect to real part; otherwise the order in 'l1' corresponds to the order in 'l0' and an element −1 in 'n1' signals that no convergence was reached for the corresponding root in 'l0' and the last computed iterate is stored in 'l1'. The collection of uncorrected roots presents more accurate yet less robust information than the collection of approximate roots, see section 10. For periodic solutions only approximations to the Floquet multipliers are provided in a field 'mu' (in order of decreasing modulus). As the characteristic matrix is not analytically available, DDE-BIFTOOL does not offer an additional correction. 7.4. Method parameters To compute a single steady state, fold, Hopf, periodic or connecting orbit solution point, several method parameters have to be passed to the appropriate routines. These parameters are collected into a structure with the fields given in Table 6. For the computation of periodic solutions, additional fields are necessary, marked with and asterisk (∗ ) in Table 6. The meaning of the different fields in Table 6 is explained in section 10. Parameters controlling the pseudo-arclength continuation (using secant approximations for tangents) are stored in a structure of the form given in Table 8. Similarly, for the approximation and correction of roots of the characteristic equation respectively for the computation of the Floquet multipliers method parameters are passed using a structure of the form given in table 7. 18 field 'newton_max_iterations' 'newton_nmon_iterations' 'halting_accuracy' 'minimal_accuracy' 'extra_condition' 'print_residual_info' ∗ 'phase_condition' ∗ 'collocation_parameters' ∗ 'adapt_mesh_before_correct' ∗ 'adapt_mesh_after_correct' content default value N0 N R+ 5, 5, 5, 5, 10 1 1e-10, 1e-9, 1e-9, 1e-8, 1e-8 R0+ {0, 1} {0, 1} 1e-8, 1e-7, 1e-7, 1e-6, 1e-6 0 0 {0, 1} [0, 1]d or empty N N 1 empty 0 3 Table 6: Point method structure: fields and possible values. When different, default values are given in the order 'stst','fold','hopf','psol', 'hcli'. Fields marked with and asterisk (∗ ) are needed and present for points of type 'psol' and 'hcli' only. 7.5. Branch structures A branch consists of an ordered array of points (all of the same type), and three method structures containing point method parameters, continuation parameters respectively stability computation field content default value Rk Rk R0+ N0 R0+ R0+ N0 R or empty N R0+ {0, 1} R0− time_lms('bdf',4) time_lms('bdf',4) time_saf(alpha,beta,0.01,0.01) 4 0.01 0.1 100 empty 6 1e-6 1 -1e-8 For steady state, fold and Hopf 'lms_parameter_alpha' 'lms_parameter_beta' 'lms_parameter_rho' 'interpolation_order' 'minimal_time_step' 'maximal_time_step' 'max_number_of_eigenvalues' 'minimal_real_part' 'max_newton_iterations' 'root_accuracy' 'remove_unconverged_roots' 'delay_accuracy' For periodic orbit 'collocation_parameters' [0, 1]d or empty 'max_number_of_eigenvalues' N 'minimal_modulus' R+ 'delay_accuracy' R0− empty 100 0.01 -1e-8 Table 7: Stability method structures: fields and possible values for the approximation and correction of roots of the characteristic equation (top), or for the approximation Floquet multipliers (bottom). The LMSparameters are default set to the fourth order backwards differentiation LMS-method. The last row in both parts is only used for sd-DDEs. 19 field 'steplength_condition' 'plot' 'prediction' 'steplength_growth_factor' 'plot_progress' 'plot_measure' 'halt_before_reject' content default value {0, 1} {0, 1} {1} R0+ {0, 1} struct or empty {0, 1} 1 1 1 1.2 1 empty 0 Table 8: Continuation method structure: fields and possible values. field subfield content 'point' 'method' 'point' 'stability' 'continuation' 'parameter' 'free' 'min_bound' 'max_bound' 'max_step' array of points (s. Table 4) point method struct stability method struct continuation method struct (s. Table 6) (s. Table 7) (s. Table 8) Np f [N R] pi [N R] p a [N R] p s Table 9: Branch structure: fields and possible values. Here, p f is the number of free parameters; pi , p a and ps are the number of minimal parameter values, maximal parameter values respectively maximal parameter steplength values. If any of these values are zero, the corresponding subfield is empty. parameters, see table 9. The branch structure has three fields. One, called 'point', which contains an array of point structures, one, called 'method', which is itself a structure containing three subfields and a third, called 'parameter' which contains four subfields. The three subfields of the method field are again structures. The first, called 'point', contains point method parameters as described in Table 6. The second, called 'stability', contains stability method parameters as described in Table 7 and the third, called 'continuation', contains continuation method parameters as described in Table 8. Hence the branch structure incorporates all necessary method parameters which are thus automatically kept when saving a branch variable to file. The parameter field contains a list of free parameter numbers which are allowed to vary during computations, and a list of parameter bounds and maximal steplengths. Each row of the bound and steplength subfields consists of a parameter number (first element) and the value for the bound or steplength limitation. Examples are given in demo neuron (see ../demos/index.html). A default, empty branch structure can be obtained by passing a list of free parameters and the point kind (as 'stst', 'fold', 'hopf', 'psol' or 'hcli') to the function df_brnch. A minimal bound zero is then set for each constant delay if the function sys_tau is defined as in section 6.2 (i.e. for DDEs). The method contains default parameters (containing appropriate point, stability and continuation fields) obtained from the function df_mthod with as only argument the type of solution point. 20 field content meaning 'field' {'parameter','x','v','omega', . . . 'profile','period','stability' . . . } 'subfield' {’ ’,'l0','l1','mu'} 'row' N or {'min','max','mean','ampl','all'} 'col' N or {'min','max','mean','ampl','all'} 'func' {'','real','imag','abs'} first field to select from a point struct empty string or 2nd field to select row index column index function to apply Table 10: Measure structure: fields, content and meaning of a structure describing a measure of a point. 7.6. Scalar measure structure After a branch has been computed some possibilities are offered to plot its content. For this a (scalar) measure structure is used which defines what information should be taken and how it should be processed to obtain a measure of a given point (such as the amplitude of the profile of a periodic solution, etc. . . ); see Table 10. The result applied to a variable point is to be interpreted as scalar_measure = func ( point . field . subfield ( row , col )); where 'field' presents the field to select, 'subfield' is empty or presents the subfield to select, ’row’ presents the row number or contains one of the functions mentioned in table 10. These functions are applied columnwise over all rows. The function 'all' specifies that the all rows should be returned. The meaning of 'col' is similar to 'row' but for columns. To avoid ambiguity it is required that either 'row' or 'col' contains a number or that both contain the function 'all'. If nonempty, the function ’func’ is applied to the result. Note that 'func' can be a standard Matlab function as well as a user written function. Note also that, when using the value 'all' in the fields 'col' and/or 'row' it is possible to return a non-scalar measure (possibly but not necessarily further processed by 'func'). 8. Point manipulation Several of the point manipulation routines have already been used in the previous section. Here we outline their functionality and input and output parameters. A brief description of parameters is also contained within the source code and can be obtained in Matlab using the help command. Note that a vector of zero elements corresponds to an empty matrix (written in Matlab as []). function [ point , success ]= p_correc (... funcs , point0 , free_par , step_cnd , method , adapt , previous , d_nr , tz ) Function p_correc corrects a given point. • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • point0: initial, approximate solution point as a point structure (see table 4). • free_par: a vector of zero, one or more free parameters. • step_cnd: a vector of zero, one or more linear steplength conditions. Each steplength condition is assumed fulfilled for the initial point and hence only the coefficients of the condition with respect to all unknowns are needed. These coefficients are passed as a point structure (see table 4). This means that for, e.g., a steady state solution point p the i-th steplength condition enforces that 21 step_cnd ( i ). parameter *( p . parameter - point0 . parameter ) ' +... step_cnd ( i ). x ' *( p .x - point0 . x ) is zero. Similar formulas hold for the other solution types. • method: a point method structure containing the method parameters (see table 6). • adapt (optional): if zero or absent, do not use adaptive mesh selection (for periodic solutions); if one, correct, use adaptive mesh selection and recorrect. • previous (optional): for periodic solutions and connecting orbits: if present and not empty, minimize phase shift with respect to this point. Note that this argument should always be present when correcting solutions for sd-DDEs, since in that case the argument d_nr always needs to be specified. In the case of steady state, fold or Hopf-like points, one can just enter an empty vector. • d_nr: (only for equations with state-dependent delays) if present, number of a negative statedependent delay. • tz: (only for equations with state-dependent delays and periodic solutions) if present, a periodic solution is computed such that τtz = 0 and dτtz /dt = 0, where τ is a negative state-dependent delay with number d_nr. For steady state solutions, a solution corresponding to τ = 0 is computed. • point: the result of correcting point0 using the method parameters, steplength condition(s) and free parameter(s) given. Stability information present in point0 is not passed onto point. If divergence occurred, point contains the final iterate. • success: nonzero if convergence was detected (that is, if the requested accuracy has been reached). function stability = p_stabil ( funcs , point , method ) Function p_stabil computes stability of a given point by approximating its stability-determining eigenvalues. • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • point: a solution point as a point structure (see table 4). • method: a stability method structure (see table 7). • stability: the computed stability of the point through a collection of approximated eigenvalues (as a structure described in table 5). For steady state, fold and Hopf points both approximations and corrections to the rightmost roots of the characteristic equation are provided. For periodic solutions approximations to the dominant Floquet multipliers are computed. function p_splot ( point ) Function p_splot plots the characteristic roots respectively Floquet multipliers of a given point (which should contain nonempty stability information). Characteristic root approximations and Floquet multipliers are plotted using ’×’, corrected characteristic roots using ’∗’. 22 function function function function function function function stst_point = p_tostst ( funcs , point ) fold_point = p_tofold ( funcs , point ) hopf_point = p_tohopf ( funcs , point , excludefreqs ) [ psol_point , stepcond ]= p_topsol ( funcs , point , ampl , degree , nr_int ) [ psol_point , stepcond ]= p_topsol ( funcs , point , ampl , coll_points ) [ psol_point , stepcond ]= p_topsol ( funcs , hcli_point ) hcli_point = p_tohcli ( funcs , point ) The functions p_tostst, p_tofold, p_tohopf, p_topsol and p_tohcli convert a given point into an approximation of a new point of the kind indicated by their name. They are used to switch from a steady state point to a Hopf point or fold point, from a Hopf point to a fold point or vice versa, from a (nearby double) Hopf point to the second Hopf point, from a Hopf point to the emanating branch of periodic solutions, from a periodic solution near a period doubling bifurcation to the period-doubled branch and from a periodic solution near a homoclinic orbit to this homoclinic orbit. The function p_tostst is also capable of extracting the initial and final steady states from a connecting orbit. The additional argument excludefreqs of function p_tohopf controls which Hopf frequency is chosen if several are possible. Function p_tohopf calls p_correc with an initial guess for the Hopf freqency equal to the eigenvalue closest to the imaginary axis. For each entry in the vector excludefreqs (of non-negative real numbers) the eigenvalue with imaginary part closest to this entry will be removed from consideration. This becomes useful in systems with large delays. Then equilibria have a large number of eigenvalues close to the imaginary axis, and correspondingly, the system experiences many Hopf bifurcations over short parameter intervals. These Hopf bifurcations can be picked up by calling p_tohopf repeatedly, and excluding frequencies one after another. When starting a periodic solution branch from a Hopf point, an equidistant mesh is produced with nr_int intervals and piecewise polynomials of degree col_degree and a steplength condition stepcond is returned which should be used (together with a corresponding free parameter) in correcting the returned point. This steplength condition (normally) prevents convergence back to the steady state solution (as a degenerate periodic solution of amplitude zero). When jumping to a period-doubled branch, a period-doubled solution profile is produced using coll_points for collocation points and a mesh which is the (scaled) concatenation of two times the original mesh. A steplength condition is returned which (normally) prevents convergence back to the single period branch. When jumping from a homoclinic orbit to a periodic solution, the steplength condition prevents divergence, by keeping the period fixed. When extracting the steady states from a connecting orbit, an array is returned in which the first element is the initial steady state, and the second element is the final steady state. function rm_point = p_remesh ( point , new_degree , new_mesh ) Function p_remesh changes the piecewise polynomial representation of a given periodic solution point. • point: initial point, containing old mesh, old degree and old profile. • new_degree: new degree of piecewise polynomials. • new_mesh: mesh for new representation of periodic solution profile either as a (non-scalar) row vector of mesh points (both interval and representation points, with the latter chosen equidistant between the former, see section 7) or as the new number of intervals. In the latter case the new mesh is adaptively chosen based on the old profile. • rm_point: returned point containing new degree, new mesh and an appropriately interpolated (but uncorrected!) profile. function tau_eva = p_tau ( funcs , point , d_nr , t ) 23 Function p_tau evaluates state-dependent delay(s) with number(s) d_nr. • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • point: a solution point as a point structure. • d_nr: number(s) of delay(s) (in increasing order) to evaluate. • t (absent for steady state solutions and optional for periodic solutions): mesh (a time point or a number of time points). If present, delay function(s) are evaluated at the points of t, otherwise at the point.mesh (if point.mesh is empty, an equidistant mesh is used). • tau_eva: evaluated values of delays (at t). The following routines are used within branch routines but are less interesting for the general user. function sc_measure = p_measur (p , measure ) Function p_measur computes the (scalar) measure measure of the given point p (see table 10). function p = p_axpy (a ,x , y ) Function p_axpy performs the axpy-operation on points. That is, it computes p=ax+y where a is a scalar, and x and y are two point structures of the same type. p is the result of the operation on all appropriate fields of the given points. If x and y are solutions on different meshes, interpolation is used and the result is obtained on the mesh of x. Stability information, if present, is not passed onto p. function n = p_norm ( point ) Function p_norm computes some norm of a given point structure. function normalized_p = p_normlz ( p ) Function p_normlz performs some normalization on the given point structure p. In particular, fold, Hopf and connecting orbit determining eigenvectors are scaled to norm 1. function [ delay_nr , tz ]= p_tsgn ( point ) Function p_tsgn detects a first negative state-dependent delay. • point: a solution point as a point structure. • delay_nr: number of the first (and only the first !) detected negative delay τ. • tz (only for periodic solutions): tz ∈ [0, 1] is a (time) point such that the delay function τ (t) has its minimal value near this point. To compute tz, a refined mesh is used in the neighbourhood of the minimum of the delay function. This point is later used to compute a periodic solution such that τtz = 0 and dτtz /dt = 0. 9. Branch manipulation Usage of most of the branch manipulation routines is illustrated in the demos neuron and sd demo (see ../demos/index.html). Here we outline their functionality and input and output variables. As for all routines in the package, a brief description of the parameters is also contained within the source code and can be obtained in Matlab using the help command. function [ c_branch , succ , fail , rjct ]= br_contn ( funcs , branch , max_tries ) The function br_contn computes (or rather extends) a branch of solution points. 24 • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • branch: initial branch containing at least two points and computation, stability and continuation method parameter structures and a free parameter structure as described in table 9. • max_tries: maximum number of corrections allowed. • c_branch: the branch returned contains a copy of the initial branch plus the extra points computed (starting from the end of the point array in the initial branch). • succ: number of successful corrections. • fail: number of failed corrections. • rjct: number of rejected points. Note also that successfully computed points are normalized using the procedure p_normlz (see section 8). function br_plot ( branch , x_measure , y_measure , line_type ) Function br_plot plots a branch (in the current figure). • branch: branch to plot (see table 9). • x_measure: (scalar) measure to produce plotting quantities for the x-axis (see table 10). If empty, the point number is used to plot against. • y_measure: (scalar) measure to produce plotting quantities for the y-axis (see table 10). If empty, the point number is used to plot against. • line_type (optional): line type to plot with. function [ x_measure , y_measure ]= df_measr ( stability , branch ) function [ x_measure , y_measure ]= df_measr ( stability , par_list , kind ) Function df_measr returns default measures for plotting. • stability: nonzero if measures are required to plot stability information. • branch: a given branch (see table 9) for which default measures should be constructed. • par_list: a list of parameters for which default measures should be constructed. • kind: a point type for which default measures should be constructed. • x_measure: default scalar measure to use for the x-axis. x_measure is chosen as the first parameter which varies along the branch or as the first parameter of par_list. • y_measure: default scalar measure to use for the y-axis. If stability is zero, the following choices are made for y_measure. For steady state solutions, the first component which varies along the branch; for fold and Hopf bifurcations the first parameter value (different from the one used for x_measure) which varies along the branch. For periodic solutions, the amplitude of the fist varying component. If stability is nonzero, y_measure selects the real part of the characteristic roots (for steady state solutions, fold and Hopf bifurcations) or the modulus of the Floquet multipliers (for periodic solutions). 25 function st_branch = br_stabl ( funcs , branch , skip , recompute ) Function br_stabl computes stability information along a previously computed branch. • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • branch: given branch (see table 9). • skip: number of points to skip between stability computations. That is, computations are performed and stability field is filled in every skip+1-th point. • recompute: if zero, do not recompute stability information present. If nonzero, discard and recompute old stability information present (for points which were not skipped). • st_branch: a copy of the given branch whose (non-skipped) points contain a non-empty stability field with computed stability information (using the method parameters contained in branch). function t_branch = br_rvers ( branch ) To continue a branch in the other direction (from the beginning instead of from the end of its point array), br_rvers reverses the order of the points in the branches point array. function recmp_branch = br_recmp ( funcs , branch , point_numbers ) Function br_recmp recomputes part of a branch. • funcs: structure of user-defined functions, defining the problem (created, for example, using set_funcs). • branch: initial branch (see table 9). • point_numbers (optional): vector of one or more point numbers which should be recomputed. Empty or absent if the complete point array should be recomputed. • recmp_branch: a copy of the initial branch with points who were (successfully) recomputed replaced. If a recomputation fails, a warning message is given and the old value remains present. This routine can, e.g., be used after changing some method parameters within the branch method structures. function [ col , lengths ]= br_measr ( branch , measure ) Function br_selec computes a measure along a branch. • branch: given branch (see table 9). • measure: given measure (see table 10). • col: the collection of measures taken along the branch (over its point array) ordered row-wise. Thus, a column vector is returned if measure is scalar. Otherwise, col contains a matrix. • lengths: vector of lengths of the measures along the branch. If the measure is not scalar, it is possible that its length varies along the branch (e.g. when plotting rightmost characteristic roots). In this situation col is a matrix with number of columns equal to the maximal length of the measures encountered. Extra elements of col are automatically put to zero by Matlab. lengths can then be used to prevent plotting of extra zeros. 26 10. Numerical methods This section contains short descriptions of the numerical methods for DDEs and the method parameters used in DDE-BIFTOOL. More details on the methods can be found in the articles [35, 18, 17, 16, 19, 39] or in [15]. For details on applying these methods to bifurcation analysis of sd-DDEs see [34]. 10.1. Determining systems Below we state the determining systems used to compute and continue steady state solutions, steady state fold and Hopf bifurcations, periodic solutions and connecting orbits of systems of delay differential equations. For each determining system we mention the number of free parameters necessary to obtain (generically) isolated solutions. In the package, the necessary number of free parameters is further raised by the number of steplength conditions plus the number of extra conditions used. This choice ensures the use of square Jacobians during Newton iteration. If, on the other hand, the number of free parameters, steplength conditions and extra conditions are not appropriately matched Newton iteration solves systems with a non-square Jacobian (for which Matlab uses an over- or underdetermined least squares procedure). If possible, it is better to avoid such a situation. A steady state solution x ∗ ∈ Rn is determined from the following ndimensional determining system with no free parameters. Steady state solutions f ( x ∗ , x ∗ , . . . , x ∗ , η ) = 0. (12) Fold bifurcations, ( x ∗ ∈ Rn , v ∈ Rn ) are determined from the following 2n + 1-dimensional determining system using one free parameter. Steady state fold bifurcations 0 = f (x∗ , x∗ , . . . , x∗ , η ) 0 = ∆( x ∗ , η, 0)v (13) T 0 = c v−1 (see (4) for the definition of the characteristic matrix ∆). Here, cT v − 1 = 0 presents a suitable T normalization of v. The vector c ∈ Rn is chosen as c = v(0) /(v(0) v(0) ), where v(0) is the initial value of v. Hopf bifurcations, ( x ∗ ∈ Rn , v ∈ Cn , ω ∈ R) are determined from the following 2n + 1-dimensional partially complex (and by this fact more properly called a 3n + 2dimensional) determining system using one free parameter. Steady state Hopf bifurcations 0 = f (x∗ , x∗ , . . . , x∗ , η ) 0 = ∆( x ∗ , η, iω )v H 0 = c v−1 27 (14) Periodic solutions Periodic solutions are found as solutions (u(s), s ∈ [0, 1]; T ∈ R) of the following (n( Ld + 1) + 1-dimensional system with no free parameters. h h τm τ1 , . . . , u ci,j − , u̇(ci,j ) = T f u(ci,j ), u ci,j − T mod [0,1] T mod [0,1] i = 0, . . . , L − 1, j = 1, . . . , d (15) u (0) = u (1), p(u) = 0. Here the notation t| condition mod [0,1] refers to t − max{k ∈ Z : k ≤ t}, and p represents the integral phase Z 1 0 u̇(s)∆u(s)ds = 0, (16) where u is the current solution and ∆u its correction. The collocation points are obtained as ci,j = ti + c j (ti+1 − ti ), i = 0, . . . , L − 1, j = 1, . . . , d, from the interval points ti , i = 0, . . . , L − 1 and the collocation parameters c j , j = 1, . . . , d. The profile u is discretized as a piecewise polynomial as explained in section 7. This representation has a discontinuous derivative at the interval points. If ci,j coincides with ti the right derivative is taken in (15), if it coincides with ti+1 the left derivative is taken. In other words the derivative taken at ci,j is that of u restricted to [ti , ti+1 ]. Connecting orbits Connecting orbits can be found as solutions of the following determining system with s+ − s− + 1 free parameters, where s+ and s− denote the number of unstable eigenvalues of x + and x − respectively. u̇(ci,j ) = T f (u(ci,j ), u(ci,j − u(c̃) = x − + e s− τ1 τm ), . . . , u(ci,j − ), η ) = 0, T T − ∑ αk v−k eλk Tc̃ , (i = 0, . . . , L − 1, j = 1, . . . , d) c̃ < 0 k =1 − 0 = f (x− , x , η ) 0 = f (x+ , x+ , η ) 0 =∆( x − (17) − , λ− k , η )vk 0 =ckH v− k − 1, (k = 1, . . . , s− ) + 0 =∆ H ( x + , λ+ k , η ) wk 0 =dkH wk+ − 1, (k = 1, . . . , s+ ) G + θ H H 0 =w2k (u(1) − x + ) + ∑ gi wk+ e−λk (θi +τ ) A1 ( x + , η ) u(1 + i ) − x + , T i =1 (k = 1, . . . , s+ ) s− u (0) = x − + e ∑ α k v − k i =1 s− 1 = ∑ | α k |2 i =1 0 = p(u, η ) Again, all arguments of u are taken modulo [0, 1]. We choose the same phase condition as for periodic + solutions and similar normalization of v− k and w + k as in (14). 28 Point method parameters The point method parameters (see table 6) specify the following options. • newton_max_iterations: maximum number of Newton iterations. • newton_nmon_iterations: during a first phase of newton_nmon_iterations+1 Newton iterations the norm of the residual is allowed to increase. After these iterations, corrections are halted upon residual increase. • halting_accuracy: corrections are halted when the norm of the last computed residual is less than or equal to halting_accuracy is reached. • minimal_accuracy: a corrected point is accepted when the norm of the last computed residual is less than or equal to minimal_accuracy. • extra_condition: this parameter is nonzero when extra conditions are provided in a routine sys cond.m which should border the determining systems during corrections. The routine accepts the current point as input and produces an array of condition residuals and corresponding condition derivatives (as an array of point structures) as illustrated below (§10.2). • print_residual_info: when nonzero, the Newton iteration number and resulting norm of the residual are printed to the screen during corrections. For periodic solutions and connecting orbits, the extra mesh parameters (see table 6) provide the following information. • phase_condition: when nonzero the integral phase condition (16) is used. • collocation_parameters: this parameter contains user given collocation parameters. When empty, Gauss-Legendre collocation points are chosen. • adapt_mesh_before_correct: before correction and if the mesh inside the point is nonempty, adapt the mesh every adapt_mesh_before_correct points. E.g.: if zero, do not adapt; if one, adapt every point; if two adapt the points with odd point number. • adapt_mesh_after_correct: similar to adapt_mesh_before_correct but adapt mesh after successful corrections and correct again. 10.2. Extra conditions When correcting a point or computing a branch, it is possible to add one or more extra conditions and corresponding free parameters to the determining systems presented earlier. These extra conditions should be implemented using a function sys_cond and setting the method parameter extra_condition to 1 (cf. table 6). The function sys_cond accepts the current point as input and produces a residual and corresponding condition derivatives (as a point structure) per extra condition. As an example, suppose we want to compute a branch of periodic solutions of system (10) subject to the following extra conditions T = 200, 0 = a212 + a221 − 1, (18) that is, we wish to continue a branch with fixed period T = 200 and parameter dependence a212 + a221 = 1. The routine shown in Listing 5 implements these conditions by evaluating and returning each residual for the given point and the derivatives of the conditions w.r.t. all unknowns (that is, w.r.t. to all the components of the point structure). 29 function [ resi , condi ]= sys_cond ( point ) % kappa beta a12 a21 tau1 tau2 tau_s if point . kind == ' psol ' % fix period at 200: resi (1)= point . period -200; % derivative of first condition wrt unknowns : condi (1)= p_axpy (0 , point ,[]); condi (1). period =1; % parameter condition : resi (2)= point . parameter (3)^2+ point . parameter (4)^2 -1; % derivative of second condition wrt unknowns : condi (2)= p_axpy (0 , point ,[]); condi (2). parameter (3)=2* point . parameter (3); condi (2). parameter (4)=2* point . parameter (4); else error ( ' SYS_COND : point is not psol . ' ); end end Listing 5: Implementation extra conditions (18) using a routine sys_cond. 10.3. Continuation During continuation, a branch is extended by a combination of predictions and corrections. A new point is predicted based on previously computed points using secant prediction over an appropriate steplength. The prediction is then corrected using the determining systems (12), (13), (14), (15) or (17) bordered with a steplength condition which requires orthogonality of the correction to the secant vector. Hence one extra free parameter is necessary compared to the numbers mentioned in the previous section. The following continuation and steplength determination strategy is used. If the last point was successfully computed, the steplength is multiplied with a given, constant factor greater than 1. If corrections diverged or if the corrected point was rejected because its accuracy was not acceptable, a new point is predicted, using linear interpolation, halfway between the last two successfully computed branch points. If the correction of this point succeeds, it is inserted in the point array of the branch (before the previously last computed point). If the correction of the interpolated point fails again, the last successfully computed branch point is rejected (for fear of branch switch) and the interpolation procedure is repeated between the (new) last two branch points. Hence, if, after a failure, the interpolation procedure succeeds, the steplength is approximately divided by a factor two. Test results indicate that this procedure is quite effective and proves an efficient alternative to using only (secant) extrapolation with steplength control. The reason for this is mainly that the secant extrapolation direction is not influenced by halving the steplength but it is by inserting a newly computed point in between the last two computed points. Continuation method parameters The continuation method parameters (see table 8) have the fol- lowing meaning. • plot: if nonzero, plot predictions and corrections during continuation. • prediction: this parameter should be 1, indicating that secant prediction is used (being currently the only alternative). 30 • steplength_growth_factor: grow the steplength with this factor in every step except during interpolation. • plot_progress: if nonzero, plotting is visible during continuation process. If zero, only the final result is drawn. • plot_measure: if empty use default measures to plot. Otherwise plot_measure contains two fields, ’x’ and ’y’, which contain measures (see table 10) for use in plotting during continuation. • halt_before_reject: If this parameter is nonzero, continuation is halted whenever (and instead of) rejecting a previously accepted point based on the above strategy. 10.4. Roots of the characteristic equation Roots of the characteristic equation are approximated using a linear multi-step (LMS-) method applied to (2). Consider the linear k-step formula k k j =0 j =0 ∑ α j y L+ j = h ∑ β j f L+ j . (19) Here, α0 = 1, h is a (fixed) step size and y j presents the numerical approximation of y(t) at the mesh point t j := jh. The right hand side f j := f (y j , ỹ(t j − τ1 ), . . . , ỹ(t j − τm )) is computed using approximations ỹ(t j − τ1 ) obtained from yi in the past, i < j. In particular, the use of so-called Nordsieck interpolation, leads to s ỹ(t j + eh) = ∑ Pl (e)y j+l , e ∈ [0, 1). (20) l =−r using s Pl (e) := ∏ k =−r, k 6=l e−k . l−k The resulting method is explicit whenever β 0 = 0 and min τi > sh. That is, y L+k can then directly be computed from (19) by evaluating y L+k = − k −1 k j =0 j =0 ∑ α j y L+ j + h ∑ β j f L+ j . whose right hand side depends only on y j , j < L + k. For the linear variational equation (2) around a steady state solution x ∗ (t) ≡ x ∗ we have m f j = A0 y j + ∑ Ai ỹ(t j − τi ) (21) i =0 where we have omitted the dependency of Ai on x ∗ . The stability of the difference scheme (19), (21) can be evaluated by setting y j = µ j− Lmin , j = Lmin , . . . , L + k where Lmin is the smallest index used, taking the determinant of (19) and computing the roots µ. If the roots of the polynomial in µ all have modulus smaller than unity, the trajectories of the LMS-method converge to zero. If roots exist with modulus greater than unity then trajectories exist which grow unbounded. Since the LMS-method forms an approximation of the time integration operator over the time step h, so do the roots µ approximate the eigenvalues of S(h, 0). The eigenvalues of S(h, 0) are exponential transforms of the roots λ of the characteristic equation (5), µ = exp(λh). 31 Hence, once µ is found, λ can be extracted using, Re(λ) = ln(|µ|) . h (22) The imaginary part of λ is found modulo π/h, using Im(λ) ≡ arcsin( Im(µ) ) |µ| h (mod π ). h (23) For small h, 0 < h 1, the smallest representation in (23) is assumed the most accurate one (that is, we let arcsin map into [−π/2, π/2]). The parameters r and s (from formula (20)) are chosen such that r ≤ s ≤ r + 2 (see [28]). The choice of h is based on the related heuristic outlined in [19]. Approximations for the rightmost roots λ obtained from the LMS-method using (22), (23) can be corrected using a Newton process on the system, 0 = ∆(λ)v (24) 0 = cT v − 1 A starting value for v is the eigenvector of ∆(λ) corresponding to its smallest eigenvalue (in modulus). Note that the collection of successfully corrected roots presents more accurate yet less robust information than the set of uncorrected roots. Indeed, attraction domains of roots of equations like (24) can be very small and hence corrections may diverge, or different roots may be corrected to a single ’exact’ root thereby missing part of the spectrum. The latter does not occur when computing the (full) spectrum of a discretization of S(h, 0). Stability information is kept in the structure of table 5 (left). The time step used is kept in field h. Approximate roots are kept in field l0, corrected roots in field l1. If unconverged corrected roots are discarded, field n1 is empty. Otherwise, the number of Newton iterations used is kept for each root in the corresponding position of n1. Here, −1 signals that convergence to the required accuracy was not reached. Stability method parameters The stability method parameters (see table 7 (top)) now have the following meaning. • lms_parameter_alpha: LMS-method parameters α j ordered from past to present, j = 0, 1, . . . , k. • lms_parameter_beta: LMS-method parameters β j ordered from past to present, j = 0, 1, . . . , k. • lms_parameter_rho: safety radius ρLMS,e of the LMS-method stability region. For a precise definition, see [15, §III.3.2]. • interpolation_order: order of the interpolation in the past, r + s = interpolation_order. • minimal_time_step: minimal time step relative to maximal delay, h τ ≥ minimal_time_step. • maximal_time_step: maximal time step relative to maximal delay, h τ ≤ minimal_time_step. • max_number_of_eigenvalues: maximum number of rightmost eigenvalues to keep. • minimal_real_part: choose h such as the discretized system approximates eigenvalues with Re(λ) ≥ minimal_real_part well, discard eigenvalues with Re(λ) < minimal_real_part. If h is smaller than its minimal value, it is set to the minimal value and a warning is given. If it is larger than its maximal value it is reduced to that number without warning. If minimal and maximal value coincide, h is set to this value without warning. If minimal_real_part is empty, the value minimal_real_part = τ1 is used. 32 • max_newton_iterations: maximum number of Newton iterations during the correction process (24). • root_accuracy: required accuracy of the norm of the residual of (24) during corrections. • remove_unconverged_roots: if this parameter is zero, unconverged roots are discarded (and stability field n1 is empty). • delay_accuracy (only for state-dependent delays): if the value of a state-dependent delay is less than delay_accuracy, the stability is not computed. 10.5. Floquet multipliers Floquet multipliers are computed as eigenvalues of the discretized time integration operator S( T, 0). The discretization is obtained using the collocation equations (15) without the modulo operation (and without phase and periodicity condition). From this system a discrete, linear map is obtained between the variables presenting the segment [−τ/T, 0] and those presenting the segment [−τ/T + 1, 1]. If these variables overlap, part of the map is just a time shift. Stability information is kept in the structure of table 5 (right). Approximations to the Floquet multipliers are kept in field mu. Stability method parameters The stability method parameters (see table 7 (bottom)) have the fol- lowing meaning. • collocation_parameters: user given collocation parameters or empty for Gauss-Legendre collocation points. • max_number_of_eigenvalues: maximum number of multipliers to keep. • minimal_modulus: discard multipliers with |µ| < minimal_modulus. • delay_accuracy (only for state-dependent delays): if the value of a state-dependent delay is less than delay_accuracy, the stability is not computed. 11. Concluding comments The first aim of DDE-BIFTOOL is to provide a portable, user-friendly tool for numerical bifurcation analysis of steady state solutions and periodic solutions of systems of delay differential equations of the kinds (1) and (7). Part of this goal was fulfilled through choosing the portable, programmerfriendly environment offered by Matlab. Robustness with respect to the numerical approximation is achieved through automatic steplength selection in approximating the rightmost characteristic roots and through collocation using piecewise polynomials combined with adaptive mesh selection. Although the package has been successfully tested on a number of realistic examples, a word of caution may be appropriate. First of all, the package is essentially a research code (hence we accept no reliability) in a quite unexplored area of current research. In our experience up to now, new examples did not fail to produce interesting theoretical questions (e.g., concerning homoclinic or heteroclinic solutions) many of which remain unsolved today. Unlike for ordinary differential equations, discretization of the state space is unavoidable during computations on delay equations. Hence the user of the package is strongly advised to investigate the effect of discretization using tests on different meshes and with different method parameters; and, if possible, to compare with analytical results and/or results obtained using simulation. Although there are no ’hard’ limits programmed in the package (with respect to system and/or mesh sizes), the user will notice the rapidly increasing computation time for increasing system 33 dimension and mesh sizes. This is most notable in the stability and periodic solution computations. Indeed, eigenvalues are computed from large sparse matrices without exploiting sparseness and the Newton procedure for periodic solutions is implemented using direct methods. Nevertheless the current version is sufficient to perform bifurcation analysis of systems with reasonable properties in reasonable execution times. Furthermore, we hope future versions will include routines which scale better with the size of the problem. 11.1. Existing extensions • Extension debiftool extra psol continues the three local codimension-one bifurcations of periodic orbits, the fold bifurcation, the period doubling and the torus bifurcation for constant and state-dependent delays. • Extension debiftool extra rotsym continues relative equilibria and relative periodic orbits and their local codimension-one bifurcations for constant delays in systems with rotational symmetry (that is, there exists a matrix A ∈ Rn×n such that A T = − A and exp( At) f ( x0 , . . . , xm ) = f (exp( At) x0 , . . . , exp( At) xm )). • Extension debiftool extra nmfm computes normal form coefficients of Hopf bifurcations, Hopf-Hopf interactions, generalized Hopf (Bautin) bifurcations, and zero-Hopf interactions (Gavrilov-Guckenheimer bifurcations) for equations with constant delays. Other possible future developments include • a graphical user interface; • incorporation of the numerical core routines into a general continuation framework such as Coco [8] (which would permit the user to grow higher-dimensional solution families and wrap other continuation algorithms around the core DDE routines), • the extension to other types of delay equations such as distributed delay and neutral functional differential equations. See also Barton et al [3] for a demonstration of how to extend DDEBIFTOOL to neutral functional differential equations. • determination of more normal-form coefficients to detect other co-dimension-two bifurcations. Acknowledgements DDE-BIFTOOL v. 2.03 is a result of the research project OT/98/16, funded by the Research Council K.U.Leuven; of the research project G.0270.00 funded by the Fund for Scientific Research - Flanders (Belgium) and of the research project IUAP P4/02 funded by the programme on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister’s Office for Science, Technology and Culture. K. Engelborghs is a Postdoctoral Fellow of the Fund for Scientific Research - Flanders (Belgium). J. Sieber’s contribution to the revision leading to version 3.0 was supported by EPSRC grant EP/J010820/1. References [1] J. Argyris, G. Faust, and M. Haase. An Exploration of Chaos — An Introduction for Natural Scientists and Engineers. North Holland Amsterdam, 1994. [2] N. V. Azbelev, V. P. Maksimov, and L. F. Rakhmatullina. Introduction to the Theory of Functional Differential Equations. Nauka, Moscow, 1991. (in Russian). 34 [3] D.A.W. Barton, B. Krauskopf, and R.E. Wilson. Collocation schemes for periodic solutions of neutral delay differential equations. Journal of Difference Equations and Applications, 12(11):1087– 1101, 2006. [4] R. Bellman and K. L. Cooke. Differential-Difference Equations, volume 6 of Mathematics in science and engineering. Academic Press, 1963. [5] D. Breda, S. Maset, and R. Vermiglio. TRACE-DDE: a tool for robust analysis and characteristic equations for delay differential equations. In J.J. Loiseau, W. Michiels, S.-I. Niculescu, and R. Sipahi, editors, Topics in Time Delay Systems: Analysis, Algorithms, and Control, volume 388 of Lecture Notes in Control and Information Sciences, pages 145–155, New York, 2009. Springer. [6] S.-N. Chow and J. K. Hale. Methods of Bifurcation Theory. Springer-Verlag, 1982. [7] S. P. Corwin, D. Sarafyan, and S. Thompson. DKLAG6: A code based on continuously imbedded sixth order Runge-Kutta methods for the solution of state dependent functional differential equations. Applied Numerical Mathematics, 24(2–3):319–330, 1997. [8] H. Dankowicz and F. Schilder. Recipes for Continuation. Computer Science and Engineering. SIAM, 2013. [9] A Dhooge, W Govaerts, and Y A Kuznetsov. MatCont: A Matlab package for numerical bifurcation analysis of ODEs. ACM Transactions on Mathematical Software, 29(2):141–164, 2003. [10] O. Diekmann, S. A. van Gils, S. M. Verduyn Lunel, and H.-O. Walther. Delay Equations: Functional-, Complex-, and Nonlinear Analysis, volume 110 of Applied Mathematical Sciences. Springer-Verlag, 1995. [11] 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 for Dynamical Systems: Path following and boundary value problems, pages 1–49. Springer-Verlag, Dordrecht, 2007. [12] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. AUTO97: Continuation and bifurcation software for ordinary differential equations; available by FTP from ftp.cs.concordia.ca in directory pub/doedel/auto. [13] R. D. Driver. Ordinary and Delay Differential Equations, volume 20 of Applied Mathematical Science. Springer-Verlag, 1977. [14] L. E. El’sgol’ts and S. B. Norkin. Introduction to the Theory and Application of Differential Equations with Deviating Arguments, volume 105 of Mathematics in science and engineering. Academic Press, 1973. [15] K. Engelborghs. Numerical Bifurcation Analysis of Delay Differential Equations. PhD thesis, Department of Computer Science, Katholieke Universiteit Leuven, Leuven, Belgium, May 2000. [16] K. Engelborghs and E. Doedel. Stability of piecewise polynomial collocation for computing periodic solutions of delay differential equations. Numerische Mathematik, 2001. Accepted. [17] K. Engelborghs, T. Luzyanina, K. J. in ´t Hout, and D. Roose. Collocation methods for the computation of periodic solutions of delay differential equations. SIAM J. Sci. Comput., 22:1593– 1609, 2000. [18] K. Engelborghs and D. Roose. Numerical computation of stability and detection of Hopf bifurcations of steady state solutions of delay differential equations. Advances in Computational Mathematics, 10(3–4):271–289, 1999. 35 [19] K. Engelborghs and D. Roose. On stability of LMS-methods and characteristic roots of delay differential equations. SIAM J. Num. Analysis, 2001. Accepted. [20] W. H. Enright and H. Hayashi. A delay differential equation solver based on a continuous Runge-Kutta method with defect control. Numer. Algorithms, 16:349–364, 1997. [21] B. Ermentrout. XPPAUT3.91 - The differential equations tool. University of Pittsburgh, Pittsburgh, (http://www.pitt.edu/∼phase/) 1998. [22] W.J.F. Govaerts. Numerical Methods for Bifurcations of Dynamical Equilibria. Miscellaneous Titles in Applied Mathematics Series. SIAM, 2000. [23] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer-Verlag New York, 1983. [24] N. Guglielmi and E. Hairer. Stiff delay equations. Scholarpedia, 2(11):2850, 2007. [25] J. K. Hale. Theory of Functional Differential Equations, volume 3 of Applied Mathematical Sciences. Springer-Verlag, 1977. [26] J. K. Hale and S. M. Verduyn Lunel. Introduction to Functional Differential Equations, volume 99 of Applied Mathematical Sciences. Springer-Verlag, 1993. [27] F. Hartung, T. Krisztin, H.-O. Walther, and J. Wu. Functional differential equations with statedependent delays: theory and applications. In P. Drábek, A. Cañada, and A. Fonda, editors, Handbook of Differential Equations: Ordinary Differential Equations, volume 3, chapter 5, pages 435–545. North-Holland, 2006. [28] T. Hong-Jiong and K. Jiao-Xun. The numerical stability of linear multistep methods for delay differential equations with many delays. SIAM Journal of Numerical Analysis, 33(3):883–889, June 1996. [29] The MathWorks Inc. MATLAB. Natick, Massachusetts, United States. [30] V. Kolmanovskii and A. Myshkis. Applied Theory of Functional Differential Equations, volume 85 of Mathematics and Its Applications. Kluwer Academic Publishers, 1992. [31] V. B. Kolmanovskii and A. Myshkis. Introduction to the theory and application of functional differential equations, volume 463 of Mathematics and its applications. Kluwer Academic Publishers, 1999. [32] V. B. Kolmanovskii and V. R. Nosov. Stability of functional differential equations, volume 180 of Mathematics in Science and Engineering. Academic Press, 1986. [33] Y. A Kuznetsov. Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences. Springer-Verlag, New York, third edition, 2004. [34] T. Luzyanina, K. Engelborghs, and D. Roose. Numerical bifurcation analysis of differential equations with state-dependent delay. Internat. J. Bifur. Chaos, 11(3):737–754, 2001. [35] T. Luzyanina and D. Roose. Numerical stability analysis and computation of Hopf bifurcation points for delay differential equations. Journal of Computational and Applied Mathematics, 72:379– 392, 1996. [36] J. Mallet-Paret and R. D. Nussbaum. Stability of periodic solutions of state-dependent delaydifferential equations. Journal of Differential Equations, 250(11):4085 – 4103, 2011. 36 [37] C. A. H. Paul. A user-guide to Archi - an explicit Runge-Kutta code for solving delay and neutral differential equations. Technical Report 283, The University of Manchester, Manchester Center for Computational Mathematics, December 1997. [38] D. Roose and R. Szalai. Continuation and bifurcation analysis of delay differential equations. In B Krauskopf, H M Osinga, and J Galán-Vioque, editors, Numerical Continuation Methods for Dynamical Systems: Path following and boundary value problems, pages 51–75. Springer-Verlag, Dordrecht, 2007. [39] G. Samaey, K. Engelborghs, and D. Roose. Numerical computation of connecting orbits in delay differential equations. Report TW 329, Department of Computer Science, K.U.Leuven, Leuven, Belgium, October 2001. [40] R. Seydel. Practical Bifurcation and Stability Analysis — From Equilibrium to Chaos, volume 5 of Interdisciplinary Applied Mathematics. Springer-Verlag Berlin, 2 edition, 1994. [41] L. F. Shampine and S. Thompson. Solving delay differential equations with dde23. Submitted, 2000. [42] L. P. Shayer and S. A. Campbell. Stability, bifurcation and multistability in a system of two coupled neurons with multiple time delays. SIAM J. Applied Mathematics, 61(2):673–700, 2000. [43] J. Sieber. Finding periodic orbits in state-dependent delay differential equations as roots of algebraic equations. Discrete and Continuous Dynamical Systems A, 32(8):2607–2651, 2012. [44] R. Szalai, G. Stépán, and S.J. Hogan. Continuation of bifurcations in periodic delay differential equations using characteristic matrices. SIAM Journal on Scientific Computing, 28(4):1301–1317, 2006. [45] K. Verheyden, T. Luzyanina, and D. Roose. Efficient computation of characteristic roots of delay differential equations using lms methods,. Journal of Computational and Applied Mathematics, 214:209–226, 2008. [46] Bram Wage. Normal form computations for delay differential equations in DDE-BIFTOOL. Master’s thesis, Utrecht University, 2014. supervised by Y.A. Kuznetsov. 37 A. Jacobians of tutorial examples neuron and sd_demo The meaning of input arguments to sys_deri is explained in section 6.2.3. The analysis of tutorial example neuron is demonstrated in demo neuron step by step (see ../demos/index.html). function J = neuron_sys_deri ( xx , par , nx , np , v ) % % Parameter vector % par =[\ kappa , \ beta , a_ {12} , a_ {21} ,\ tau_1 ,\ tau_2 , \ tau_s ]. J =[]; if length ( nx )==1 && isempty ( np ) && isempty ( v ) % % First - order derivatives wrt to state nx +1 if nx ==0 % derivative wrt x ( t ) J (1 ,1)= - par (1); J (2 ,2)= - par (1); elseif nx ==1 % derivative wrt x (t - tau1 ) J (2 ,1)= par (4)*(1 - tanh ( xx (1 ,2))^2); J (2 ,2)=0; elseif nx ==2 % derivative wrt x (t - tau2 ) J (1 ,2)= par (3)*(1 - tanh ( xx (2 ,3))^2); J (2 ,2)=0; elseif nx ==3 % derivative wrt x (t - tau_s ) J (1 ,1)= par (2)*(1 - tanh ( xx (1 ,4))^2); J (2 ,2)= par (2)*(1 - tanh ( xx (2 ,4))^2); end elseif isempty ( nx ) && length ( np )==1 && isempty ( v ) % % First order derivatives wrt parameters if np ==1 % derivative wrt kappa J (1 ,1)= - xx (1 ,1); J (2 ,1)= - xx (2 ,1); elseif np ==2 % derivative wrt beta J (1 ,1)= tanh ( xx (1 ,4)); J (2 ,1)= tanh ( xx (2 ,4)); elseif np ==3 % derivative wrt a12 J (1 ,1)= tanh ( xx (2 ,3)); J (2 ,1)=0; elseif np ==4 % derivative wrt a21 J (2 ,1)= tanh ( xx (1 ,2)); elseif np ==5 || np ==6 || np ==7 % derivative wrt tau J = zeros (2 ,1); end ; elseif length ( nx )==1 && length ( np )==1 && isempty ( v ) % % Mixed state , parameter derivatives if nx ==0 % derivative wrt x ( t ) if np ==1 % derivative wrt beta J (1 ,1)= -1; J (2 ,2)= -1; else J = zeros (2); end ; elseif nx ==1 % derivative wrt x (t - tau1 ) if np ==4 % derivative wrt a21 J (2 ,1)=1 - tanh ( xx (1 ,2))^2; 38 J (2 ,2)=0; else J = zeros (2); end ; elseif nx ==2 % derivative wrt x (t - tau2 ) if np ==3 % derivative wrt a12 J (1 ,2)=1 - tanh ( xx (2 ,3))^2; J (2 ,2)=0; else J = zeros (2); end ; elseif nx ==3 % derivative wrt x (t - tau_s ) if np ==2 % derivative wrt beta J (1 ,1)=1 - tanh ( xx (1 ,4))^2; J (2 ,2)=1 - tanh ( xx (2 ,4))^2; else J = zeros (2); end ; end ; elseif length ( nx )==2 && isempty ( np ) && ~ isempty ( v ) % % Second order derivatives wrt state variables if nx (1)==0 % first derivative wrt x ( t ) J = zeros (2); elseif nx (1)==1 % first derivative wrt x (t - tau1 ) if nx (2)==1 th = tanh ( xx (1 ,2)); J (2 ,1)= -2* par (4)* th *(1 - th * th )* v (1); J (2 ,2)=0; else J = zeros (2); end ; elseif nx (1)==2 % derivative wrt x (t - tau2 ) if nx (2)==2 th = tanh ( xx (2 ,3)); J (1 ,2)= -2* par (3)* th *(1 - th * th )* v (2); J (2 ,2)=0; else J = zeros (2); end elseif nx (1)==3 % derivative wrt x (t - tau_s ) if nx (2)==3 th1 = tanh ( xx (1 ,4)); J (1 ,1)= -2* par (2)* th1 *(1 - th1 * th1 )* v (1); th2 = tanh ( xx (2 ,4)); J (2 ,2)= -2* par (2)* th2 *(1 - th2 * th2 )* v (2); else J = zeros (2); end end end % Raise error if the requested derivative does not exist 39 if isempty ( J ) error ([ ' SYS_DERI : requested derivative nx =% d , np =% d , ' ,... ' size ( v )=% d could not be computed ! ' ] , nx , np , size ( v )); end end Listing 6: Jacobians of right-hand side neuron_sys_rhs in section 6.2.1 for neuron_sys_rhs. The meaning of input arguments to sys_dtau is explained in section 6.3.4. The analysis of tutorial example sd_demo is demonstrated step by step in demo sd demo (see ../demos/index.html). function dtau = sd_dtau ( ind , xx , par , nx , np ) % p1 p2 p3 p4 p5 p6 p7 p8 p9 p10 p11 dtau =[]; if length ( nx )==1 && isempty ( np ) % % first order derivatives wrt state variables : if nx ==0 % derivative wrt x ( t ) if ind ==3 dtau (1:5)=0; dtau (2)= par (5)* par (10)* xx (2 ,2); elseif ind ==4 dtau (1)= xx (2 ,3)/(1+ xx (1 ,1)* xx (2 ,3))^2; dtau (2:5)=0; elseif ind ==5 dtau (1:5)=0; dtau (4)=1; elseif ind ==6 dtau (5)=1; else dtau (1:5)=0; end ; elseif nx ==1 % derivative wrt x (t - tau1 ) if ind ==3 dtau (1:5)=0; dtau (2)= par (5)* par (10)* xx (2 ,1); else dtau (1:5)=0; end ; elseif nx ==2 % derivative wrt x (t - tau2 ) if ind ==4 dtau (1:5)=0; dtau (2)= xx (1 ,1)/(1+ xx (1 ,1)* xx (2 ,3))^2; else dtau (1:5)=0; end ; else dtau (1:5)=0; end ; elseif isempty ( nx ) && length ( np )==1 , % % First order derivatives wrt parameters if ind ==1 && np ==10 dtau =1; elseif ind ==2 && np ==11 40 dtau =1; elseif ind ==3 && np ==5 dtau = par (10)* xx (2 ,1)* xx (2 ,2); elseif ind ==3 && np ==10 dtau = par (5)* xx (2 ,1)* xx (2 ,2); else dtau =0; end ; elseif length ( nx )==2 && isempty ( np ) , % % Second order derivatives wrt state variables dtau = zeros (5); if ind ==3 if ( nx (1)==0 && nx (2)==1) || ( nx (1)==1 && nx (2)==0) dtau (2 ,2)= par (5)* par (10); end elseif ind ==4 if nx (1)==0 && nx (2)==0 dtau (1 ,1)= -2* xx (2 ,3)* xx (2 ,3)/(1+ xx (1 ,1)* xx (2 ,2))^3; elseif nx (1)==0 && nx (2)==2 dtau (1 ,2)=(1 - xx (1 ,1)* xx (2 ,3))/(1+ xx (1 ,1)* xx (2 ,2))^3; elseif nx (1)==2 && nx (2)==0 dtau (2 ,1)=(1 - xx (1 ,1)* xx (2 ,3))/(1+ xx (1 ,1)* xx (2 ,2))^3; elseif nx (1)==2 && nx (2)==2 dtau (2 ,2)= -2* xx (1 ,1)* xx (1 ,1)/(1+ xx (1 ,1)* xx (2 ,2))^3; end end elseif length ( nx )==1 && length ( np )==1 , % % Mixed state parameter derivatives dtau (1:5)=0; if ind ==3 if nx ==0 && np ==5 dtau (2)= par (10)* xx (2 ,2); elseif nx ==0 && np ==10 dtau (2)= par (5)* xx (2 ,2); elseif nx ==1 && np ==5 dtau (2)= par (10)* xx (2 ,1); elseif nx ==1 && np ==10 dtau (2)= par (5)* xx (2 ,1); end end end % % Otherwise raise error % Raise error if the requested derivative does not exist if isempty ( dtau ) error ([ ' SYS_DTAU , delay % d : requested derivative ' ,... ' nx =% d , np =% d does not exist ! ' ] , ind , nx , np ); end end Listing 7: Jacobians of the delay function sd_tau in Listing 4 for sd-DDE tutorial example (11). 41 B. Octave compatibility considerations nargin incompatibility The core DDE-BifTool code of version v2.0x was likely octave compatible. The changes to 3.1.1, replacing function names by function handles, broke this compatibility initially, because, for example, the call nargin(sys_tau) gives an error message in octave (version 3.2.3) if sys_tau is a function handle. To remedy this problem the additional field tp_del is attached to the structure funcs defining the problem. The field funcs.tp_del is set in set_funcs (see Section 7.1 and Table 3). As of version 3.8.1, octave also gives an error message when one loads function handles as created using set_funcs from a data file. Function handles have to be re-created after a clear or a restart of the session. For an up-to-date list of known differences in syntax and semantics between matlab and octave see http://www.gnu.org/software/octave. Output The gradual updating of plots using drawnow slows down for the gnuplot3 -based plot interface of octave (as of version 3.2.3) as points get added to the plot. Setting the field continuation.plot to 0.5 (that is, less than 1 but larger than 0), prints the values on the screen instead of updating the plot. The fltk-based plotting interface of octave does not appear to experience this slow-down. Useful options to be set in octave: • graphics_toolkit('fltk') sets the graphics toolkit to the fltk-based interface (faster plotting); • page_output_immediately(true); prints out the results of any fprintf or disp commands immediately; • page_screen_output(false); stops paging the terminal output. 3 http://www.gnuplot.info/ 42
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 44 Page Mode : UseOutlines Author : Title : Subject : Creator : LaTeX with hyperref package Producer : pdfTeX-1.40.14 Create Date : 2017:04:05 01:29:42+01:00 Modify Date : 2017:04:05 01:29:42+01:00 Trapped : False PTEX Fullbanner : This is pdfTeX, Version 3.1415926-2.5-1.40.14 (TeX Live 2013/Debian) kpathsea version 6.1.1EXIF Metadata provided by EXIF.tools