\% ********** Class ODEs ********** \topic ode_Class {ordinary differential equations} The class {\bf generic autonomous ODEs} comprises systems of the first order differential equations which can by written in the following form: \par \qc{$\dot x=F(x,\alpha),$} \par where $x\in {\bf R}^n$ is a vector of phase variables, $\alpha\in {\bf R}^m$ is a vector of numerical parameters, and $F$ is a vector-function (right-hand sides). $F$ is assumed to be sufficiently smooth but has no other special properties (like symmetry, etc.). \par A system of ODEs is specified or edited in the \jump kh_RHSWin {System specification dialog box} which is opened by the {\bf Select|System} command from the \jump kh_MainWin {Main window.} \par There are following class-specific fields in the dialog box.\par {\desc 3 {\bf Coordinates}\>is where you list names of phase (state) variables. See \jump kh_RHSWin {System specification dialog box} on how to fill in the fields. The derivative with respect to {\bf time} in the left-hand side of the system is denoted by prime ($\prime$). \par {\bf Parameters}\>is where you list names of parameters your system depends on.\par {\bf Time}\>is where you type in the name of time variable which must be a simple name.\par } The following constants are defined and may be used in RHS, its derivatives and local functions:\par {\desc 2 {\bf COORDINATESDIM}\>the number of coordinates (counting array components).\par {\bf PARAMETERSDIM}\>the number of parameters.\par } \par {\bf Example 1:} For the scalar ODE \par \qc{$\dot{u}=\lambda e^{u},$} \par the following specification is valid. \par {\li 5 \desc 4 {\bf Coordinates:} \>$U$\par {\bf Parameters:} \>$lambda$\par {\bf Time:}\> $t$\par {\bf RHS:} \>$U\prime=lambda*exp(U);$\par } \par {\bf Example 2:} For the system of two ODEs \par \qc{$ \system 2 rcl 1.2 1 {\dot{\xi}_1}{=}{\xi_2} {\dot{\xi}_2}{=}{a\xi^{2}_{1}+ b \xi_{1} \xi_{2}} $} \par the following specification is valid. \par {\li 5 \desc 5 {\bf Coordinates:} \>$x[2]$\par {\bf Parameters:} \>$a,b$\par {\bf Time:}\> $t$\par {\bf RHS:} \>$x\prime [0]=x[1];$\par {} \>$x\prime [1]=a*x[0]*x[0]+b*x[0]*x[1];$\par } \end \% ode_Class \topic ode_Types {types of curves} CONTENT currently supports the continuation of the following curves: \par \par \jump ode_or {Orbit} \par \jump ode_eq {Equilibrium curve} \par \jump ode_lc {Limit cycle curve} \par \jump ode_lp {Limit point curve} \par \jump ode_H {Hopf curve} \par \jump ode_DH {double Hopf curve} \end \% ode_Types \topic ode_eq {equilibrium curve} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot{x}=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^1, $ } \qr {(1)} \par an {\bf equilibrium curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+1}$ endowed with coordinates $(x,\alpha)$ defined by\par \qc{ $ f(x,\alpha)=0. $ }\qr {(2)}\par If one introduces $y=(x,\alpha) \in {\bf R}^{n+1},$ the manifold (2) is given by\par \qc{ $ F(y)=0, $ } \qr {(3)}\par where $F:\: {\bf R}^{n+1} \rightarrow {\bf R}^n,$ $F(y)=f(x,\alpha),$ are the \jump kh_DefFunc {defining functions}. The Jacobian matrix of (3) involved in the continuation has the form \par \qc{$ F_y(y)=(f_x(x,\alpha)|f_{\alpha}(x,\alpha)). $}\par No special properties of this matrix are assumed and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_eq {test functions} can be computed along the equilibrium curve to detect and process \jump test_eq {equilibrium singularities}. \par \par An equilibrium curve can be continued from a user-supplied \jump point_eq {point}, a \jump eq_eq {equilibrium point} (including \jump LP_eq {limit point} and \jump Hopf_eq {Hopf point}) or a \jump branch_eq {branching point} on the computed equilibrium curve. To \jump conti {continue} of the equilibrium curve, one needs to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_eqp {a point}\par $\;\;\;$ \jump start_eqe {an equilibrium}\par $\;\;\;$ \jump start_eqb {a branching point} \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_eq \topic start_eqp {starter parameters (equilibrium curve from a point)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_eq {equilibrium curve} from a user-supplied {\bf point} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 18 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Corrector data}\>\par {\li 1 MaxIter} \> Maximal number of corrections to locate the first point.\par {\li 1 MaxNewtonIter} \> Maximal number of corrections with recomputation of the Jacobian matrix.\par {\li 1 VarTolerance} \> Tolerance with respect to phase variables and parameters.\par {\li 1 FuncTolerance} \> Tolerance with respect to defining functions.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_eq {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the starter before computing the first point. Otherwise, CONTENT uses the values listed in the window. \par } \par \include ode_init \end \% start_eqp \topic start_eqe {starter parameters (equilibrium curve from an equilibrium)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_eq {equilibrium curve} from an \jump eq_eq {equilibrium} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 13 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_eq {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the starter before computing the first point. Otherwise, CONTENT uses the values listed in the window. \par } \par \include ode_init \end \% start_eqe \topic start_eqb {starter parameters (equilibrium curve from a branching point)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} an \jump ode_eq {equilibrium curve} from a {\bf branching point} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par { \desc 13 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted. \par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch switch data}\> \par {\li 1 Branch} \> Switch toggle ({\bf primary}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_eq {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par } \par \end \% start_eqb \topic test_eq {test functions and singularities along equilibrium curve} The following \jump kh_TestFunc {test-functions} can be computed along the \jump ode_eq {equilibrium curve} $F(y)=f(x,\alpha)=0$: \par \qc{ $ \table 3 rcl 0 0 {\psi_1}{\:=\:}{det \matrix 2 c 0 0 {F_y}{v^T},} {\psi_2}{\:=\:}{v_{n+1},} {\psi_3}{\:=\:}{det (2f_x \otimes I_n).} $} \par Here $v=(v_1,v_2,\ldots,v_{n+1})\in {\bf R}^{n+1}$ is the normalized tangent vector to the curve at a point $y=(x,\alpha)$ on it; $I_n$ is the unit $n\times n$ matrix; and $\otimes$ means the \jump bialt {bialternate product} of two matrices. \par \par The following {\bf singularities} can be detected and located as regular zeroes of the equilibrium test-functions: \par {\li 2 \desc 3 \jump branch_eq {Branching point}:\> $\psi_1=0.$ \par \jump LP_eq {Limit point}:\> $\psi_2=0,\; \psi_1 \neq 0.$ \par \jump Hopf_eq {Hopf bifurcation}:\> $\psi_3=0.$ \par } \end \% sing_eq \topic point_eq {point} You can compute \par {\li 2 \desc 3 \jump ode_eq {an equilibrium curve}\> \par \jump ode_or {an orbit}\> \par \jump ode_lc {a limit cycle}\> \par } \par starting from a {\bf point} in the phase-parameter space. In the former case, the point should be rather close to an \jump eq_eq {equilibrium}, while in the latter case it should be close to the limit cycle of known period. \end \topic ode_or {orbit of ODEs} Given a system of ODEs, \par \qc{ $ \dot{x}=F(x),\: \: x \in {\bf R}^n, $ } \par an {\bf orbit} starting at $x_0 \in {\bf R}^n$ is a curve $x=\varphi(t),\; -\infty\; <\; t\; < \; +\infty,\; \varphi(0)=x_0$, such that \par \qc{ $ \dot{\varphi}(t) = F(\varphi(t)). $ } \par The vector $F(x)$ is tangent to the orbit at each nonequilibrium point. \par Orbits are generated numerically by \jump integr {integrator}. \end \topic eq_eq {equilibrium} Given a system of ODEs, \par \qc{ $ \dot{x}=F(x),\: \: x \in {\bf R}^n, $ } \par an {\bf equilibrium} is a point $x_0 \in {\bf R}^n$ such that $F(x_{0})=0.$ \par {\bf Eigenvalues} $\lambda_{1},\lambda_{2},\ldots,\lambda_{n}$ of an equilibrium $x_0$ are roots of the {\bf characteristic equation} \par \qc{ $ h(\lambda) = det (A-\lambda I_n) = 0, $ } \par where \par \qc{$A=F_x(x_0)$} \par is the {\bf Jacobian matrix} and $I_n$ is the $n\; \times \;n$ unit matrix. \par The equilibrium is {\bf hyperbolic} if there is no eigenvalue with zero real part. A hyperbolic equilibrium has invariant {\bf stable and unstable manifolds} of dimension equal to the number of eigenvalues with negative and positive real part, respectively. If there is an eigenvalue with zero real part, the equilibrium is called {\bf nonhyperbolic}. A nonhyperbolic equilibrium has a {\bf center manifold} of dimension equal to the number of eigenvalues on the imaginary axis. The invariant manifolds are tangent to the corresponding generalized eigenspaces at $x_0$. \end %\eq_eq \topic branch_eq {branching point} Generically, two branches $M_{1,2}$ of the \jump ode_eq {equilibrium curve} intersect transversally at a {\bf branching point}: \par \qc{\picture 0.5 fig10_10.fig} \par Let $v \in {\bf R}^{n+1}$ be the normalized tangent vector to the primary (traced) branch of the equilibrium curve defined by \par \qc{$ F(y)=f(x,\alpha)=0. $} \par Consider the $(n+1)\;\times\; (n+1)$ matrix \par \qc{$ D=\matrix 2 c 0 0 {F_y(y_0)}{v^T} . $} \par A tangent vector $V$ to the secondary branch of the equilibrium curve at a branching point $y_0$ can be computed as following. A unit vector $q\in {\bf R}^{n+1}$ that is orthogonal to $v$ is a normalized null-eigenvector of $D$: \par \qc{$ Dq=0,\,\,\, \langle q,q\rangle = 1. $} \par A normalized null-eigenvector $p\in {\bf R}^{n+1}$ of the transposed matrix $D^T$, \par \qc{$ D^{T}p=0,\,\,\, \langle p,p\rangle = 1, $} \par has the form \par \qc{$ p=\matrix 2 c 0 0 {\phi}{0}, $} \par where $\phi \in {\bf R}^{n}$. A tangent vector $V$ to the secondary branch is then given by $V=q+kv$, where \par \qc{$ k=-\frac{b_{22}}{2b_{12}}=-\frac{\langle \phi,B(q,q)\rangle}{2\langle \phi,B(v,q)\rangle}, $} \par provided $b_{12}\neq 0$. Here the bilinear function $B(v,q)$ is given by \par \qc{$ B_i(v,q)= \Sum {j,k=1} {n+1} {{\brackets.| {\frac{\partial^2 F_i(y_0)}{\partial y_j\partial y_k}}}_{y=y_0}} \; v_j q_k, \, i=1,2,\ldots,n. $} \par \par The vector $V$ is then normalized and stored together with the tangent vector $v$. This allows to switch branches. If $b_{12}=0$, the branching point is {\bf non-simple} and the switching is not supported. The number $k=ctg\, \varphi$, where $\varphi$ is the angle between the branches, is displayed in the \jump kh_MainWin {message field}.\par \end \% branch_eq \topic LP_eq {limit point} At a {\bf limit point}, the tangent vector $v$ to the \jump ode_eq {equilibrium curve} has the form \par \qc{$v=\matrix 2 c 0 0 {q}{0}, $} \par where $q\in {\bf R}^{n}$ is a null-eigenvector, $Aq=0$, of the Jacobian matrix $A=f_x(x_0,\alpha_0)$. Generically, two equilibria collide and disappear at the critical parameter value $\alpha_0$. \par \par \qc{\picture 0.5 fig10_4.fig} \par \par The restricted to the one-dimensional \jump eq_eq {center manifold} system at $\alpha=\alpha_0$ has the form \par \qc{$ \dot \xi=\frac{1}{2}a \xi^2 + O(|\xi|^3),\; \xi \in {\bf R}^1. $}\par where the quadratic normal form coefficient $a$ is computed by\par \qc{$ a=\langle p,B(q,q) \rangle, $}\par where the eigenvectors $q$ and $p$ satisfy the equations\par \qc{$ Aq=0,\; A^T p=0,\; \langle q,q\rangle = 1,\; \langle p,q\rangle = 1, $}\par and\par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, \, i=1,2,\ldots,n. $}\par If $a \neq 0$ and the system depends generically on the parameter, two equilibria collide at the bifurcation point. \par The normalized eigenvector $q$ corresponding to zero eigenvalue is stored for two-parameter continuation of limit points. The value of $a$ can be read in the \jump kh_MainWin {message field}.\par \end \% LP_eq \topic Hopf_eq {Hopf bifurcation} If the critical equilibrium $x_0$ has a simple pair of purely imaginary eigenvalues $\lambda_{1,2}=\pm i\omega_0,\: \omega_0>0,$ then generically, a unique small \jump lc_lc {limit cycle} bifurcates from the equilibria as the parameter crosses the critical value $\alpha_0$ in one or another direction, while the equilibrium point changes stability. \par The restricted to the two-dimensional \jump eq_eq {center manifold} the \jump ode_Class {ODE system} has the complex form\par \qc{$ \dot z=i\omega_0 z+\frac{1}{2}g_{20}z^2 + g_{11}z\bar z + \frac{1}{2}g_{02}{\bar z}^2 + \frac{1}{2}g_{21}z^2\bar z + \ldots,\; z \in {\bf C}^1. $}\par This equation can be transformed by an invertible smooth change of coordinates to the Poincare normal form \par \qc{$ \dot z=i\omega_0 z + c_1 z^2\bar z + O(|z|^4), $}\par where\par \qc{$ Re\, c_1=\frac{1}{2\omega_0}Re(ig_{20}g_{11}+\omega_0 g_{21}). $}\par The {\bf first Lyapunov coefficient} $l_1=\frac{1}{\omega_0}Re\, c_1$ is computed at the critical parameter value by the formula\par \qc{$ l_1(0)=\frac{1}{2\omega_0} Re\brackets[. { \langle p,C(q,q,\bar q)\rangle - 2\langle p,B(q,A^{-1}B(q,\bar q))\rangle } \, + \, \brackets.] { \langle p,B(\bar q,{(2i\omega_0 I_n - A)}^{-1}B(q,q))\rangle }, $}\par where the complex eigenvectors $q,p \in {\bf C}^n$ satisfy\par \qc{$ Aq=i\omega_0 q,\; A^T p=-i\omega_0 p,\; \langle Re\,q,Im\, q\rangle = 0, \; \langle p,q\rangle = 1, $}\par the multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. \par Provided the critical eigenvalues cross the imaginary axis at a nonzero velocity with respect to the parameter, a unique stable \jump lc_lc {limit cycle} bifurcates if $l_1 <0$, \par \qc{\picture 0.9 supHopf.fig} \par while a unique unstable limit cycle exists for close parameter value if $l_1 >0$: \par \qc{\picture 0.9 subHopf.fig} \par If the critical equilibrium has a simple real eigenvalues $\lambda_{1,2}=\pm \omega_0,\: \omega_0>0,$ then it is a {\bf neutral saddle} which does not bifurcate. The equilibrium has two real linearly independent eigenvectors $q,p \in {\bf R}^n$ satisfy \par \qc{$ Aq=\omega_0 q,\; Ap=-\omega_0 p,\; \langle q,q,\rangle = \langle p,p\rangle = 1, $}\par \par The critical value $\omega_0$ and the normalized eigenvectors $p,q$ are stored for two-parameter continuation of the Hopf (neutral saddle) point, as well as for starting continuation of the bifurcating \jump ode_lc {limit cycle curve}. The starting periodic solution for the cycle continuation is taken in the form \par \qc{ $ x(t)=\varepsilon (q_{R} cos\:\omega_{0}t\: -\: q_{I} sin\:\omega_{0}t), $ } \par where $q_{R}+iq_{I}=\:q$ and $\epsilon$ is a user-defined small number. \par The values $\omega_0$ (and $l_1$) can be read in the \jump kh_MainWin {message field}.\par \end \topic ode_lp {limit point curve} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)} \par a {\bf limit point curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{N}$ endowed with coordinates $(x,y,\ldots,\alpha)$, whose projection onto the $(x,\alpha)$-space gives \jump LP_eq {limit points} of (1) having zero eigenvalue $\lambda_1=0.$ \par \par The \jump kh_DefFunc {defining functions} used in CONTENT for the limit point curve continuation can be specified by \par \par $\;\;\;$\jump ode_lp0 {the standard augmented system,} \par $\;\;\;$\jump ode_lp1 {the minimal bordered system.} \par \end %ode_lp \topic ode_lp0 {limit point curve (standard)} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)}\par a {\bf standard limit point curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{2n+2}$ endowed with coordinates $y=(x,v,\alpha)$ defined by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 rcl 1.2 1 {f(x,\alpha)}{=}{0} {f_{x}(x,\alpha)v}{=}{0} {\langle v,v\rangle - 1}{=}{0} $ }\qr {(2)}\par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_lp0 {test functions} can be computed along the limit point curve to detect and process \jump test_lp0 {limit point singularities}. \par \par A limit point curve can be continued from a regular \jump LP_eq {limit point} (including \jump cusp_eq {cusp}, \jump BT {Bogdanov-Takens}, and \jump ZH {zero-Hopf points}). To start \jump conti {continuation} of the limit point curve, one needs to set relevant \jump start_lp0 {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_lp0 \topic ode_lp1 {limit point curve (bordered)} This method to define the \jump ode_lp {limit point curve} is implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)}\par a {\bf boredered limit point curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+2}$ endowed with coordinates $y=(x,\alpha)$ defined by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 2 rcl 1.2 1 {f(x,\alpha)}{=}{0} {g(x,\alpha)}{=}{0,} $ }\qr {(2)} \par where the scalar function $g=g(x,\alpha)$ is obtained from solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W} {V^{T}}{0} \matrix 2 c 2 0 {v(x,\alpha)}{g(x,\alpha)}$ $=$ $\matrix 2 c 2 0 {0}{1}, $} \par where $A(x,\alpha)=f_x(x,\alpha)$ and the vectors $V,W \in {\bf R}^n$ are adapted along the curve to make the above system nonsingular. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_lp1 {test functions} can be computed along the limit point curve to detect and process \jump test_lp1 {limit point singularities}. \par \par A limit point curve can be continued from a regular \jump LP_eq {limit point} (including \jump cusp_eq {cusp}, \jump BT {Bogdanov-Takens}, and \jump ZH {zero-Hopf points}). To start \jump conti {continuation} of the limit point curve, one needs to set relevant \jump start_lp1 {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_lp1 \topic test_lp0 {test functions and singularities along the standard limit point curve} The \jump kh_TestFunc {test functions} that can be computed along the \jump ode_lp0 {standard limit point curve} are the following: \par \par \qc{ $ \table 3 rcl 0 0 {\psi_1}{\:=\:}{\langle w,B(v,v)\rangle,} {\psi_2}{\:=\:}{\langle \tilde{w},v\rangle,} {\psi_3}{\:=\:}{det (2f_x \otimes I_n).} $} \par Here $A=f_x(x,\alpha)$ and eigenvectors $v$ and $w(\tilde{w})$ satisfy the equations \par \qc{$ Av=0,\; A^T w=A^T \tilde{w}=0,\; \langle v,v\rangle = \langle \tilde{w},\tilde{w}\rangle = 1,\; \langle w,v\rangle = 1, $}\par and\par \qc{$ B_i(v,w)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; v_j w_k, \, i=1,2,\ldots,n; $}\par and $\otimes$ means the \jump bialt {bialternate product} of two matrices. \par \par The following {\bf singularities} can be detected and located as regular zeroes of the limit point test-functions: \par {\li 2 \desc 3 \jump cusp_eq {Cusp}:\> $\psi_1=0,\; \psi_3 \neq 0.$ \par \jump BT {Bogdanov-Takens}:\> $\psi_2=\psi_3=0.$ \par \jump ZH {Zero-Hopf}:\> $\psi_3=0,\; \psi_2 \neq 0.$ \par } \end \% test_lp0 \topic test_lp1 {test functions and singularities along the bordered limit point curve} The \jump kh_TestFunc {test functions} that can be computed along the \jump ode_lp1 {bordered limit point curve} are the following: \par \par \qc{ $ \table 3 rcl 0 0 {\psi_1}{\:=\:}{\langle w,B(v,v)\rangle,} {\psi_2}{\:=\:}{\langle \tilde{w},v\rangle,} {\psi_3}{\:=\:}{det (2f_x \otimes I_n).} $} \par Here $A=f_x(x,\alpha)$ and eigenvectors $v$ and $w(\tilde{w})$ satisfy the equations \par \qc{$ Av=0,\; A^T w=A^T \tilde{w}=0,\; \langle v,v\rangle = \langle \tilde{w},\tilde{w}\rangle = 1,\; \langle w,v\rangle = 1, $}\par and\par \qc{$ B_i(v,w)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; v_j w_k, \, i=1,2,\ldots,n; $}\par and $\otimes$ means the \jump bialt {bialternate product} of two matrices. \par \par The following {\bf singularities} can be detected and located as regular zeroes of the limit point test-functions: \par {\li 2 \desc 3 \jump cusp_eq {Cusp}:\> $\psi_1=0,\; \psi_3 \neq 0.$ \par \jump BT {Bogdanov-Takens}:\> $\psi_2=\psi_3=0.$ \par \jump ZH {Zero-Hopf}:\> $\psi_3=0,\; \psi_2 \neq 0.$ \par } \end \% test_lp0 \topic start_lp0 {starter parameters (standard limit point curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_lp0 {standard limit point curve} has the following {\bf parameters}: \par \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lp0 {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_lp0 \topic start_lp1 {starter parameters (limit point curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_lp1 {bordered limit point curve} has the following {\bf parameters}: \par \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lp1 {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_lp1 \topic ode_H {Hopf curve} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)} \par a {\bf Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{N}$ endowed with coordinates $(x,y,\ldots,\alpha)$, whose projection onto the $(x,\alpha)$-space gives \jump Hopf_eq {Hopf points} of (1) having a pair of \jump eq_eq {eigenvalues} with zero sum: $\lambda_1+\lambda_2=0.$ If the critical eigenvalues are purely imaginary $(\lambda_{1,2}=\pm i\omega_0,\; \omega_0>0)$, the equilibrium is a \jump eq_eq {nonhyperbolic} {\bf neutral focus}. If the eigenvalues are real $(\lambda_{1,2}=\pm \gamma_0,\; \gamma_0>0)$, the equilibrium is a {\bf neutral saddle}. \par \par The \jump kh_DefFunc {defining functions} used in CONTENT for Hopf curve continuation can be specified by \par \par $\;\;\;$\jump ode_H0 {the standard augmented system,} \par $\;\;\;$\jump ode_H1 {the bordered biproduct system.} \par $\;\;\;$\jump ode_H2 {the bordered squared Jacobian system.} \par \par Supported {\bf methods} of Hopf curve computation differ by \jump kh_DefFunc {defining}, \jump kh_TestFunc {test}, and \jump kh_ProcFunc {processing functions}. \end %H \topic ode_GH {generalized Hopf curve} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)} \par a {\bf generalized Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{N}$ endowed with coordinates $(x,y,\ldots,\alpha)$, whose projection onto the $(x,\alpha)$-space gives \jump GH {generalized Hopf (Bautin) points} of (1) having a pair of purely imaginary \jump eq_eq {eigenvalues} $(\lambda_{1,2}=\pm i\omega_0,\; \omega_0>0)$, and zero \jump Hopf_eq {first Lyapunov coefficient} $l_1=0$. \par \par The \jump kh_DefFunc {defining functions} used in CONTENT for the generalized Hopf continuation can be specified by \par \par \par $\;\;\;$\jump ode_GH1 {the minimally extended system.} \par $\;\;\;$\jump ode_GH2 {the maximally extended system.} \par \par Supported {\bf methods} of the generalized Hopf curve computation differ by \jump kh_DefFunc {defining functions}, \end %GH \topic ode_H0 {Hopf curve (standard)} Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)}\par a {\bf standard} \jump ode_H {Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{2n+3}$ endowed with coordinates $y=(x,v,\kappa,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 4 rcl 1.2 1 {f(x,\alpha)}{=}{0} {([f_{x}(x,\alpha)]^{2}+\kappa I)v}{=}{0} {\langle v,v\rangle - 1}{=}{0} {\langle w,v\rangle}{=}{0,} $ }\qr {(2)} \par where $w\in {\bf R}^n$ is a vector that is adapted along the curve to be {\bf not orthogonal} to the null-eigenspace of the matrix $[f_{x}(x,\alpha)]^{2}+\kappa I$. Namely, to compute the next point on the curve, vector $w$ is taken to be orthogonal to the vector $v$ at the found point and such that $\{v,w\}$ span the eigenspace $T^{c}$ corresponding to the zero-sum eigenpair. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par Several \jump test_H {test functions} can be computed along the standard Hopf curve to detect and process \jump test_H {Hopf singularities}. \par \par The standard Hopf curve can be continued from a regular \jump Hopf_eq {Hopf point} (including \jump BT {Bogdanov-Takens}, \jump GH {generalized Hopf}, and \jump ZH {zero-Hopf points}) or from a \jump DH {double Hopf point}. To \jump conti {continue} of the standard Hopf curve, one has to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_H0 {a Hopf point}\par $\;\;\;$ \jump start_HDH {a double Hopf point} \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_H0 \topic ode_H1 {Hopf curve (biproduct)} This method to define the \jump ode_H {Hopf curve} is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)}\par a {\bf biproduct} \jump ode_H {Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+2}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 2 rcl 1.2 1 {f(x,\alpha)}{=}{0} {det\; G(x,\alpha)}{=}{0,} $ }\qr {(2)}\par where the components $g_{ij}$ of the $(2\times 2)-$matrix $G$ are obtained from solving the double bordered $(m+2)-$dimensional system with $2m =n(n-1)$ \par \qc{$ \matrix 3 lrr 2 1 {B}{W_{1}}{W_{2}} {V^{T}_{1}}{d_{11}}{d_{12}} {V^{T}_{2}}{d_{21}}{d_{22}}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}. $} \par Here the matrix $B=2f_{x}(x,\alpha) \otimes I_{n}$ is a $(m\;\times\;m)-$matrix , where $\otimes$ stands for the \jump bialt {bialternate product} of two $(n\;\times\;n)$ matrices. The vectors $V_1,V_2,W_1,W_2\in {\bf R}^m$ and scalars $d_{11}, d_{12}, d_{21}, d_{22}$ are adapted along the curve to make the above system {\bf nonsingular}. The defining system (2) actually determines Jacobian matrices with a zero-sum pair of eigenvalues. So the solution branch may contain Hopf points, Bogdanov-Takens points and neutral saddle points. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par Several \jump test_H1 {test functions} can be computed along the biproduct Hopf curve to detect and process \jump test_H1 {Hopf singularities}. \par \par The biproduct Hopf curve can be continued from a regular \jump Hopf_eq {Hopf point} (including \jump BT {Bogdanov-Takens}, \jump GH {generalized Hopf}, and \jump ZH {zero-Hopf points}) or from a \jump DH {double Hopf point}. To \jump conti {continue} of the biproduct Hopf curve, one has to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_H1 {a Hopf point}\par $\;\;\;$ \jump start_H1DH {a double Hopf point} \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_H1 \topic ode_H2 {Hopf curve (bordered squared Jacobian)} This method to define the \jump ode_H {Hopf curve} is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^2, $ } \qr {(1)}\par a {\bf squared Jacobian} \jump ode_H {Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha,\kappa)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 rcl 1.2 1 {f(x,\alpha)}{=}{0} {g_{i_{1},j_{1}}(x,\alpha)}{=}{0,} {g_{i_{2},j_{2}}(x,\alpha)}{=}{0,} $ }\qr {(2)}\par where the components $g_{ij}$ of the $(2\times 2)-$matrix $G$ are obtained from solving the double bordered $(n+2)-$dimensional system \par \qc{$ \matrix 3 lrr 2 1 {B}{W_{1}}{W_{2}} {V^{T}_{1}}{0}{0} {V^{T}_{2}}{0}{0}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}. $} \par Here the matrix $B=[f_{x}(x,\alpha)]^2 + \kappa I_{n}$. The vectors $V_1,V_2,W_1,W_2\in {\bf R}^n$ are adapted along the curve to make the above system {\bf nonsingular}. The defining system (2) actually determines Jacobian matrices with a zero-sum pair of eigenvalues. So the solution branch may contain Hopf points, Bogdanov-Takens points and neutral saddle points. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par Several \jump test_H2 {test functions} can be computed along the biproduct Hopf curve to detect and process \jump test_H2 {Hopf singularities}. \par \par The bordered squared Hopf curve can be continued from a regular \jump Hopf_eq {Hopf point} (including \jump BT {Bogdanov-Takens}, \jump GH {generalized Hopf}, and \jump ZH {zero-Hopf points}) or from a \jump DH {double Hopf point}. To \jump conti {continue} of the bordered squared Hopf curve, one has to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_H2 {a Hopf point}\par $\;\;\;$ \jump start_H2DH {a double Hopf point} \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_H1 \topic test_H {test functions and singularities along the standard Hopf curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_H0 {standard Hopf curve} are the following: \par \qc{ $ \table 4 rcl 0 0 {\psi_1}{\:=\:}{\frac{1}{2\omega_0} Re\brackets[. { \langle p,C(q,q,\bar q)\rangle - 2\langle p,B(q,A^{-1}B(q,\bar q))\rangle } \, + \, \brackets.] { \langle p,B(\bar q,{(2i\omega_0 I_n - A)}^{-1}B(q,q))\rangle },} {\psi_2}{\:=\:}{\kappa,} {\psi_3}{\:=\:}{det(A),} {\psi_4}{\:=\:}{det(2A|_{T^{c}}\otimes I_{n-2}).} $} \par Here $A = f_x(x,\alpha)$, the complex vectors $q,p \in {\bf C}^n$ satisfy\par \qc{$ Aq=i\omega_0 q,\; A^T p=-i\omega_0 p,\; \langle Re\,q,Im\, q\rangle = 0, \; \langle p,q\rangle = 1, $}\par $T^{c} \subset {\bf R}^{n}$ is the critical two-dimensional eigenspace corresponding to the pair of eigenvalues of $A$ with zero sum, while the multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. The test function $\psi_1$ is only defined when the critical eigenvalues are {\bf imaginary}, i.e. $\kappa=\omega^{2}_{0} > 0$. \par The following {\bf singularities} can be detected and located as regular zeros of the above defined test-functions: \par {\li 2 \desc 4 \jump GH {Generalized Hopf}:\> $\psi_1=0.$ \par \jump BT {Bogdanov-Takens}:\> $\psi_2=\psi_3=0.$ \par \jump ZH {Zero-Hopf}:\> $\psi_3=0,\; \psi_2 \neq 0.$ \par \jump DH {Double Hopf}:\> $\psi_4=0.$ \par } \end \% test_H \topic test_H1 {test functions and singularities along the biproduct Hopf curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_H1 {biproduct Hopf curve} are the following: \par \qc{ $ \table 4 rcl 0 0 {\psi_1}{\:=\:}{\frac{1}{2\omega_0} Re\brackets[. { \langle p,C(q,q,\bar q)\rangle - 2\langle p,B(q,A^{-1}B(q,\bar q))\rangle } \, + \, \brackets.] { \langle p,B(\bar q,{(2i\omega_0 I_n - A)}^{-1}B(q,q))\rangle },} {\psi_2}{\:=\:}{\frac{\langle v,Av\rangle \; \langle w,Aw\rangle - \langle w,Av\rangle \; \langle v,Aw\rangle} {\langle v,v\rangle \; \langle w,w\rangle - (\langle v,w\rangle)^{2}},} {\psi_3}{\:=\:}{det(A),} {\psi_4}{\:=\:}{g_{22}.} $} \par Here $A = f_x(x,\alpha)$; two real vectors $v,w \in {\bf R}^n$ are such that $Q=v \wedge w \in {\bf R}^m$ is the singular vector of $2f_{x}(x,\alpha) \otimes I_{n}$ (where $\wedge$ denotes the \jump wedge {wedge product} of two vectors, and $\otimes$ stands for the \jump bialt {bialternate product} of two matrices); the complex vectors $q,p \in {\bf C}^n$ satisfy\par \qc{$ Aq=i\omega_0 q,\; A^T p=-i\omega_0 p,\; \langle Re\,q,Im\, q\rangle = 0, \; \langle p,q\rangle = 1. $}\par The multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. The test function $\psi_1$ is only defined when the critical eigenvalues are {\bf imaginary}, i.e $\psi_{2}=\omega^{2}_{0}>0$. \par The following {\bf singularities} can be detected and located as regular zeros of the above defined test-functions: \par {\li 2 \desc 4 \jump GH {Generalized Hopf}:\> $\psi_1=0.$ \par \jump BT {Bogdanov-Takens}:\> $\psi_2=\psi_3=0.$ \par \jump ZH {Zero-Hopf}:\> $\psi_3=0,\; \psi_2 \neq 0.$ \par \jump DH {Double Hopf}:\> $\psi_4=0.$ \par } \end \% test_H1 \topic test_H2 {test functions and singularities along the squared Hopf curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_H1 {squared Hopf curve} are the following: \par \qc{ $ \table 4 rcl 0 0 {\psi_1}{\:=\:}{\frac{1}{2\omega_0} Re\brackets[. { \langle p,C(q,q,\bar q)\rangle - 2\langle p,B(q,A^{-1}B(q,\bar q))\rangle } \, + \, \brackets.] { \langle p,B(\bar q,{(2i\omega_0 I_n - A)}^{-1}B(q,q))\rangle },} {\psi_2}{\:=\:}{\kappa,} {\psi_3}{\:=\:}{det(A),} {\psi_4}{\:=\:}{det(2A|_{T^{c}}\otimes I_{n-2}).} $} \par Here $A = f_x(x,\alpha)$, the complex vectors $q,p \in {\bf C}^n$ satisfy\par \qc{$ Aq=i\omega_0 q,\; A^T p=-i\omega_0 p,\; \langle Re\,q,Im\, q\rangle = 0, \; \langle p,q\rangle = 1, $}\par $T^{c} \subset {\bf R}^{n}$ is the critical two-dimensional eigenspace corresponding to the pair of eigenvalues of $A$ with zero sum, while the multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. The test function $\psi_1$ is only defined when the critical eigenvalues are {\bf imaginary}, i.e. $\kappa=\omega^{2}_{0} > 0$. \par The following {\bf singularities} can be detected and located as regular zeros of the above defined test-functions: \par {\li 2 \desc 4 \jump GH {Generalized Hopf}:\> $\psi_1=0.$ \par \jump BT {Bogdanov-Takens}:\> $\psi_2=\psi_3=0.$ \par \jump ZH {Zero-Hopf}:\> $\psi_3=0,\; \psi_2 \neq 0.$ \par \jump DH {Double Hopf}:\> $\psi_4=0.$ \par } \end \% test_H2 \topic start_H0 {starter parameters (standard Hopf curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H0 {standard Hopf curve} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_H0 \topic start_H1 {starter parameters (biproduct Hopf curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H1 {biproduct Hopf curve} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H1 {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_H1 \topic start_H2 {starter parameters (squared Hopf curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H2 {squared Hopf curve} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H2 {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_H2 \topic start_HDH {starter parameters (standard Hopf curve from DH)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H0 {standard Hopf curve} starting at a \jump DH {double Hopf point} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch switch data}\> \par {\li 1 Branch} \> Switch toggle ({\bf primary}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_HDH \topic start_H1DH {starter parameters (biproduct Hopf curve from DH)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H1 {biproduct Hopf curve} starting at a \jump DH {double Hopf point} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch switch data}\> \par {\li 1 Branch} \> Switch toggle ({\bf primary}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_H1DH \topic start_H2DH {starter parameters (squared Hopf curve from DH)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_H2 {squared Hopf curve} starting at a \jump DH {double Hopf point} has the following {\bf parameters}: \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch switch data}\> \par {\li 1 Branch} \> Switch toggle ({\bf primary}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_H2 {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_H2DH \topic ode_GH1 {Generalized Hopf curve (minimally extended)} This method to define the \jump ode_GH {generalized Hopf curve} is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf minimally extended generalized Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 rcl 1.2 1 {f(x,\alpha)}{=}{0} {det\; G(x,\alpha)}{=}{0,} {L_0}{=}{0,} $ }\qr {(2)}\par where the components $g_{ij}$ of the $(2\times 2)-$matrix $G$ are obtained from solving the double bordered $(m+2)-$dimensional system with $2m =n(n-1)$ \par \qc{$ \matrix 3 lrr 2 1 {B}{W_{1}}{W_{2}} {V^{T}_{1}}{d_{11}}{d_{12}} {V^{T}_{2}}{d_{21}}{d_{22}}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}. $} \par Here the matrix $B=2f_{x}(x,\alpha) \otimes I_{n}$ is a $(m\;\times\;m)-$matrix , where $\otimes$ stands for the \jump bialt {bialternate product} of two $(n\;\times\;n)$ matrices. The vectors $V_1,V_2,W_1,W_2\in {\bf R}^m$ and scalars $d_{11}, d_{12}, d_{21}, d_{22}$ are adapted along the curve to make the above system {\bf nonsingular}. Furthermore, $L_0$ is proportional to the first Lyapunov coefficient and is computed by the formula\par \qc{$ L_0=Re\brackets[. { \langle p,C(q,q,\bar q)\rangle - 2\langle p,B(q,A^{-1}B(q,\bar q))\rangle } \, + \, \brackets.] { \langle p,B(\bar q,{(2i\omega_0 I_n - A)}^{-1}B(q,q))\rangle }, $}\par where the complex eigenvectors $q,p \in {\bf C}^n$ satisfy\par \qc{$ Aq=i\omega_0 q,\; A^T p=-i\omega_0 p,\; \langle Re\,q,Im\, q\rangle = 0, \; \langle p,q\rangle = 1, $}\par and the multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par No test functions are computed along the generalized Hopf curve. \par \par The minimally extended generalized Hopf curve can be continued from a \jump GH {generalized Hopf point}. To \jump conti {continue} of the generalized Hopf curve, one has to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_GH1 {a generalized Hopf point}\par \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_GH1 \topic start_GH1 {starter parameters (generalized Hopf curve (minimal))} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_GH1 {minimally extended generalized Hopf curve} has the following {\bf parameters}: \par { \desc 13 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_GH1 \topic ode_GH2 {Generalized Hopf curve (maximally extended)} This method to define the \jump ode_GH {generalized Hopf curve} is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf maximally extended generalized Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{8n+6}$ endowed with coordinates $y=(x,\alpha,q,p,v,w,\omega_0,\lambda)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 8 rcl 1.2 1 {f(x,\alpha)}{=}{0} {Aq - i\omega_0 q}{=}{0,} {A^T p + \lambda p}{=}{0,} {\langle q_0, q \rangle -1}{=}{0,} {\langle p, q \rangle -1}{=}{0,} {Av - B(q,\bar{q})}{=}{0,} {(2i\omega_0 I_n - A)w - B(q,q)}{=}{0,} {Re\langle p,C(q,q,\bar{q})-2B(q,v)+B(\bar{q},w)\rangle}{=}{0.} $ }\qr {(2)}\par Here $A=f_{x}(x,\alpha)$ and $q_0$ is the eigenvector at a previously computed point at the curve. The vectors $q,q_0,p,w$ are complex. The complex variable $\lambda$ is introduced artificially to regularize the system; formally $\lambda=i\omega_0$ along the GH curve. The multilinear functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k, $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $}\par where $i=1,2,\ldots,n$. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par No test functions are computed along the generalized Hopf curve. \par \par The maximally extended generalized Hopf curve can be continued from a \jump GH {generalized Hopf point}. To \jump conti {continue} of the generalized Hopf curve, one has to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_GH2 {a generalized Hopf point}\par \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_GH1 \topic start_GH2 {starter parameters (generalized Hopf curve (maximal))} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_GH2 {maximally extended generalized Hopf curve} has the following {\bf parameters}: \par { \desc 13 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_GH2 \topic wedge {wedge product} The {\bf wedge product} $v \wedge w$ of two vectors $v,w \in {\bf R}^{n}$ is a vector $u \in {\bf R}^m$ with $m=\frac{n(n-1)}{2},$ whose components are enumerated by the multiindex $(p,q)$ and defined by the formula\par \qc{$ u_{(p,q)} = w_{q}v_{p} - w_{p}v_{q}, $} \par where $p=2,3,\ldots,\,n;$ $q=1,2,\ldots,\,p-1$. \par The wedge product $v \wedge w$ is linear with respect to both $v$ and $w$ and vanishes if and only if $v,w$ are linearly dependent. It is antisymmetric $(v \wedge w + w \wedge v = 0$) and is determined up to a scalar multiple by the two-dimensional space that contains $v,w$. Conversely, if it is nonzero, then it defines this space completely. \end \% wedge \topic bialt {bialternate product} The {\bf bialternate product} $(A \otimes B)$ of two $n\, \times\, n$ matrices $A$ and $B$ is an $m\, \times\, m$ matrix, $m=\frac{n(n-1)}{2},$ whose elements are enumerated by the multiindex $\{(p,q),(r,s)\}$, where $p=2,3,\ldots,\,n;$ $q=1,2,\ldots,\,p-1$, and $r=2,3,\ldots,\,n;$ $s=1,2,\ldots,\,r-1$, and are given by the formula \par \qc{ $ (A\otimes B)_{(p,q),(r,s)}=\frac{1}{2} {\brackets() {\brackets||\table 2 lr 0 0 {a_{pr}}{\;a_{ps}} {b_{qr}}{\:b_{qs}} \; + \; {\brackets||\table 2 cc 0 0 {b_{pr}}{\:b_{ps}} {a_{qr}}{\:a_{qs}}} }, $ } \par where $a_{jk}$ are elements of $A$, and $b_{jk}$ are elements of $B,\:j,k=1,2,\ldots,n$. There are two important special cases of the bialternate product:\par \qc{$ (2A\otimes I)_{(p,q),(r,s)}= \system 6 cl 1.2 1 {-a_{ps}} {{\f * if}\: r=q,} {a_{pr}} {{\f * if}\: r\neq p \: {\f * and}\: s=q,} {a_{pp}+\:a_{qq}} {{\f * if}\: r=p \: {\f * and}\: s=q,} {a_{qs}} {{\f * if}\: r=p \: {\f * and}\: s \neq q,} {-a_{qr}} {{\f * if}\: s=p,} {0} {\f * otherwise,} $} \par and \par \qc{ $ (A\otimes A)_{(p,q),(r,s)}={\brackets||\table 2 lr 0 0 {a_{pr}}{\:a_{ps}}{a_{qr}}{\:a_{qs}}}= a_{pr}a_{qs}-\:a_{ps}a_{qr}. $ } \par One has \par \qc{$ det (2A \otimes I_n) = \prod {j>k}{}{(\lambda_j+\:\lambda_k)}, $}\par and \par \qc{$ det (A \otimes A \:-\: I_m) = \prod {j>k}{}{(\lambda_{j}\lambda_{k} \: - \: 1)}, $}\par where $\lambda_k$ are the eigenvalues of matrix $A$. \end \% bialt \topic cusp_eq {cusp bifurcation} Suppose, at the critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the system has a simple zero eigenvalue and no other eigenvalues with zero real part, while the quadratic coefficient in the normal form of the \jump LP_eq {limit point} vanishes. At $\alpha=\alpha_0 \in {\bf R}^2$ the restricted to the one-dimensional \jump eq_eq {center manifold} generic system has the form \par \qc{$ \dot \xi=c \xi^3 + O(|\xi|^4),\; \xi \in {\bf R}^1. $} \par If $c\neq 0$ and the system depends generically on active parameters, its restriction to the parameter-dependent center manifold is locally topologically equivalent to the normal form \par \qc{$\dot{\xi}=\beta_1 + \beta_2 \xi + \sigma \xi^3,$} \par where $\sigma = sign(c)=\pm 1$. Therefore, generically, three equilibria collide at the cusp point and the projection of the \jump ode_lp {limit point curve} onto the parameter plane forms a {\bf semicubic parabola}. \par The {\bf cusp normal form coefficient} can be computed at the critical parameter value by the formula\par \qc{$ c=\frac{1}{6}\langle p,C(q,q,q)\rangle - \frac{1}{2} \langle p,B(q,s)\rangle, $}\par where the null-vectors $q,p \in {\bf R}^n$ are defined by \par \qc{$ Aq=0,\;\; A^T p=0, \;\; \langle q,q \rangle = \langle p,q \rangle = 1,\; $} \par while the vector $s \in {\bf R}^n$ can be extracted from the solution of the bordered $(n+1)$-dimensional system \par \qc{$ \matrix 2 lr 2 1 {A}{q}{p^{T}}{0} \matrix 2 c 2 0 {s}{u}$ $=$ $\matrix 2 c 2 0 {-\langle p, B(q,q) \rangle}{0}. $} \par Here the multilinear vector-functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ in terms of the system right-hand sides $f_i(x,\alpha),\; i=1,\ldots,n$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $} \par where $i=1,\ldots,n$. For relevant bifurcation diagrams, see \jump Ku95 {[Kuznetsov, 1995,1998]}. \par \par The null-eigenvector $q$ is stored for starting two-parameter continuation of the passing \jump ode_lp {limit point curve}, as well a for the three-parameter continuation of the \jump ode_cp {cusp bifurcation curve}. The value of the normal form coefficient $c$ can be read in the \jump kh_MainWin {message field}. \end %cusp_eq \topic ode_cp {cusp curve} This method to continue the \jump cusp_eq {cusp bifurcation} in three parameters is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf cusp curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 lcl 1.2 1 {f(x,\alpha)}{=}{0} {g_1(x,\alpha)}{=}{0} {g_2(x,\alpha)}{=}{0,} $ }\qr {(2)} \par where the scalar function $g_1=g_1(x,\alpha)$ is obtained from solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 2 0 {v_1(x,\alpha)}{g_1(x,\alpha)}$ $=$ $\matrix 2 c 3 0 {0}{1}, $} \par and the scalar function $g_2=g_2(x,\alpha)$ is obtained from solving another single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 2 0 {v_2(x,\alpha)}{g_2(x,\alpha)}$ $=$ $\matrix 2 c 2.5 0 {-B(v_1(x,\alpha),v_1(x,\alpha))}{0}. $} \par Here $A(x,\alpha)=f_x(x,\alpha)$ and the multilinear vector-function $B(q,p)$ is defined at $(x_0,\alpha_0)$ in terms of the system right-hand sides $f_i(x,\alpha),\; i=1,\ldots,n:$ \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k. $} \par The vectors $V_1,W_1 \in {\bf R}^n$ are adapted along the curve to make the above systems nonsingular. If $(x,\alpha)$ is a solution to (2), the system (1) has an equilibrium $x$ exibiting a \jump cusp_eq {cusp bifurcation} at the corresponding parameter values $\alpha$. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_cp {test functions} can be computed along the curve of cusp points to detect and process \jump test_cp {cusp singularities}. \par A cusp curve can be continued from a generic \jump cusp_eq {cusp point}, as well as from a \jump ZA {triple degenerate BT point}. To start \jump conti {continuation} of the cusp curve, one needs to set relevant \jump start_bt {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_cp \topic test_cp {test functions and singularities along the cusp curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_cp {cusp curve} are the following: \par \qc{ $ \table 2 rcl 0 0 {\psi_1}{\:=\:}{\frac{1}{6}\langle p,C(q,q,q)\rangle - \frac{1}{2} \langle p,B(q,s)\rangle,} {\psi_2}{\:=\:}{g_{\lambda}.} $} \par Here the null-vectors $q=v_1,p \in {\bf R}^n$ are defined by \par \qc{$ Aq=0,\;\; A^T p=0, \;\; \langle q,q \rangle = \langle p,q \rangle = 1,\; $} \par while the vector $s=v_2(x,\alpha) \in {\bf R}^n$. The multilinear vector-functions $B(q,p)$ and $C(p,q,r)$ are defined at $(x_0,\alpha_0)$ in terms of the system right-hand sides $f_i(x,\alpha),\; i=1,\ldots,n$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k $} \par and \par \qc{$ C_i(p,q,r)= \Sum {j,k,l=1} {n} {{\brackets.| { \frac{\partial^3 f_i(x,\alpha_0)}{\partial x_j\partial x_k\partial x_l}}}_{x=x_0}} p_j q_k r_l, $} \par where $i=1,\ldots,n$. \par The scalar function $g_{\lambda}=g_{\lambda}(x,\alpha)$ is obtained by solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 2 0 {v_{\lambda}(x,\alpha)} {g_{\lambda}(x,\alpha)}$ $=$ $\matrix 2 c 2.5 0 {v_1(x,\alpha)}{0}. $} \par The following {\bf singularities} can be detected and located along the \jump ode_cp {cusp curve} as regular zeros of the above defined test-functions: \par {\li 2 \desc 2 \jump ST {Swallow tail}:\> $\psi_1=0.$ \par \jump ZA {Triple degenerate BT}:\> $\psi_2=0$ \par } \end \% test_cp \topic start_cp {starter parameters (cusp curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_cp {cusp curve} has the following {\bf parameters}: \par \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_cp {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_cp \topic ST {Swallow tail bifurcation} The {\bf Swallow tail poiny} is a codim 3 point on the \jump ode_cp {cusp curve}, where the cubic coefficient of the normal form for the \jump ode_cp {cusp bifurcation} vanishes: \par \qc{$ c=0. $} \end \%ST \topic BT {Bogdanov-Takens bifurcation} Suppose, at the critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the system has an equilibrium with a double zero eigenvalue while the corresponding block of in the Jordan canonical form of the Jacobian matrix is nonzero. At $\alpha=\alpha_0 \in {\bf R}^2$ the restricted to the two-dimensional \jump eq_eq {center manifold} system is locally smoothly orbitally equivalent near the critical equilibrium $x_0$ to \par \qc{$ \system 2 rcl 1.2 1 {\dot{\xi}_1}{=}{\xi_2} {\dot{\xi}_2}{=}{a\xi^{2}_{1}+ b \xi_{1} \xi_{2} + O(||\xi||^3),} $} \par where $\xi = (\xi_{1},\xi_{2})^{T} \in {\bf R}^2$. If $ab\neq 0$ and the system depends generically on active parameters, its restriction to the parameter-dependent center manifold is locally topologically equivalent to the normal form \par \qc{$ \system 2 rcl 1.2 1 {\dot{\xi}_1}{=}{\xi_2} {\dot{\xi}_2}{=}{\beta_1 + \beta_{2} \xi_{1}+ \xi^{2}_{1} + s \xi_{1} \xi_{2}.} $} \par where $s = sign(ab)=\pm 1$. Therefore, generically, a Bogdanov-Takens point is a tangency point for a \jump ode_lp {limit point curve} and a \jump ode_H {Hopf curve}, as well as the origin of a {\bf homoclinic curve}. \par {\bf Bogdanov normal form coefficients} can be computed at the critical parameter values by the formulas\par \qc{$ a=\frac{1}{2}\langle r,B(q,q)\rangle,\;\;b=\langle s,B(q,q)\rangle + \langle r,B(q,p)\rangle, $}\par where the (generalized) eigenvectors $q,p,r,s \in {\bf R}^n$ are defined by\par \qc{$ Aq=0,\; Ap=q,\; A^T s=0,\; A^T r=s,\; $}\par and satisfy \qc{$ \langle q,q \rangle = \langle s,p \rangle = \langle r,q \rangle = 1,\; \langle q,p \rangle = \langle s,q \rangle = \langle r,p \rangle = 0,\; $} \par Here the multilinear function $B(q,p)$ is defined at $(x_0,\alpha_0)$ in terms of the system right-hand sides $f_i(x,\alpha),\; i=1,\ldots,n$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k. $} \par For relevant bifurcation diagrams, see \jump Ku95 {[Kuznetsov, 1995,1998]}. \par \par The generalized eigenvectors $q,p$ are stored for starting two-parameter continuation of the passing \jump ode_lp {limit point} and \jump ode_H {Hopf curve}, as well as for the three-parameter continuation of the \jump ode_bt {Bogdanov-Takens curve}. The value of the normal form coefficients $a,b$ can be read in the \jump kh_MainWin {message field}. \end %BT \topic ode_bt {Bogdanov-Takens curve} This method to continue the \jump BT {Bogdanov-Takens bifurcation} in three parameters is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf Bogdanov-Takens curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 lcl 1.2 1 {f(x,\alpha)}{=}{0} {g(x,\alpha)}{=}{0} {g_{\lambda}(x,\alpha)}{=}{0,} $ }\qr {(2)} \par where the scalar function $g=g(x,\alpha)$ is obtained from solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 3 0 {v(x,\alpha)}{g(x,\alpha)}$ $=$ $\matrix 2 c 3 0 {0}{1}, $} \par and the scalar function $g_{\lambda}=g_{\lambda}(x,\alpha)$ is obtained from solving another single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 2 0 {v_{\lambda}(x,\alpha)}{g_{\lambda}(x,\alpha)}$ $=$ $\matrix 2 c 3 0 {v(x,\alpha)}{0}. $} \par Here $A(x,\alpha)=f_x(x,\alpha)$ and the vectors $V_1,W_1 \in {\bf R}^n$ are adapted along the curve to make the above systems nonsingular. If $(x,\alpha)$ is a solution to (2), the system (1) has an equilibrium $x$ exibiting a \jump BT {Bogdanov-Takens bifurcation} at the corresponding parameter values $\alpha$. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_bt {test functions} can be computed along the BT-curve to detect and process \jump test_bt {BT singularities}. \par A BT-curve can be continued from a generic \jump BT {Bogdanov-Takens point}, nonlinearly degenerate Bogdanov-Takens points (i.e. \jump ZA {triple degenerate BT points} and \jump ZB {double degenerate BT points}), as well as from and \jump HBT {Hopf-BT points} and \jump TZ {triple zero points}. To start \jump conti {continuation} of the Bogdanov-Takens curve, one needs to set relevant \jump start_bt {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_bt \topic test_bt {test functions and singularities along the BT-curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_bt {Bogdanov-Takens curve} are the following: \par \qc{ $ \table 4 rcl 0 0 {\psi_1}{\:=\:}{g_{\lambda \lambda},} {\psi_2}{\:=\:}{\frac{1}{2}\langle r,B(q,q)\rangle,} {\psi_3}{\:=\:}{\langle s,B(q,q)\rangle + \langle r,B(q,p)\rangle,} {\psi_4}{\:=\:}{g_{22}.} $} \par Here the (generalized) eigenvectors $q=v,p=v_{\lambda},r,s \in {\bf R}^n$ are defined by \par \qc{$ Aq=0,\; Ap=q,\; A^T s=0,\; A^T r=s,\; $}\par and satisfy \par \qc{$ \langle q,q \rangle = \langle s,p \rangle = \langle r,q \rangle = 1,\; \langle q,p \rangle = \langle s,q \rangle = \langle r,p \rangle = 0.\; $} \par The multilinear function $B(q,p)$ is defined at $(x_0,\alpha_0)$ in terms of the system right-hand sides $f_i(x,\alpha),\; i=1,\ldots,n$ by \par \qc{$ B_i(q,p)= \Sum {j,k=1} {n} {{\brackets.| {\frac{\partial^2 f_i(x,\alpha_0)}{\partial x_j\partial x_k}}}_{x=x_0}} \; q_j p_k. $} \par The scalar function $g_{\lambda \lambda}=g_{\lambda \lambda}(x,\alpha)$ is obtained by solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_1} {V^{T}_1}{0} \matrix 2 c 2 0 {v_{\lambda \lambda}(x,\alpha)} {g_{\lambda \lambda}(x,\alpha)}$ $=$ $\matrix 2 c 2.5 0 {2v_{\lambda}(x,\alpha)}{0}, $} \par and the scalar function $g_{22}=g_{22}(x,\alpha)$ is obtained by solving the double bordered $(m+2)-$dimensional system with $2m=n(n-1)$: \par \qc{$ \matrix 3 lrr 2 1 {B(x,\alpha)}{W_{1}}{W_{2}} {V^{T}_{1}}{d_{11}}{d_{12}} {V^{T}_{2}}{d_{21}}{d_{22}}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}, $} \par with the matrix $B(x,\alpha)=2f_{x}(x,\alpha) \otimes I_{n}$ is a $(m\;\times\;m)-$matrix , where $\otimes$ stands for the \jump bialt {bialternate product} of two $(n\;\times\;n)$ matrices, while the vectors $V_1,V_2,W_1,W_2\in {\bf R}^m$ and scalars $d_{11}, d_{12}, d_{21}, d_{22}$ are selected to make the above system nonsingular. \par The following {\bf singularities} can be detected and located along the \jump ode_bt {Bogdanov-Takens curve} as regular zeros of the above defined test-functions: \par {\li 2 \desc 4 \jump TZ {Triple zero eigenvalue}:\> $\psi_1=0$ \par \jump ZA {Triple degenerate BT}:\> $\psi_2=0$ \par \jump ZB {Double degenerate BT}:\> $\psi_3=0$ \par \jump HBT {Hopf-BT}:\> $\psi_4=0.$ \par } \end \% test_bt \topic start_bt {starter parameters (Bogdanov-Takens bifurcation curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_bt {cusp curve} has the following {\bf parameters}: \par \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_bt {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_bt \topic ZA {triple degenerate BT point} The {\bf triple degenerate BT point} is a codim 3 point on the \jump ode_bt {Bogdanov-Takens curve}, where the coefficients of the normal form for the \jump BT {Bogdanov-Takens singularity} satisfy: \par \qc{$ a=0, \;\;\; b \neq 0. $} \end \%ZA \topic ZB {double degenerate BT point} The {\bf double degenerate BT point} is a a codim 3 point on the \jump ode_bt {Bogdanov-Takens curve}, where the coefficients of the normal form for the \jump BT {Bogdanov-Takens singularity} satisfy: \par \qc{$ a \neq 0, \;\;\; b = 0. $} \end \%ZB \topic TZ {triple zero point} The {\bf triple zero point} is a a codim 3 point on the \jump ode_bt {Bogdanov-Takens curve}, where the corresponding critical equilibrium has a triple zero eigenvalue. \end \% TZ \topic HBT {Hopf-BT point} The {\bf Hopf-BT point} is a a codim 3 point on the \jump ode_bt {Bogdanov-Takens curve}, where the corresponding critical equilibrium has a zero-sum pair of eigenvalues and a double zero eigenvalue. \end \%HBT \topic ZH {zero-Hopf bifurcation} Suppose, at the critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the system has a simple pair of complex eigenvalues with zero real part, as well as simple zero eigenvalue. Near these critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the restricted to the three-dimensional \jump eq_eq {center manifold} generic system is locally smoothly orbitally equivalent to \par \qc{$ \system 2 rcl 1.2 1 {\dot{u}}{=}{\beta_{1}+Bu^{2}+C|z|^{2}+O(||(u,z,\bar{z})||^{4})} {\dot{z}}{=}{(\beta_{2}+i\omega)z + D uz + Eu^{2}z + O(||(u,z,\bar{z})||^{4}),} $} \par where $u \in {\bf R}^{1}, z \in {\bf C}^{1},$ while $\omega,B,C,E \in {\bf R}^1$ and $D \in {\bf C}^1$ are smooth functions of $\beta \in {\bf R}^2.$ Generically, a zero-Hopf point is a tangency point for a \jump ode_lp {limit point curve} and a \jump ode_H {Hopf curve}. "Chaotic" dynamics may also be present for nearby parameter values. For computation of the normal form coefficients $(\omega,B,C,D,E)$ and for approximate bifurcation diagrams, see \jump Ku95 {[Kuznetsov, 1995,1998]}. \par \par The Hopf frequency $\omega_0$ and the critical eigenvectors are stored for starting two-parameter continuation of the passing \jump ode_lp {limit point}, \jump ode_H {Hopf}, as well as \jump ode_zh {zero-Hopf curves}. The value of $\omega_0$ can be read in the \jump kh_MainWin {message field}. \end %ZH \topic ode_zh {zero-Hopf curve} This method to continue the \jump ZH {zero-Hopf bifurcation} in three parameters is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf zero-Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 lcl 1.2 1 {f(x,\alpha)}{=}{0} {g(x,\alpha)}{=}{0} {g_11(x,\alpha)g_22(x,\alpha)-g_12(x,\alpha)g_21(x,\alpha)}{=}{0,} $ }\qr {(2)} \par where the scalar function $g=g(x,\alpha)$ is obtained from solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_{11}} {V^{T}_{11}}{0} \matrix 2 c 3 0 {v(x,\alpha)}{g(x,\alpha)}$ $=$ $\matrix 2 c 3 0 {0}{1}, $} \par and the functions $g_{ij}=g_{ij}(x,\alpha)$ are obtained by solving the double bordered $(m+2)-$dimensional system with $2m=n(n-1)$: \par \qc{$ \matrix 3 lrr 2 1 {B(x,\alpha)}{W_{1}}{W_{2}} {V^{T}_{1}}{d_{11}}{d_{12}} {V^{T}_{2}}{d_{21}}{d_{22}}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}, $} \par with the matrix $B(x,\alpha)=2f_{x}(x,\alpha) \otimes I_{n}$ is a $(m\;\times\;m)-$matrix , where $\otimes$ stands for the \jump bialt {bialternate product} of two $(n\;\times\;n)$ matrices. The vectors $V_{11},W_{11} \in {\Bf R}^n$, $V_1,V_2,W_1,W_2\in {\bf R}^m$ and the scalars $d_{11}, d_{12}, d_{21}, d_{22}$ are selected to make the above systems nonsingular. If $(x,\alpha)$ is a solution to (2), the system (1) has an equilibrium $x$ exibiting a \jump ZH {zero-Hopf bifurcation} at the corresponding parameter values $\alpha$. \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear system of equations. \par Several \jump test_zh {test functions} can be computed along the BT-curve to detect and process \jump test_zh {zero-Hopf singularities}. \par A zero-Hopf curve can be continued from generic \jump ZH {zero-Hopf points}, \jump HBT {Hopf-BT points} and \jump TZ {triple zero points}. To start \jump conti {continuation} of the zero-Hopf curve, one needs to set relevant \jump start_zh {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_zh \topic test_zh {test functions and singularities along the zero-Hopf curve} The \jump kh_TestFunc {test-functions} that can be computed along the \jump ode_zh {zero-Hopf curve} are the following: \par \qc{ $ \table 2 rcl 0 0 {\psi_1}{\:=\:}{g_{22},} {\psi_2}{\:=\:}{g_{\lambda},} $} \par the scalar function $g_{\lambda}=g_{\lambda}(x,\alpha)$ is obtained by solving the single bordered $(n+1)-$dimensional system \par \qc{$ \matrix 2 cc 2 1 {A(x,\alpha)}{W_{11}} {V^{T}_{11}}{0} \matrix 2 c 2 0 {v_{\lambda}(x,\alpha)} {g_{\lambda}(x,\alpha)}$ $=$ $\matrix 2 c 3 0 {v(x,\alpha)}{0}. $} \par The following {\bf singularities} can be detected and located along the \jump ode_zh {zero-Hopf curve} as regular zeros of the above defined test-functions: \par {\li 2 \desc 3 \jump ZDH {Zero-double Hopf}:\> $\psi_1=0,\; \psi_2\neq 0$ \par \jump HBT {Hopf-BT}:\> $\psi_1=\psi_2=0$ \par \jump TZ {Triple zero eigenvalue}:\> $\psi_2=0,\; \psi_1\neq 0.$ \par } \end \% test_zh \topic ZDH {zero-double Hopf point} The {\bf zero-double Hopf point} is a a codim 3 point on the \jump ode_zh {zero-Hopf curve}, where the corresponding critical equilibrium has two pairs of zero-sum eigenvalues and a zero eigenvalue. \end \%ZDH \topic start_zh {starter parameters (zero-Hopf bifurcation curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_zh {zero-Hopf curve} has the following {\bf parameters}: \par \par { \desc 15 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Singular Solver data}\>\par {\li 1 SingTolerance} \> Tolerance for the singularity of a matrix.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_zh {Singularity types}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_zh \topic GH {generalized Hopf bifurcation} Suppose, at the critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the system has a simple pair of complex eigenvalues with zero real part, as well as zero \jump Hopf_eq {first Lyapunov coefficient}. At $\alpha=\alpha_0 \in {\bf R}^2$ the restricted to the two-dimensional \jump eq_eq {center manifold} system is locally smoothly orbitally equivalent to \par \qc{$ \dot{z} = i\omega z + L_{2}z|z|^4 + O(|z|^6),\; z \in {\bf C}^1, $} \par where $L_{2} \in {\bf R}^1$. If $L_{2} \neq 0$ and the system depends generically on active parameters, its restriction to the parameter-dependent center manifold is locally topologically equivalent to the normal form \par \qc{$ \dot{z}=(\beta_1+i)z + \beta_{2} z|z|^2 + s z|z|^4, $} \par where $s = sign(L_{2})=\pm 1$. Therefore, generically, a generalized Hopf point is the origin of a {\bf cycle limit point curve}, which is tangent to the \jump ode_H {Hopf curve} passing through this point. For computation of the normal form coefficient $L_{2}$ and for relevant bifurcation diagrams, see \jump Ku95 {[Kuznetsov, 1995,1998]}. \par \par The Hopf frequency $\omega_0$ and the critical eigenvector are stored for starting two-parameter continuation of the passing \jump ode_H {Hopf curve}. The value of $\omega_0$ can be read in the \jump kh_MainWin {message field}. \end %GH \topic DH {double Hopf bifurcation} Suppose, at the critical parameter values $\alpha=\alpha_0 \in {\bf R}^2$ the system has two simple pairs of eigenvalues with zero real part. Near these critical parameter values, the restricted to the four-dimensional \jump eq_eq {center manifold} system is locally smoothly orbitally equivalent to \par \qc{$ \system 4 rcl 1.2 1 {\dot{r}_{1}}{=}{r_{1}(\mu_{1}+\;p_{11}r^{2}_{1}+\;p_{12}r^{2}_{2}+\;s_{1}r^{4}_{2}+O((r^2_1+r^2_2)^2)} {\dot{r}_2}{=}{r_{2}(\mu_2+\;p_{21}r^{2}_{1}+\;p_{22}r^{2}_{2}+\;s_{2}r^{4}_{1}+O((r^2_1+r^2_2)^2)} {\dot{\varphi}_1}{=}{\omega_1+\;O(1)} {\dot{\varphi}_2}{=}{\omega_2+\;O(1).} $} \par where $\omega_{i},\;p_{ij},\;s_i,\; i,j=1,2$ are smooth real functions of small parameters $\mu_{1}$ and $\mu_{2}$. Generically, a double Hopf point is the intersection of two \jump ode_H {Hopf curves}. Quasiperiodic and "chaotic" dynamics may also be present for nearby parameter values. For computation of the normal form coefficients and for approximate bifurcation diagrams, see \jump Ku95 {[Kuznetsov, 1995,1998]}. \par \par The Hopf frequencies $\omega_{1,2}$ and the critical eigenvectors are stored for starting two-parameter continuation of the passing \jump ode_H {Hopf curves}, as well as to start the continuation of the \jump ode_DH {double Hopf curve} in three active parameters. The values of $\omega_{1,2}$ can be read in the \jump kh_MainWin {message field}. \end %DH \topic ode_DH {double Hopf curve} The continuation of this curve is developed and implemented in cooperation with {\bf W. Govaerts} and {\bf B. Sijnave} (Vakgroep Toegepaste Wiskunde en Informatica, Universiteit Gent, Belgium). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot x=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^3, $ } \qr {(1)}\par a {\bf double Hopf curve} is a one-dimensional manifold $\Gamma$ in ${\bf R}^{n+3}$ endowed with coordinates $y=(x,\alpha)$ specified by the following \jump kh_DefFunc {defining functions} \par \qc{ $ \system 3 rcl 1.2 1 {f(x,\alpha)}{=}{0} {g_{i_{1}j_{1}}(x,\alpha)}{=}{0,} {g_{i_{2}j_{2}}(x,\alpha)}{=}{0,} $ }\qr {(2)}\par where the scalar functions $g_{ij},\; i=1,2,$ are extracted from the solution of the bordered $(m+2)$ dimensional system with $2m=n(n-1)$ \par \qc{$ \matrix 3 lrr 2 1 {B}{W_{1}}{W_{2}} {V^{T}_{1}}{d_{11}}{d_{12}} {V^{T}_{2}}{d_{21}}{d_{22}}\matrix 3 cc 2 0 {Q_{1}}{\;Q_{2}} {g_{11}}{\;g_{12}} {g_{21}}{\;g_{22}}$ $=$ $\matrix 3 cc 3 0 {0}{\;0} {1}{\;0} {0}{\;1}. $} \par Here $(m\;\times\;m)$ matrix $B=2f_{x}(x,\alpha) \otimes I_{n}$, where $\otimes$ stands for the \jump bialt {bialternate product} of two $(n\;\times\;n)$ matrices. The vectors $V_{i},W_{i}\in {\bf R}^m$ and scalars $d_{11}, d_{12}, d_{21}, d_{22}$ are adapted along the curve to make the above system {\bf nonsingular}. Also, the index pairs $(i_1,j_1)$ and $(i_2,j_2)$ are adapted along the double Hopf curve. \par The projection of the curve $\Gamma$ onto the $(x,\alpha)$-space defines equilibria of (1) having two pairs of eigenvalues with zero sum: $\lambda_1+\lambda_2=0,\;\lambda_3+\lambda_4=0.$ Such pairs can be Hopf pairs ($\lambda_i$ purely imaginary), Bogdanov-Takens pairs ($\lambda$ = 0) or neutral saddle pairs ($\lambda_i$ real). \par No special properties of the Jacobian matrix of (2) are assumed in the \jump conti {continuation} and the {\bf generic linear algebra library} is used to solve appearing linear systems. \par \par To start the \jump conti {continuation} of a double Hopf curve from a \jump DH {double Hopf point}, one has to set relevant \jump start_DH {starter parameters} via the \jump kh_StarterWin {starter window}. \end \% ode_DH \topic start_DH {starter parameters (double Hopf curve)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_DH {double Hopf curve} has the following {\bf parameters}: \par { \desc 11 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Eigenvalues}\> \par {\li 1 Compute} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_DH \topic integr {integrator} An \jump ode_or {orbit} starting at $t=t_0$ from a point $x_0 \in {\bf R}^n$ of a system of ODEs is approximated by the sequence of points $x^k,k=0,1,2,\ldots$ corresponding to time moments $t^k,k=0,1,2,\ldots$ using the following methods: \par $\;\;\;$ \jump method_EU {Euler}\par $\;\;\;$ \jump method_RK {Runge-Kutta}\par $\;\;\;$ \jump method_RD5 {RADAU5}\par \par Scalar functions $g_m(x,t)$ can be {\bf monitored} along the orbit. If such a function changes sign between $(x^{k},t^{k})$ and $(x^{k+1},t^{k+1})$, its zero can be {\bf detected} and computed by bi-sections to withing given tolerance. \par \par Parameters of the integrators can be modified via the corresponding \jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be set in the corresponding \jump kh_StarterWin {starter window}. \end \topic method_EU {Euler method} An \jump ode_or {orbit} starting at $t=t_0$ from a point $x_0 \in {\bf R}^n$ of a system of ODEs \par \qc{ $\dot{x}=F(t,x),\;\; x \in {\bf R}^n,$ } \par is approximated in the {\bf explicit Euler method with constant stepsize} by the sequence of points $x^k,k=0,1,2,\ldots$ corresponding to time moments $t^k,k=0,1,2,\ldots$ as following. \par {\bf Set}: $x^0=x_0,\; t^0=t_0$ \par {\bf Iterate for} $k=0,1,2,\ldots$ \par {\li 4 {$x^{k+1}=x^{k}+F(t^{k},x^{k})H_0,\;\; t^{k+1}=t^{k}+H_0$} } \par The \jump integrpar_EU {parameters of the integrator} can be modified via the corresponding \jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be set in the corresponding \jump kh_StarterWin {starter window}. \end \topic method_RK {Runge-Kutta method} An \jump ode_or {orbit} starting at $t=t_0$ from a point $x_0 \in {\bf R}^n$ of a system of ODEs \par \qc{ $\dot{x}=F(t,x),\;\; x \in {\bf R}^n,$ } \par is approximated in the explicit order 4 {\bf Runge-Kutta method with Merson's modification} by the sequence of points $x^k,k=0,1,2,\ldots,$ corresponding to time moments $t^k,k=0,1,2,\ldots$ as following. \par Introduce the matrix and the vector \par \qc{$ a=\matrix 4 rrrr 2 1 {\frac{1}{3}}{0}{0}{0} {\frac{1}{6}}{\frac{1}{6}}{0}{0} {\frac{1}{8}}{0}{\frac{3}{8}}{0} {\frac{1}{2}}{0}{-\frac{3}{2}}{2},\;\;\; c=\matrix 4 c 2 1 {0}{\frac{1}{3}}{0}{\frac{1}{2}} $} \par (their elements are $a_{ij},c_j,$ where $i=1,2,3,4;\; j=0,1,2,3$). \par {\bf Set}: $x^0=x_0,\; t^0=t_0,\; h^0=H_0$ \par {\bf Iterate for} $k=0,1,2,\ldots$\par AGAIN: {\li 2 \par {\bf If} $h^k < H_{min}$ {\bf or} $|t^k-t_0|>T_{max}$ {\bf then Stop}\par {\li 4 \par {$X_0=x^k$} \par {\bf Iterate for} $i=1,2,3,4:\; X_i=x^k+h^k\Sum {j=0} {j\leq i} {a_{ij} F(t^k+c_{j} h^k,X_j)},$ \par {$X_5=x^k+[F(t^k,x^k)+4F(t^k+c_{3} h^{k},X_3)+F(t^k+c_{4}h^k,X_4)]\frac{h^k}{6},$} \par {$E=||x_4-x_5||$} } \par {\bf If} $E > \varepsilon$ {\bf then} \par {\li 4 {$h^{k+1}=\frac{h^k}{2}$ {\bf goto} AGAIN} } \par {\bf else} \par {\li 4 {$x^{k+1}=X_4,\;\; t^{k+1}=t^k + h^k$} } \par {\bf endif} \par {\bf If} $E<\frac{\varepsilon}{50}$ {\bf then} $h^{k+1}=max(2h^k,H_{max})$ } \par \par The \jump integrpar_RK {parameters of the integrator} can be modified via the corresponding \jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be set in the corresponding \jump kh_StarterWin {starter window}. \end \topic method_RD5 {Radau IIA method} An \jump ode_or {orbit} starting at $t=t_0$ from a point $x_0 \in {\bf R}^n$ of a system of ODEs \par \qc{ $\dot{x}=F(t,x),\;\; x \in {\bf R}^n,$ } \par is approximated in the implicit order 5 {\bf Radau IIA method by Ehle} as a sequence of points $x^k,k=0,1,2,\ldots,$ corresponding to time moments $t^k,k=0,1,2,\ldots$, constructed as follows. \par Introduce one matrix and two vectors: \par \qc{$ a=\matrix 3 ccc 2 1 {\frac{88-7\sqrt{6}}{360}}{\frac{296-169\sqrt{6}}{1800}}{\frac{-2+3\sqrt{6}}{225}} {\frac{296+169\sqrt{6}}{1800}}{\frac{88+7\sqrt{6}}{360}}{\frac{-2-3\sqrt{6}}{225}} {\frac{16-\sqrt{6}}{36}}{\frac{16+\sqrt{6}}{36}}{\frac{1}{9}},\;\;\; b=\matrix 3 c 2 1 {\frac{16-\sqrt{6}}{36}}{\frac{16+\sqrt{6}}{36}}{\frac{1}{9}},\;\;\; c=\matrix 3 c 2 1 {\frac{4-\sqrt{6}}{10}}{\frac{4+\sqrt{6}}{10}}{1}. $} \par Denote by $\omega_{ij}$ the elements of the inverse matrix $a^{-1}$. \par {\bf Set}: $x^0=x_0,\; t^0=t_0,\; h^0=H_0$ \par {\bf Iterate for} $k=0,1,2,\ldots$\par AGAIN: {\li 2 \par {\bf If} $h^k < H_{min}$ {\bf or} $|t^k-t_0|>T_{max}$ {\bf then Stop}\par {\li 4 \par {\bf Iterate for} $i=1,2,3:\;$ \par {\bf Solve by Newton iterations for $X_i\;:$} $X_i=x^k+h^k\Sum {j=1} {3} {a_{ij} F(t^k+c_{j} h^k,X_j)},$ \par {\bf If} no convergence {\bf then} $h^{k}:=\frac{h^k}{2}$ {\bf goto} AGAIN \par $x^{k+1}=\brackets() { 1-\Sum {i,j=1} {3} {b_i \omega_{ij}} } x^k+ \Sum {i,j=1} {3} {b_i \omega_{ij} X_j},$ } \par {\bf Compute local error and new step size} $h^{k+1}$ \par $t^{k+1}=t^k + h^k$ } \par Exact step-control mechanism is described in Section IV.8 of the book E. Hairer and G. Wanner "Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems". Springer Series in Computational Mathematics 14, Springer-Verlag 1991. \par The \jump integrpar_RD5 {parameters of the integrator} can be modified via the corresponding \jump kh_GeneratorWin {generator window}. \jump start_or {Initial data} can be set in the corresponding \jump kh_StarterWin {starter window}. \end \topic integrpar_EU {integrator parameters (Euler method)} The \jump method_EU {Euler integrator} has the following parameters: \par \par \desc 6 {\bf Integration data}\>\par {\li 1 Interval} \> Integration interval length $T_{max}$\par {\li 1 InitStepsize} \> Initial stepsize $H_0$\par {\bf User function zero}\>\par {\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par {\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par \par The parameters can be modified via the \jump kh_GeneratorWin {generator window}. \end \topic integrpar_RK {integrator parameters (Runge-Kutta method)} The \jump method_RK {Runge-Kutta integrator} has the following parameters: \par \par \desc 10 {\bf Integration data}\>\par {\li 1 Interval} \> Integration interval length $T_{max}$\par {\li 1 InitStepsize} \> Initial stepsize $H_0$\par {\li 1 CurrentStepsize} \> Size of the current step along the orbit\par {\li 1 MinStepsize} \> Minimal stepsize $H_{min}$\par {\li 1 MaxStepsize} \> Maximal stepsize $H_{max}$\par {\li 1 Tolerance} \> Absolute step tolerance $\varepsilon$\par {\bf User function zero}\>\par {\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par {\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par \par The parameters can be modified via the \jump kh_GeneratorWin {generator window}. \end \topic integrpar_RD5 {integrator parameters (RADAU5)} The \jump method_RD5 {RADAU5 integrator} has the following parameters: \par \par \desc 11 {\bf Integration data}\>\par {\li 1 Interval} \> Integration interval length $T_{max}$\par {\li 1 InitStepsize} \> Initial stepsize $H_0$\par {\li 1 MinStepsize} \> Minimal stepsize $H_{min}$\par {\li 1 MaxStepsize} \> Maximal stepsize $H_{max}$\par {\li 1 Tolerance} \> Absolute step tolerance $\varepsilon_a$\par {\li 1 RelTolerance} \> Relative step tolerance $\varepsilon_r$\par {\li 1 MaxNewtonIter} \> Maximum number of Newton iterations $M_{Newt}$\par {\bf User function zero}\>\par {\li 1 UserFuncTolerance} \> Accuracy of zero location of user functions\par {\li 1 MaxIterFunc} \> Maximal number of bi-sections for zero location\par \par The parameters can be modified via the \jump kh_GeneratorWin {generator window}. \end \topic start_or {starter parameters for orbit} The \jump kh_Starter {starter} which sets necessary data to \jump kh_Generator {generate} the \jump ode_or {orbit} has the following {\bf parameters}: \par { \desc 8 {\bf Initial point}\>\par {\li 1 $t$} \> A name of the time variable and its initial value.\par {\li 1 $x$} \> A list of phase coordinate names and their values at the initial point.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point.\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function names}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}).\par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables and parameters. When it is specified and activated, it is called by the started before computing the first point. \par } \par The parameters can be modified via the \jump kh_StarterWin {starter window}. \end \% start_or \topic Ku95 {book on bifurcations} Kuznetsov, Yu.A. "Elements of Applied Bifurcation Theory", Springer-Verlag, New York, 1995,1998. \end \topic ode_init {-} Here is an example of a function SetInitPoint: \par \verbatim 2 X=A; Y=A*A; \par In this example, X and Y are phase variables and A is a parameter, as specified in the \jump kh_RHSWin {system specification dialog box.} \par In general, you write a body of a C function that assigns values to the coordinates of the initial point and/or the corresponding parameter values. \end \topic ode_lc {limit cycle curve} The continuation of this curve is developed and implemented in cooperation with E. Doedel (Computer Science, Concordia University, Montreal, Canada). \par Given a \jump ode_Class {system of ODEs}, \par \qc{ $ \dot{x}=f(x,\alpha),\: \: x \in {\bf R}^n,\: \alpha \in {\bf R}^1, $ } \qr {(1)} \par a {\bf limit cycle curve} is a one-dimensional manifold in the space of smooth vector functions on the unit interval, satisfying the following boundary-value problem with the periodic boundary conditions and an integral phase condition: \par \qc{ $ \system 3 rcl 1.2 1 {\dot{u}(\tau)-T_{0}f(u(\tau),\:\alpha)}{=}{0,\:\: \tau \in (0,1)} {u(0)-u(1)}{=}{0,} {\Int {0} {1} {\langle u(\tau),\dot{v}(\tau)\rangle \: d\tau}}{=}{0.} $ }\qr {(2)} \par Each solution of (2) corresponds to a \jump lc_lc {limit cycle} of period $T_0$ at parameter value $\alpha$ in (1), that has the minimal $L_2$-distance from a reference solution $v(\tau),$ \par \qc{ $ \rho(\sigma)=\Int {0} {1} ||u(\tau+\sigma)-v(\tau)||^2\: d\tau, $ } \par over possible phase shifts $\sigma$. \par The boundary-value problem (1) is approximated in CONTENT by orthogonal collocation as following. Introduce a partitioning of the interval $[0,1]$ by $N-1$ internal mesh points ($N \geq 2$), \par \qc{ $ 0=\tau_0<\tau_1<\ldots<\tau_{N-1}<\tau_N=1, $ } \par and approximate the solution to (2) by a piecewise-differentiable continuous function that is a vector-polynomial $u^{(j)}(\tau)$ of degree $m$ within each subinterval $[\tau_{j},\tau_{j+1}], j=0,1,\ldots,N-1$. "Collocation" consists in requiring the approximate solution to satisfy exactly the first equation in (2) at $m$ ($2 \leq m \leq 7$) {\bf collocation points} within each subinterval: \par \qc{ $ \tau_{j}<\zeta_{j,1}<\zeta_{j,2}<\ldots<\zeta_{j,m}<\tau_{j+1}. $ } \par Namely, \par \qc{ $ {\brackets.|\frac{du^{(j)}}{d\tau}}_{\tau=\zeta_{j,i}}-T_{0}f(u^{(j)}(\zeta_{j,i}))=0, $ }\qr{(3)} \par for $i=1,\ldots,m$, whenever $j=0,1,\ldots,N-1$. \par Characterize each polynomial $u^{(j)}(\tau),\: j=0,1,\ldots,N-1,$ by vectors \par \qc{ $ u^{j,k}=u^{(j)}(\tau_{j,i})\in {\bf R}^n,\: k=0,1,\ldots,m, $ } \par which represent the (unknown) solution at the equidistant points \par \qc{ $ \tau_{j}=\tau_{j,0}<\tau_{j,1}<\tau_{j,2}<\ldots<\tau_{j,m}=\tau_{j+1}, $ } \par given by \par \qc{ $ \tau_{j,k}=\tau_{j}+\frac{k}{m}(\tau_{j+1}-\tau_{j}),\: \: j=0,1,\ldots,N-1,\: k=0,1,\ldots,m. $ } \par The polynomial $u^{(j)}(\tau)$ are represented via the interpolation formula \par \qc{ $ u^{(j)}(\tau)=\Sum {i=0} {m} {u^{j,i}l_{j,i}(\tau)}, $ } \par where $l_{j,i}(\tau)$ are the {\bf Lagrange basis polynomialss}: \par \qc{ $ l_{j,i}(\tau)=\Prod {k=0,\:k\neq i} {m} {\frac{\tau-\tau_{j,k}}{\tau_{j,i}-\tau_{j,k}}}. $ } \par Equations (3) now can be considered as equations for $u^{j,i}$. The periodicity and the phase conditions in (2) are also substituted by their discrete approximations, \par \qc{ $ u^{0,0}=u^{N-1,m} $ }\qr{(4)} \par and \par \qc{ $ \Sum {j=0} {N-1} {\Sum {i=0} {m} {\omega_{j,i} \langle u^{j,i},\dot{v}^{j,i}\rangle}} = 0, $ }\qr{(5)} \par where $\dot{v}^{j,i}$ are the values of the derivative of the reference periodic solution at the points $\tau_{j,i}$, and $\omega_{j,i}$ are the {\bf Lagrange quadrature coefficients}. \par Equations (3),(4), and (5) compose a \jump kh_DefFunc {defining system} for the limit cycle curve. This system has $nmN+n+1$ scalar equations for the unknown components of $u^{j,i}$, the period $T_0$, and the parameter $\alpha$. The number of equations is one less than the number of unknowns. The Jacobian matrix of the system (3),(4),(5) has special structure, therefore a {\bf local elimination method} is used to solve linear equations appearing in the continuation [Doedel,1986]. \par The optimal choice of the collocation points $\{\zeta_{j,i}\}$ is to place them at the {\bf Gauss points}, which are the roots $\zeta_{j,i},\: i=1,2,\ldots,m,$ of the $m$th degree {\bf Legendre polynomial} relative to the subinterval $[\tau_j,\tau_{j+1}]$. In this case the solution of (3),(4),(5) will be an extremely accurate approximation to the solution of (2) at the {\bf mesh points}, namely \par \qc{ $ ||u(\tau_j)-u^{j,0}||=O(h^{2m}), $ } \par as $h=\max {0\leq j \leq N-1} {\Delta_j} \rightarrow 0$. Along the limit cycle curve, the mesh points are {\bf adapted} according to the requirement that the approximation error be uniformly distributed over the time interval. \par Several \jump test_lc {test functions} can be computed along the limit cycle curve to detect and process \jump test_lc {limit cycle singularities}. \par \par A limit cycle curve can be continued from a user-supplied \jump point_eq {point}, \jump lc_lc {cycle}, a \jump Hopf_eq {Hopf point}, as well as from a \jump branch_lc {branching point}, \jump LPC_lc {limit point}, \jump PD_lc {period-doubling point}, or a \jump NS_lc {Neimark-Sacker point} on the computed limit cycle curve. To \jump conti {continue} a limit cycle curve, one needs to set relevant {\bf parameters to start from} \par \par $\;\;\;$ \jump start_lcp {a point}\par $\;\;\;$ \jump start_lclc {a limit cycle}\par $\;\;\;$ \jump start_lcH {a Hopf bifurcation point}\par $\;\;\;$ \jump start_lcBP {a cycle branching point} $\;\;\;$ \jump start_lclc {a cycle limit point} $\;\;\;$ \jump start_lcPD {a period-doubling point}\par $\;\;\;$ \jump start_lclc {a Neimark-Sacker bifurcation point} \par \par via the \jump kh_StarterWin {starter window}. \end \% ode_lc \topic lc_lc {limit cycle} Given a system of ODEs, \par \qc{ $ \dot{x}=F(x),\: \: x \in {\bf R}^n, $ }\qr{(1)} \par a {\bf limit cycle} is an isolated closed orbit $L_0 \subset {\bf R}^n$ that is defined by a $T_0-$periodic solution to (1) \par \qc{ $ x_{0}(t+T_{0})=x_{0}(t). $ } \par \qc{\picture 0.4 fig1_13.fig} \par {\bf Multipliers} $\mu_{1},\mu_{2},\ldots,\mu_{n}$ of a limit cycle $L_0$ are the eigenvalues of the {\bf monodromy matrix} $M(T_{0})$. The matrix $M(t)$ satisfies the periodic {\bf variational equation}:\par \qc{ $ \dot{M}(t)=F_{x}(x_{0}(t))M(t) $ } \par with the initial condition $M(0)=I_{n}$, where $I_{n}$ is the $n\; \times \;n$ unit matrix. \par Matrix $M(T_{0})$ always has a trivial multiplier $\mu_{n}=1$. The {\bf Poincare map} $\xi \rightarrow P(\xi)$ is defined on a local crossection $\Sigma$ by orbits of (1) passing near $L_{0}$. The intersection point $\xi_0$ of the cycle with the hyper-plane $\Sigma$ is a fixed point of the Poincare map, $P(\xi_{0})=\xi_{0}$. The nontrivial multipliers coincide with the eigenvalues of the $(n-1)\:\times\:(n-1)$ matrix $P_{\xi}(\xi_{0})$. \par The limit cycle is {\bf hyperbolic} if there is no other multipliers with unit module. A hyperbolic limit cycle has invariant {\bf stable and unstable manifolds} of dimension equal to the number of multipliers with $|\mu|<1$ and $|\mu|>1$ plus one, respectively. If there is more than one multiplier with $|\mu|=1$, the cycle is called {\bf nonhyperbolic}. \par One can compute \par {\li 2 \desc 1 \jump ode_lc {a limit cycle curve}\> \par } \par starting from a known {\bf cycle} with a given period. In CONTENT, the time interval $[0,T_0]$ is automatically scaled to the unit interval $[0,1]$, while the period $T_0$ is treated as an extra active parameter. \end \% lc_lc \topic test_lc {test functions and singularities along the limit cycle curve} A special decomposition process (elimination of local variables) applied to the Jacobian matrix of the equations defining the \jump ode_lc {limit cycle curve} gives an approximation to the \jump lc_lc {monodromy matrix} $M=M(T_{0})$ in the form \par \qc{ $ M=-P^{-1}_{1}P_{0}, $ } \par where $P_{0,1}$ are $n\; \times \;n$ matrices. It also produces a $(2n+2)\: \times\: (2n+2)$ matrix $E$ of the reduced system. \par The following \jump kh_TestFunc {test-functions} can be computed along the limit cycle curve: \par \qc{ $ \table 4 rcl 0 0 {\psi_1}{\:=\:}{det(E),} {\psi_2}{\:=\:}{v_{K+1},} {\psi_3}{\:=\:}{det (M\:+\:I_{n})\:=\:\frac{1}{det (P_{1})} det (-P_{0}\:+\:P_{1}),} {\psi_4}{\:=\:}{det (M \otimes M\:-\:I_{J})\:=\:\frac{1}{det (P_{1}\otimes P_{1})} det (P_{0}\otimes P_{0}\: - \:P_{1}\otimes P_{1}) .} $ } \par Here $I_k$ is the unit $k\:\times\: k$ matrix; $v_{K+1}$ is the $\alpha-$component of the normalized tangent vector to the limit cycle curve $v=(v_1,v_2,\ldots,v_{K+1})\in {\bf R}^{K+1}$, $K=nmN+n+1, \: J=\frac{n(n-1)}{2}$; and $\otimes$ means the \jump bialt {bialternate product} of two matrices. \par The following {\bf singularities} can be detected and located as regular zeroes of the limit cycle test-functions: \par {\li 2 \desc 4 \jump branch_lc {Cycle branching point}:\> $\psi_1=0.$ \par \jump LPC_lc {Cycle limit point}:\> $\psi_2=0,\; \psi_1 \neq 0.$ \par \jump PD_lc {Period-doubling}:\> $\psi_3=0.$ \par \jump NS_lc {Neimark-Sacker bifurcation}:\> $\psi_4=0,\: \psi_{1} \psi_{2} \neq 0.$ \par } \end \%test_lc \topic branch_lc {cycle branching point} Generically, two branches $M_{1,2}$ of the \jump ode_lc {limit cycle curve} intersect transversally at a {\bf cycle branching point}: \par \qc{\picture 0.5 fig10_10.fig} \par CONTENT computes a vector that is orthogonal to vector $v$ (tangent to the primary branch $M_{1}$) and belongs to the plane spanned by $v$ and $V$ (vector tangent to the secondary branch $M_{2}$). This allows to switch branches. \end \% LPC_lc \topic LPC_lc {cycle limit point} At the critical parameter value, the system has a limit cycle $L_{0}$ with a nontrivial \jump lc_lc {multiplier} $\mu_1=1$. When the active parameter crosses the critical value, generically, two limit cycles, $L_1$ and $L_2$, collide and disappear. \par \qc{\picture 0.9 fig5_13.fig} \par The figure illustrates the bifurcation in the three-dimensional case. \par Generically, the Poincare map has a one-dimensional center manifold and its restriction to this manifold is locally topologically equivalent to the normal form \par \qc {$ \xi \:\rightarrow\: \beta\: + \:\xi\: \pm\: \xi^2 $} \par (see, for example \jump Ku95 {[Kuznetsov, 1995,1998]}). \end \% LPC_lc \topic PD_lc {period-doubling bifurcation} At the critical parameter value, the system has a limit cycle $L_{0}$ with a nontrivial \jump lc_lc {multiplier} $\mu_1=-1$. When the active parameter crosses the critical value, generically, a {\bf double-period cycle} $L_1$ bifurcates from the {\bf primary cycle} $L_0$ which changes its stability. \par \qc{\picture 0.9 fig5_14.fig} \par The figure illustrates the bifurcation in the three-dimensional case. \par Generically, the Poincare map has a one-dimensional center manifold and its restriction to this manifold is locally topologically equivalent to the normal form \par \qc {$ \xi\:\rightarrow\: -(1\:+\:\beta) \xi\: \pm \:\xi^3 $} \par (see, for example \jump Ku95 {[Kuznetsov, 1995,1998]}). \end \% PD_lc \topic NS_lc {Neimark-Sacker bifurcation} At the critical parameter value, the system has a limit cycle $L_{0}$ with a pair of complex \jump lc_lc {multipliers} on the unit circle \par \qc{ $ \mu_{1,2}=e^{\pm i\theta_{0}}. $ } \par When the active parameter crosses the critical value, generically, an {\bf invariant 2-torus} $T$ bifurcates from the cycle $L_0$ which changes its stability. \par \qc{\picture 0.9 fig5_15.fig} \par The figure illustrates the bifurcation in the three-dimensional case. See details, for example, in \jump Ku95 {[Kuznetsov, 1995,1998]}. \end \% NS_lc \topic start_lcp {starter parameters (limit cycle curve from a point)} This \jump kh_Starter {starter} produces necessary data to \jump kh_Generator {generate} the \jump ode_lc {limit cycle curve} from a \jump point_eq {point}. A point on the cycle and its period are assumed to be known sufficiently accurate. The starter integrates the \jump ode_or {orbit} from the point over the time interval equal to the period. Then, it produces the discretization data to continue the \jump ode_lc {limit cycle curve} using linear interpolation with current values of the discretization parameters. Only the interpolated orbit is visualized. \par It has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 17 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at a point on the cycle.\par {\li 1 $\alpha$} \> A list of parameter names and their values at a point on the cycle. The names of active parameters are highlighted.\par {\li 1 Period} \> The period of the initial cycle. Highlighted if active\par {\bf Discretization data}\> \par {\li 1 ntst} \> The number $N$ of mesh points.\par {\li 1 ncol} \> The number $m$ of collocation points.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lc {Singularity type}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Standard functions}\> \par {\li 1 \jump ode_stdfLC {Function name}} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables, parameters, and period. When it is specified and activated, it is called by the starter before computing the first point. Otherwise, CONTENT uses the values listed in the window. \par } \par \include ode_init \end \% start_lcp \topic start_lclc {starter parameters (limit cycle curve from a limit cycle)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_lc {limit cycle curve} from a \jump lc_lc {limit cycle} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 16 {\bf Initial point}\>\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\li 1 Period} \> The period of the initial cycle. Highlighted if active\par {\bf Discretization data}\> \par {\li 1 ntst} \> The number $N$ of mesh points.\par {\li 1 ncol} \> The number $m$ of collocation points.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lc {Singularity type}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Standard functions}\> \par {\li 1 \jump ode_stdfLC {Function name}} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitCycle} \> \jump kh_FuncPar {Functional parameter} assigning (approximate) initial cycle, corresponding parameters, and the period. When it is specified and activated, it is called by the starter before computing the first point on the curve. Otherwise, CONTENT uses the values listed in the window. \par } \par \include ode_initLC \end \% start_lclc \topic start_lcH {starter parameters (limit cycle curve from a Hopf point)} This \jump kh_Starter {starter} produces necessary data to \jump kh_Generator {generate} the \jump ode_lc {limit cycle curve} from a \jump Hopf_eq {Hopf bifurcation point}, using linear approximation to the bifurcating cycle. \par It has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 18 {\bf Initial point}\>\par {\li 1 $x$} \> A list of phase coordinate names and their values at a point on the cycle.\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Discretization data}\> \par {\li 1 ntst} \> The number $N$ of mesh points.\par {\li 1 ncol} \> The number $m$ of collocation points.\par {\bf Switch data}\> \par {\li 1 Amplitude} \> An approximate size of the bifurcating limit cycle.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lc {Singularity type}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Standard functions}\> \par {\li 1 \jump ode_stdfLC {Function name}} \> Compute toggle ({\bf yes}|{\bf no}). \par {\bf Set initial point}\> \par {\li 1 SetInitPoint} \> \jump kh_FuncPar {Functional parameter} assigning initial values to phase variables, parameters, and period. When it is specified and activated, it is called by the starter before computing the first point. Otherwise, CONTENT uses the values listed in the window. \par } \par \include ode_init \end \% start_lcH \topic start_lcPD {starter parameters (limit cycle curve from a period-doubling point)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_lc {limit cycle curve} from a \jump PD_lc {period-doubling point} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 17 {\bf Initial point}\>\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Discretization data}\> \par {\li 1 ntst} \> The number $N$ of mesh points. {\bf If a double-period cycle is thought to continue, multiply the old value of ntst by two}.\par {\li 1 ncol} \> The number $m$ of collocation points.\par {\bf Switch data}\> \par {\li 1 Shift} \> An approximate distance of the double-period cycle from the original limit cycle.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch to switch}\> \par {\li 1 Branch} \> Switch toggle ({\bf original}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lc {Singularity type}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Standard functions}\> \par {\li 1 \jump ode_stdfLC {Function name}} \> Compute toggle ({\bf yes}|{\bf no}). \par } \end \% start_lcPD \topic start_lcBP {starter parameters (limit cycle curve from a cycle branching point)} The \jump kh_Starter {starter} which produces necessary data to \jump kh_Generator {generate} the \jump ode_lc {limit cycle curve} from a \jump branch_lc {cycle branching point} has the following parameters that can be modified via the \jump kh_StarterWin {starter window.} \par \par { \desc 17 {\bf Initial point}\>\par {\li 1 $\alpha$} \> A list of parameter names and their values at the initial point. The names of active parameters are highlighted.\par {\bf Discretization data}\> \par {\li 1 ntst} \> The number $N$ of mesh points.\par {\li 1 ncol} \> The number $m$ of collocation points.\par {\bf Switch data}\> \par {\li 1 Shift} \> An approximate distance of the double-period cycle from the original limit cycle.\par {\bf Jacobian data}\> \par {\li 1 Increment} \> Increment to approximate partial derivatives by finite differences.\par {\bf Branch to switch}\> \par {\li 1 Branch} \> Switch toggle ({\bf primary}|{\bf secondary}).\par {\bf Monitor singularities}\>\par {\li 1 \jump test_lc {Singularity type}}\> Monitor toggle ({\bf yes}|{\bf no}).\par {\bf User defined functions}\>\par {\li 1 \jump kh_UserFunc {Function name}}\> Process toggle ({\bf ignore}|{\bf monitor}|{\bf detect}|{\bf append}).\par {\bf Standard functions}\> \par {\li 1 \jump ode_stdfLC {Function name}} \> Compute toggle ({\bf yes}|{\bf no}). \par } \end \% start_lcBP \topic ode_initLC {-} The function SetInitCycle is called by CONTENT several times. It has one input parameter {\bf ip} which specifies the purpose of the call: \par {\desc 4 {\bf ip}\>{\bf Meaning}\par -1\>initialize;\par $\:\:n$\>assign time and phase variables corresponding to n-th point on the cycle (n=0,1,...).\par -2\>terminate;\par} \par This function should return 0 if it be called more times. Return value 1 indicates the last point on the cycle. The names of the time, phase variables, and parameters should be used as specified in the \jump kh_RHSWin {system specification dialog box.} The {\bf Period} of the cycle is computed automatically as the difference of the time values at the first and the last points, which are assumed to be close in the phase space. \par In the following examples, T denotes time, X and Y are the names of the phase variables, and OMEGA is the parameter. \par The first example shows how to setup an initial cycle given by the explicit formulas:\par \qc{ $ X(T)=cos(\omega*T),\:\:Y(T)=sin(\omega*T). $ } \par on a uniform mesh. The corresponding function SetInitCycle can be specified as following: \par \verbatim 12 double PERIOD; if (ip<0) { PERIOD=8.0*atan(1.0)/OMEGA; return 0; } if (ip>=0 && ip