Computer-assisted proofs of Hopf bubbles and degenerate Hopf bifurcations

The Hopf bifurcation is generic in one-parameter vector fields of ODEs. Roughly speaking, this means it is prevalent in a topological sense among all vector fields that possess an equilibrium point with a pair of complex-conjugate imaginary eigenvalues. In the world of two-parameter vector fields, the richness of possible bifurcations is much greater. In this paper, we are primarily interested in a type of degenerate Hopf bifurcation that results from a particular failure of the eigenvalue transversality condition. The failure we are interested in is where the curve of eigenvalues has a quadratic tangency with the imaginary axis, so the eigenvalues do not cross transversally and, importantly, do not cross the imaginary axis at all.

Before surveying some literature about this bifurcation pattern, let us construct a minimal working example. Consider the ODE system

	$\displaystyle\dot{x}$	$\displaystyle=\beta x-y-x(x^{2}+y^{2}+\alpha^{2})$		(1)
	$\displaystyle\dot{y}$	$\displaystyle=x+\beta y-y(x^{2}+y^{2}+\alpha^{2}).$		(2)

The reader familiar with the normal form of the Hopf bifurcation should find this familiar, but might be be unnerved by the $\alpha^{2}$ term, which is not present in the typical normal form. The linearization at the equilibrium $(0,0)$, produces the matrix

$$\left(\begin{array}[]{cc}\beta-\alpha^{2}&-1\\ 1&\beta-\alpha^{2}\end{array}\right),$$

which has the pair of complex-conjugate eigenvalues $\lambda=\beta-\alpha^{2}\pm i$. Treating $\alpha$ as being a fixed constant, we have a supercritical Hopf bifurcation as $\beta$ passes through $\alpha^{2}$. Conversely, if $\beta>0$ is fixed and we interpret $\alpha$ as the parameter, there are two supercritical Hopf bifurcations: as $\alpha$ passes through $\pm\sqrt{\beta}$. However, something problematic happens if $\beta=0$: the eigenvalue branches $\alpha\mapsto-\alpha^{2}\pm i$ are tangent to the imaginary axis at $\alpha=0$ and do not cross at all, so the non-degeneracy condition of the Hopf bifurcation fails.

More information can be gleaned by transforming to polar coordinates. Setting $x=r\cos\theta$ and $y=r\sin\theta$ results in the equation for the radial component

$$\dot{r}=r(\beta-r^{2}-\alpha^{2}),$$

while $\dot{\theta}=1$. There is a nontrivial periodic orbit (in fact, limit cycle) if $\beta-\alpha^{2}>0$, with amplitude $\sqrt{\beta-\alpha^{2}}$. Alternatively, there is a surface of periodic orbits described by the graph of $\beta=r^{2}+\alpha^{2}$. For fixed $\beta>0$, the amplitude $r$ of the periodic orbit, as a function of $\alpha$, is $\alpha\mapsto\sqrt{\beta-\alpha^{2}}$, whose graph is half of an ellipse with semi-major axis $2\sqrt{\beta}$ and semi-minor axis $\sqrt{\beta}$.

This non-local bifurcation structure, characterized by the connection of two Hopf bifurcations by a one-parameter branch of periodic orbits, has been given several names in different scientific fields. In intracellular calcium, it is frequently called a Hopf bubble [7, 16, 22, 27]. In infectious-disease modelling, where the bifurcation typically occurs at an endemic equilibrium, the accepted term is endemic bubble [8, 17, 18, 23, 26]. A more general definition of a (parametrically) non-local structure called bubbling is given in [15]. In the present work, we will refer to the structure as a Hopf bubble, since that name is descriptive of the geometric picture, the bifurcation involved, and is sufficiently general to apply in different scenarios in a model-independent way.

Hopf bubbles can, in many instances, be understood as being generated by a codimension-two bifurcation; see Section 1.3. However, this is not to say that they are rare. Aside from the applications in calcium dynamics and infectious-disease modelling in the previous paragraph, they have been observed numerically in models of neurons [1], condensed-phase combustion [20], predator-prey models [3], enzyme-catalyzed reactions [12], and a plant-water ecosystem model [31]. A recent computer-assisted proof also established the existence of a Hopf bubble in the Lorenz-84 model [30]. We will later refer to the degenerate Hopf bifurcation that gives rise to Hopf bubbles as a bubble bifurcation. The literature seems to not have an accepted name for this bifurcation, so we have elected to give it one here.

Bubble bifurcations are one type of dgenerate Hopf bifurcation. Another is the Bautin bifurcation, which occurs at a Hopf bifurcation whose first Lyapunov coefficient vanishes [14]. Like the bubble bifurcation, the Bautin bifurcation is a codimension-two bifurcation of periodic orbits. To illustrate this bifurcation, consider the planar normal form

	$\displaystyle\dot{x}$	$\displaystyle=\beta x-y-x(x^{2}+y^{2})(\alpha-x^{2}-y^{2})$		(3)
	$\displaystyle\dot{y}$	$\displaystyle=x+\beta y-y(x^{2}+y^{2})(\alpha-x^{2}-y^{2}).$		(4)

Transforming to polar coordinates, the angular component decouples producing $\dot{\theta}=1$, while the radial component gives

$$\dot{r}=r(\beta+\alpha r^{2}-r^{4}).$$

Nontrivial periodic orbits are therefore determined by the zero set of $r\mapsto\beta+\alpha r^{2}-r^{4}$, which defines a two-dimensional smooth manifold. See later Figure 15 for a triangulation of (part of) this manifold.

The bubble bifurcation has been fully characterized by LeBlanc [17], using the center manifold reduction and normal form theory for functional differential equations [9] and a prior classification of degenerate Hopf bifurcations for ODEs by Golubitsky and Langford [11]. The study of the Bautin bifurcation goes back to 1949 with the work of Bautin [2], and a more modern derivation based on normal form theory appears in the textbook of Kuznetsov [14].

Normal form theory and centre manifold reduction are incredibly powerful, providing both the direction of the bifurcation and criticality of bifurcating periodic orbits. The drawback is that they are computationally intricate and inherently local. In the present work, we advocate for an analysis of degenerate Hopf bifurcations in ODE and delay differential equations (DDE) by way of two-parameter continuation, desingularization and computer-assited proofs. We will discuss the latter topics in the next section.

The intuition behind our continuation idea can be pictorially seen in Figure 1, and understood analytically by way of our minimal working example, (1)–(2), and its surface of periodic orbits described by the equation $\beta=\alpha^{2}+r^{2}$. There is an implicit relationship between two scalar parameters ($\alpha$ and $\beta$) and the set of periodic orbits of the dynamical system. At a Hopf bifurcation, one of these periodic orbits should retract onto a steady state. In other words, their amplitudes (relative to steady state) should become zero. Using a desingularization technique to topologically separate periodic orbits from steady states that they bifurcate from, we can use two-parameter continuation to numerically compute and continue periodic orbits as they pass through a degenerate Hopf bifurcation. With computer-assisted proof techniques, we can prove that the numerically-computed objects are close to true periodic orbits, with explicit error bounds. Further a posteriori analysis can then be used to prove the existence of a degenerate Hopf bifurcation based on the output of the computer-assisted proof of the continuation.

Computer-assisted proofs of Hopf bifurcations have been completed in [30] using a desingularization (sometimes called blow-up) approach, in conjunction with an a posteriori Newton-Kantorvich-like theorem. We lean heavily on these ideas in the present work. The former desingularization idea, which we adapt to delay equations in Section 3.1, is used to resolve the fact that, at the level of a Fourier series, a periodic orbit that limits (as a parameter varies) to a fixed point at a Hopf bifurcation is, itself, indistinguishable from that fixed point. Our methods of computer-assisted proof are based on contraction mappings, and it is critical that the objects we prove are isolated as fixed points. The desingularization idea exploits the fact that an amplitude-based change of variables can be used to develop an equivalent problem where the representative of a periodic orbit consists of a Fourier series, a fixed point and a real number. The latter real number represents a signed amplitude. This reformulation results in fixed points being spatially islolated from periodic orbits, thereby allowing contraction-based computer-assisted proofs to succeed.

The other big point of inspiration in this work is validated multi-parameter continuation. The technique was developed in [10] for continuation in general Banach spaces, and applied to some steady state problems for PDEs. We will overview this method in Section 2. It is based on a combination of a two-dimensional analogue of pseudo-arclength continuation and a posteriori, uniform validation of zeroes of nonlinear maps by a Newton-Kantorovich-like theorem applied to a Newton-like operator.

The main contribution of the present work is a general-purpose code for computing and proving two-parameter families of periodic orbits in polynomial delay differential equations. Equations of advanced or mixed-type can similarly be handled; there is no restriction whether delays are positive or negative. Ordinary differential equations can also be handled as a special case. Orbits can be proven in the vicinity of (and at) Hopf bifurcations, whether these are non-degenerate or degenerate. The first major release of the library BiValVe (Bifurcation Validation Venture, [4]) is being made in conjunction with the present work, and builds on an earlier version of the code associated to the work [29, 30]. Some non-polynomial delay differential equations can be handled using the polynomial embedding technique. The existence of Hopf bifurcation curves and degenerate Hopf bifurcations can then be completed by post-processing of the output of the computer-assisted proof.

To explore the applicability of our validated numerical methods, we explore Hopf bifurcations and degenerate Hopf bifurcations in

• •

  the extended Lorenz-84 model (ODE)

• •

  a time-delay SI model (DDE)

• •

  the FitzHugh-Nagumo equation (ODE)

• •

  an ODE with complicated branches of periodic orbits (ODE)

The first two examples replicate and extend some of the analysis appearing in [17, 30] using our computational scheme. The third examples provides, to our knowledge, the first analytically verified results on degenerate Hopf bifurcations and Hopf bubbles in that equation. The final example is a carefully designed ODE that exhibits the degenerate Hopf bifurcation in addition to folding and pinching in some projections of the manifold.

Section 2 serves as an overview of two-parameter continuation, both in the finite-dimensional and infinite-dimensional case. We introduce the continuation scheme for periodic orbits near Hopf bifurcations in Section 3 in the general case of delay differential equations. Some technical bounds for the computer-assisted proofs in Section 4. A specification to ordinary differential equations is presented in Section 5. In Section 6, we connect the computer-assisted proofs of the manifold of periodic orbits to analytical proofs of Hopf bifurcation curves, Hopf bubbles, and degenerate Hopf bifurcations. Our examples are presented in Section 7, and we complete a discussion and comment on future research directions in Section 8.

Given $n\in\mathbb{N}$, denote $\mathbb{C}^{\mathbb{Z}}_{n}$ the vector space of $\mathbb{Z}$-indexed sequences of elements of $\mathbb{C}^{n}.$ Denote $\mbox{Symm}(\mathbb{C}^{\mathbb{Z}}_{n})$ the proper subspace of $\mathbb{C}^{\mathbb{Z}}_{n}$ consisting of symmetric sequences; $z\in\mbox{Symm}(\mathbb{C}^{\mathbb{Z}}_{n})$ if and only if $z\in\mathbb{C}^{\mathbb{Z}}_{n}$ and $z_{k}=\overline{z_{-k}}$ for all $k\in\mathbb{Z}$. For any subspace $U\subset\mathbb{C}^{\mathbb{Z}}_{n}$ closed under (componentwise) complex conjugation, denote $\operatorname{Symm}(U)=U\cap\operatorname{Symm}(\mathbb{C}^{\mathbb{Z}}_{n})$. We will sometimes drop the subscript $n$ on $\mathbb{C}^{\mathbb{Z}}_{n}$ when the context is clear.

Given $\nu>1$, we denote $\ell_{\nu}^{1}(\mathbb{C}^{n})$ the subspace of $\mathbb{C}_{n}^{\mathbb{Z}}$ whose elements $z$ satisfy $||z||_{\nu}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\sum_{k\in\mathbb{Z}}|z_{k}|\nu^{|k|}<\infty$. $|\cdot|$ will always be taken to be the norm on $\mathbb{C}^{n}$ induced by the standard inner product. Denote $K_{\nu}(\mathbb{C}^{n})$ the subspace of $\mathbb{C}_{n}^{\mathbb{Z}}$ whose elements $z$ satisfy $||z||_{\nu,K}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\sum_{k\in\mathbb{Z}}(\nu^{|k|}/|k|)|z_{k}|<\infty.$ Introduce a bilinear form on $\langle\cdot,\cdot\rangle$ on $\ell_{\nu}^{1}(\mathbb{C})$ as follows:

$$\langle v,w\rangle=\sum_{k\in\mathbb{Z}}v_{k}\overline{w_{k}}.$$

For $v,w\in\mathbb{C}^{\mathbb{Z}}_{1}$, define their convolution $v*w$ by

$\displaystyle(v*w)_{k}=\sum_{k_{1}+k_{2}=k}v_{k_{1}}w_{k_{2}}$

(5)

whenever this series converges. Convolution is commutative and associative for sequences in $\ell_{\nu}^{1}(\mathbb{C})$. In this case, we define multiple convolutions (e.g. triple convolutions $a*b*c$) inductively, by associativity.

A function $g:\ell_{\nu}^{1}(\mathbb{C}^{n})\rightarrow\ell_{\nu}^{1}(\mathbb{C})$ is a convlution polynomial of degree $q$ if

$\displaystyle g(z)$

$\displaystyle=c_{0}+\sum_{k=1}^{q}\sum_{p\in\mathcal{M}_{k}}c_{p}(z_{p_{1}}*\cdots*z_{p_{k}}),$

for $c_{p}\in\mathbb{C}$, where $\mathcal{M}_{k}$ denotes the set of $k$-element multisets of $\{1,\dots,n\}$, and each multiset $p$ is identified with the unique tuple $(p_{1},\dots,p_{k})\{1,\dots,n\}^{k}$ such that $p_{i}\leq p_{i+1}$ for all $i=1,\dots,k-1$. Analogously, $g:\ell_{\nu}^{1}(\mathbb{C}^{n})\rightarrow\ell_{\nu}^{1}(\mathbb{C}^{m})$ is a convolution polynomial of degree $q$ if $z\mapsto g(z)_{\cdot,j}\in\ell_{\nu}^{1}(\mathbb{C})$ is a convolution polynomial of degree $q$, for $j=1\dots,m$.

If $X$ is a vector space, $0_{X}$ will denote the zero vector in that space. If $X$ is a metric space and $U\subset X$, we denote $U^{\circ}$ its interior, $\partial U$ its boundary, and $\overline{U}$ its closure.

An interval vector in $\mathbb{R}^{k}$ for some $k\geq 1$ is a subset of the form $v=[a_{1},b_{1}]\times\cdots\times[a_{k},b_{k}]$ for real scalars $a_{j},b_{j}\in\mathbb{R}$, $j=1,\dots,k$. If $\mathbb{R}^{k}$ is equipped with a norm $||\cdot||$ we define $||v||=\sup_{w\in v}||w||$. Similarly, an interval vector in $\mathbb{C}^{k}$ is a product $v=A_{1}\times\cdots\times A_{k}$, where each $A_{j}$ is a closed disc in $\mathbb{C}$. We define the norm of a complex interval vector as $||v||=\sup_{w\in v}||w||$ whenever $||\cdot||$ is a norm on $\mathbb{C}^{k}$. In each case, be it real or comlex, the unit interval vector is the unique interval vector $V$ with the property that $||V||=1$, and for any other interval vector $v$ with $||v||=1$, the inclusion $v\subseteq V$ is satisfied.

Let $\mathcal{X}$ and $\mathcal{Y}$ be finite-dimensional vector spaces over the field $\mathbb{R}$, with $\dim\mathcal{X}=\dim\mathcal{Y}+2$, and consider a map $\mathcal{G}:\mathcal{X}\rightarrow\mathcal{Y}$. We are interested in the zero set of $\mathcal{G}$. Given the codimension of $\mathcal{G}$, we expect zeroes to be in a two-dimensional manifold in $\mathcal{X}$. In what follows, equality with zero should be understood in the sense of numerical precision. For example, if we write $\mathcal{G}(x)=0$, what we mean is that $||\mathcal{G}(x)||$ is machine precision small.

Let $\hat{x}_{0}\in\mathcal{X}$ satisfy $G(\hat{x}_{0})=0$, and suppose $D\mathcal{G}(\hat{x}_{0})$ has two-dimensional kernel. This property is generic. Let $\{\hat{\Phi}_{1},\hat{\Phi}_{2}\}$ span the kernel. The tangent plane $\mathcal{T}_{\hat{x}_{0}}\mathcal{M}$ at $\hat{x}_{0}$ of the two-dimensional solution manifold $\mathcal{M}$ is therefore spanned by $\{\hat{\Phi}_{1},\hat{\Phi}_{2}\}$. If $\epsilon_{1},\epsilon_{2}$ are small, we have

$$\mathcal{G}(\hat{x}_{0}+\epsilon_{1}\hat{\Phi}_{1}+\epsilon_{2}\hat{\Phi}_{2})=\mathcal{G}(\hat{x}_{0})+\epsilon_{1}D\mathcal{G}(\hat{x}_{0})\hat{\Phi}_{1}+\epsilon_{2}D\mathcal{G}(\hat{x}_{0})\hat{\Phi}_{2}+o(|\epsilon|)=o(|\epsilon|)$$

by Taylor’s theorem, provided $\mathcal{G}$ is differentiable at $\hat{x}_{0}$, and $|\epsilon|$ is any Euclidean norm of $\epsilon=(\epsilon_{1},\epsilon_{2})$. If $\mathcal{G}$ is twice continuously differentiable, the error is improved to $O(|\epsilon|^{2})$. Therefore, new candidate zeroes of $\mathcal{G}$ can be computed using $\hat{x}_{0}$ and a basis for the tangent space. This idea is at the heart of the continuation.

The continuation from $\hat{x}_{0}$ is done by way of iterative triangulation of the manifold $\mathcal{M}$. First, we compute an orthonormal basis of $\mathcal{T}_{\hat{x}_{0}}\mathcal{M}$ by applying the Gram-Schmidt process to $\{\hat{\Phi}_{1},\hat{\Phi}_{2}\}$; see Section 2.2 for some technical details. Using this orthonormal basis to define a local coordinate system, six vertices of a regular hexagon are computed around $\hat{x}_{0}$ at a specified distance $\sigma$ from $\hat{x}_{0}$; see Figure 2. Let these vertices be denoted $\hat{x}_{1}$ through $\hat{x}_{6}$, arranged in counterclockwise order (relative to the local coordinate system) around $\hat{x}_{0}$. Note that this means

$$\hat{x}_{i}=\hat{x}_{0}+\epsilon_{j,1}\hat{\Phi}_{1}+\epsilon_{j,2}\hat{\Phi}_{2}$$

for some small $\epsilon_{j,1},\epsilon_{j,2}$, which means that $G(\hat{x}_{j})\approx 0$ for $j=1,\dots,6$. Each of these candidate zeroes $\hat{x}_{j}$ are then refined by applying Newton’s method to the map

$\displaystyle x\mapsto\mathcal{G}_{j}(x)=\left(\begin{array}[]{c}\mathcal{G}(x)\\ \langle\hat{\Phi}_{1},x-\hat{x}_{j}\rangle\\ \langle\hat{\Phi}_{2},x-\hat{x}_{j}\rangle\end{array}\right).$

(9)

The two added inner product equations ensure isolation of the solution (hence quadratic convergence of the Newton iterates) and that the Newton correction is perpendicular to the tangent plane. This map can similarly be used to refine the original zero $\hat{x}_{0}$.

This initial hexagonal “patch” is itself formed by six triangles; see Figure 2. The continuation algorithm proceeds by selecting one of the boundary vertices (i.e. one of the vertices $\hat{x}_{1}$ through $\hat{x}_{6}$) and attempting to “grow” the manifold further. We describe this “growth” phase below. However, first, some terminology. The vertices will now be referred to as nodes. An edge is a line segment connecting two nodes, and they will be denoted by pairs of nodes: $\{\hat{x}_{i},\hat{x}_{j}\}$ for node $\hat{x}_{i}$ connected to $\hat{x}_{j}$. Two nodes are incident if they are connected by an edge. A simplex is the convex hull of three edges that form a triangle. Once the hexagonal patch is created, the data consists of:

• •

  the nodes $\hat{x}_{0},\dots,\hat{x}_{6}$;

• •

  the “boundary edges” $\{\hat{x}_{1},\hat{x}_{2}\},\dots,\{\hat{x}_{5},\hat{x}_{6}\},\{\hat{x}_{6},\hat{x}_{1}\}$;

• •

  the six simplices formed by triangles with $\hat{x}_{0}$ as one of the nodes.

Two simplices are adjacent if they share an edge. An internal simplex is a simplex that is adjacent to three other simplices, and it is a boundary simplex otherwise. Therefore, the six simplices of the initial patch are all considered boundary since they are adjacent to exactly two others. Similarly, an edge of a simplex can be declared boundary or internal; internal edges are those that are shared with another simplex, and boundary edges are not. A boundary node is any node on a boundary edge, a frontal node is a boundary node on an edge shared by two simplices, and an internal node is a node that is not a boundary node.

Let $x$ be a chosen frontal node (see Section 2.2.3 for a general discussion on frontal node selection). By construction, $x$ is part of (at least) three distinct edges, two of which connect to boundary nodes, and at least one that connects to an internal node. The algorithm to grow a simplex is as follows. See Figure 3 for a visualization.

1. 1.

   Compute an orthonormal basis for the tangent space $\mathcal{T}_{x}\mathcal{M}$.

2. 2.

   Compute the (average) gap complement direction. Let $x^{\circ}_{1},\dots,x^{\circ}_{m}$ denote the list of internal nodes nodes incident to $x$. The gap complement direction is $y_{c}=-x+\frac{1}{m}\sum_{i=1}^{m}x_{i}^{\circ}$.

3. 3.

   Let $x_{1}$ and $x_{2}$ denote the boundary nodes incident to $x$; form the edge directions $y_{1}=x_{1}-x$ and $y_{2}=x_{2}-x$.

4. 4.

   Orthogonally project $y_{c}$, $y_{1}$ and $y_{2}$ onto the tangent plane $\mathcal{T}_{x}\mathcal{M}$. Let these projections be denoted $Py_{c}$, $Py_{1}$ and $Py_{2}$. Introduce a two-dimensional coordinate system on $\mathcal{T}_{x}\mathcal{M}$ by way of a “unitary” (see Section 2.2) transformation to $\mathbb{R}^{2}$. Let $\tilde{P}y_{c}$, $\tilde{P}y_{1}$ and $\tilde{P}y_{2}$ denote the representatives in $\mathbb{R}^{2}$.

5. 5.

   In the local two-dimensional coordinate system, compute the counter-clockwise (positive) angle $\theta_{1}$ required to complete a rotation from $\tilde{P}y_{c}$ to $\tilde{P}y_{1}$, and the counter-clockwise (positive) angle $\theta_{2}$ required for rotation from $\tilde{P}y_{c}$ to $\tilde{P}y_{2}$. The gap angle $\gamma$ and orientation $\rho$ is

   $$\gamma=\max\{\theta_{1}-\theta_{2},\theta_{2}-\theta_{1}\},\hskip 28.45274pt\rho=\left\{\begin{array}[]{ll}1,&\theta_{2}>\theta_{1}\\ 2&\theta_{2}<\theta_{1}\end{array}\right.$$

6. 6.

   If $\gamma<\frac{\pi}{6}$ then close the gap: add the the triangle formed by the nodes $\{x,x_{1},x_{2}\}$ to the list of simplices, flag it as a boundary simplex, flag $x$ as an internal node, and conclude the growth step. Otherwise, proceed to step 7.

7. 7.

   Generate the predictor “fan” in $\mathbb{R}^{2}$: let $k=\min\left\{1,\lfloor 3\theta/\pi\rfloor\right\}$ and define the predictors

   $$\tilde{P}\tilde{y}_{j}=R(j\gamma/k)\tilde{P}y_{\rho},\hskip 28.45274ptj=1,\dots,k,$$

   for $R(\theta)$ the $2\times 2$ counterclockwise rotation matrix through angle $\theta$.

8. 8.

   Invert the unitary transformation and map $\tilde{P}\tilde{y}_{j}$ into the tangent plane $\mathcal{T}_{x}\mathcal{M}$; let the result be the vectors $P\tilde{y}_{j}$, $j=1,\dots,k$.

9. 9.

   Define predictors $\hat{x}_{j}=x+\sigma P\tilde{y}_{j}$ for $j=1,\dots,k$, where $\sigma$ is a user-specified step size. Refine them using Newton’s method applied to (9), where $\hat{\Phi}_{1}$ and $\hat{\Phi}_{2}$ are now the orthonormal basis for $\mathcal{T}_{x}\mathcal{M}$.

10. 10.

    Add the triangles formed by nodes $\{x,\hat{x}_{1},\hat{x}_{2}\}$, $\{x,\hat{x}_{2},\hat{x}_{3}\}$, $\dots$, $\{x,\hat{x}_{k-1},\hat{x}_{k}\}$. To the list of simplices, flag them as boundary simplices, and flag $x$ as an internal node. Do the same with the triangles formed by $\{x,\hat{x}_{1},*\}$ and $\{x,\hat{x}_{k},*\}$, where $*$ denotes ones of $x_{1}$ and $x_{2}$, depending on orientation $\rho$.

At the end of the simplex growth algorithm, one or more boundary simplices is added to the dictionary, and one additional node will be flagged as internal. The structure of this algorithm ensures that each simplex added to the dictionary will be adjacent to exactly two others, indicating that our convention of internal vs. boundary simplices is effective. Once the growth algorithm is complete, another frontal node is selected and the growth phase is repeated. This continues until a sufficient portion of the manifold has been computed (i.e. a user-specified exit condition is reached), or until a Newton’s method correction fails.

Here we collect some remarks concerning implementation of the finite-dimensional continuation. These comments may be applicable for general continuation problems, but we will often emphasize our specific situation which is continuation of periodic orbits in delay differential equations.

Two-parameter continuation can be introduced generally in a Banach space, but here we take the perspective that $\mathcal{G}:\mathcal{X}\rightarrow\mathcal{Y}$ of Section 2.1 is a projection of a nonlinear map $G:X\rightarrow Y$, and we are actually interested in computing the zeroes of $G$, rather than those of $\mathcal{G}$. That is to say, there exist projection operators $\pi_{\mathcal{X}}:X\rightarrow\mathcal{X}$, $\pi_{\mathcal{Y}}:Y\rightarrow\mathcal{Y}$, and associated embeddings $i_{\mathcal{X}}:\mathcal{X}\hookrightarrow X$, $i_{\mathcal{Y}}:\mathcal{Y}\hookrightarrow Y$, such that $\mathcal{G}=\pi_{\mathcal{Y}}G\circ i_{\mathcal{X}}$.

While the zeroes of $G$ are what we want, we will still use the approximate zeroes of $\mathcal{G}$ to rigorously enclose them. Abstractly, let $\hat{x}_{0},\hat{x}_{1},\hat{x}_{2}\in\mathcal{X}$ be the three nodes of a simplex. Then we may introduce a coordinate system on this simplex as follows: write each element $\hat{x}_{s}$ of this simplex as a unique linear combination

$\displaystyle\hat{x}_{s}=\hat{x}_{0}+s_{1}(\hat{x}_{1}-\hat{x}_{0})+s_{2}(\hat{x}_{2}-\hat{x}_{0})$

(10)

for $s=(s_{1},s_{2})\in\Delta=\{(a,b)\in\mathbb{R}^{2}:0\leq a,b\leq 1,\hskip 2.84526pt0\leq a+b\leq 1\}$. Let $\{\hat{\Phi}_{j,1},\hat{\Phi}_{j,2}\}$ for $j=0,1,2$ denote an orthonormal basis for the kernel of $D\mathcal{G}(\hat{x}_{j})$. We can then form the interpolated kernels

$$\hat{\Phi}_{s,i}=\hat{\Phi}_{0,i}+s_{1}(\hat{\Phi}_{1,i}-\hat{\Phi}_{0,i})+s_{2}(\hat{\Phi}_{2,i}-\hat{\Phi}_{0,i})$$

for $i=1,2$. We introduce a nonlinear map $G_{s}:X\rightarrow Y\times\mathbb{R}^{2}$,

$\displaystyle G_{s}(x)=\left(\begin{array}[]{c}G(x)\\ \langle\hat{\Phi}_{s,1},\pi_{\mathcal{X}}x-\hat{x}_{s}\rangle\\ \langle\hat{\Phi}_{s,2},\pi_{\mathcal{X}}x-\hat{x}_{s}\rangle\end{array}\right).$

(14)

The objective is to prove that $G_{s}$ has a unique zero close to $\hat{x}_{s}$ for each $s\in\Delta$. If this can be proven, and in fact $G_{s}$ is $C^{k}$ for some $k\geq 1$, then one can prove [10] that the the zero set of $G$ is a $C^{k}$, two-dimensional manifold. If this same can be proven for a collection of simplices, then the zero set is globally (i.e. on the union of the cobordant simplicial patches) a $C^{k}$ manifold over all patches that can be proven; that is, the transition maps between patches are $C^{k}$.

The radii polynomial approach is essentially a Newton-Kantorovich theorem with an approximate derivative, approximate inverse, and domain parametrization. It will be used to connect the approximate zeroes of $\mathcal{G}$ with exact zeroes of $G$ from the previous section. We state it generally for a family of $s\in\Delta$ -dependent maps $F_{s}:X_{1}\rightarrow X_{2}$. We include a short proof for completeness.

For our problem, injectivity of $A_{s}$ will always follow from a successful verification of $p(r_{0})<0$. See Lemma 4.

Let $\hat{x}_{s}$ be the convex combination defined by (10) for the simplex nodes $\hat{x}_{0}$, $\hat{x}_{1}$ and $\hat{x}_{2}$. We will say this simplex has been validated if we successfully find a radius $r_{0}$ such that the conditions of the radii polynomial theorem are successful for the nonlinear map $G_{s}:X\rightarrow Y\times\mathbb{R}^{2}$.

In our implementation, we generate simplices “offline” first. This allows the workload to be distributed across several computers, since the validation step can be restricted to only a subset of the computed simplices. We do not implement a typical refinement procedure, where simplices that fail to validated are split, with more nodes added and corrected with Newton. Rather, we implement an adaptive refinement step, which can help with the validation if failure is primarily a result of interval over-estimation. See Section 4.2.1.

Theorem 1 guarantees that the map from the standard simplex to the zero set of $G_{s}(\cdot)$ is $C^{1}$. There is then a natural question as to the smoothness of the manifold obtained by gluing together the images of the $C^{1}$ maps. This is answered in the affirmative in [10], and is primarily a consequence of the numerical data being equal on cobordant simplices.

We begin by doing a “blowup” around the periodic orbit. Write $y=x+az$, for $x$ a candidate equilibrium point. Then we get the rescaled vector field

$\displaystyle\tilde{f}(z_{1},\dots,z_{J},x,a,\alpha,\beta)=\left\{\begin{array}[]{ll}a^{-1}(f(x+az_{1},\dots,x+az_{J},\alpha,\beta)-f(x,\dots,x,\alpha,\beta)),&a\neq 0\\ \sum_{j=0}^{J}d_{y_{j}}f(x,\dots,x,\alpha,\beta)z_{j},&a=0.\end{array}\right.$

The interpretation is that $||z||=1$, so that $a$ behaves like the relative norm-amplitude of the periodic orbit. The vector field above is $C^{k-1}$ provided the original function $f$ is $C^{k}$ and $x$ is an equilibrium point; that is, $f(x,\dots,x,\alpha,\beta)=0$. At this stage we can summarize by saying that our goal is to find a pair $(x,z)$ such

	$\displaystyle f(x,\dots,x,\alpha,\beta)$	$\displaystyle=0$
	$\displaystyle\dot{z}(t)$	$\displaystyle=\tilde{f}(z(t+\mu_{1}),\dots,z(t+\mu_{J}),x,a,\alpha,\beta),$
	$\displaystyle||z||$	$\displaystyle=1$

where $z$ is $\omega$-periodic for an unknown period $\omega$; equivalently, the frequency of $z$ is $\psi=\frac{2\pi}{\omega}$.

In what follows, it will be beneficial for the vector field $\tilde{f}$ to be polynomial. This is because we will make use of a Fourier spectral method, and polynomial nonlinearities in the function space translate directly to convolution-type nonlinearities on the sequence space of Fourier coefficients. While we can make due with non-polynomial nonlinearities, it is greatly simplifies the computer-assisted proof if they are polynomial. To fix this, we generally advocate the use of the polynomial embedding technique. The idea is that many analytic, non-polynomial functions are themselves solutions of polynomial ordinary differential equations. The reader may consult [van den Berg, Groothedde, Lessard] for a brief survey of this idea in the context of delay differential equations. See Section 7.2 for a specific example.

Applying the polynomial embedding procedure always introduces additional scalar differential equations. If we need to introduce $m$ extra scalar equations to get a polynomial vector field, this will also introduce $m$ natural boundary conditions that fix the initial conditions of the new components. As a consequence, we need to bring in $m$ unfolding parameters to balance the system. This is accomplished in a problem-specific way; see [van den Berg, Groothedde, Lessard] for some general guidelines and a discussion for the need of these extra unfolding parameters. By an abuse of notation, we assume $\tilde{f}$ is polynomial (that is, the embedding has already been performed), and we write it as

$$\tilde{f}(z_{1},\dots,z_{J},x,a,\alpha,\beta,\eta),$$

where $\eta\in\mathbb{R}^{m}$ is a vector representing the unfolding parameter, and we now interpret $\tilde{f}:(\mathbb{R}^{n+m})^{J+1}\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{n+m}$. We then write the natural boundary condition corresponding to the polynomial embedding as

$$\theta_{BC}(z(0),x,a,\alpha,\beta,\eta)=0\in\mathbb{R}^{m}.$$

It can also be useful to eliminate non-polynomial parameter dependence from the vector field, especially if the latter has high-order polynomial terms with respect to the state variable $z$. This can often be accomplished by introducing extra scalar variables. For example, if $\tilde{f}$ is

$$\alpha e^{-px}z_{1}-\beta z_{2}^{2}+\eta_{1}$$

and we want to eliminate the non-polynomial term $e^{-px}$ from the vector field, then we can introduce a new variable $\eta_{2}$ and impose the equality constraint $0=\eta_{2}-e^{-px}$. The result is that the vector field becomes

$$\alpha\eta_{2}z_{1}-\beta z_{2}^{2}+\eta_{1}.$$

Since this operation introduces new variables and additional constraints, we incorporate the constraints as extra components in the natural boundary condition function $\theta_{BC}$ of the polynomial embedding. Since this type of operation will introduce an equal number of additional scalar variables and boundary conditions, we will neglect them from the dimension counting.

The final thing we need to take into account is that every periodic orbit is equivalent to a one-dimensional continuum by way of phase shifts. Since our computer-assisted approach to proving periodic orbits is based on Newton’s method and contraction maps, we need to handle this lack of isolation. This can be done by including a phase condition. In this paper we will make use of an anchor condition; we select a periodic function $\hat{z}$ having the same period as $z$, and we require that

$$\int\langle z(s),\hat{z}^{\prime}(s)\rangle ds=0.$$

We are ready to write down a zero-finding problem for our rescaled periodic orbits. First, combining the work of the previous sections, we must simultaneously solve the equations

$$\begin{cases}\dot{z}=\tilde{f}(z(t+\mu_{1}),\dots,z(t+\mu_{J}),x,a,\alpha,\beta,\eta),&\quad\text{(delay differential equations)}\\ \|z\|=1,&\quad\text{(amplitude condition of scaled orbit)}\\ \int\langle z(s),\hat{z}^{\prime}(s)\rangle ds=0,&\quad\text{(anchor condition)}\\ f(x,\dots,x,\alpha,\beta)=0,&\quad\text{($x$ is a steady state)}\\ \theta_{BC}(z(0),x,a,\alpha,\beta,\eta)=0.&\quad\text{(embedding boundary condition)}\end{cases}$$

(20)

At this stage, the period $\omega$ of the periodic orbit is implicit. In passing to the spectral representation, we will make it explicit. Define the frequency $\psi=\frac{2\pi}{\omega}$ and expand $z$ in Fourier series:

$\displaystyle z(t)$

$\displaystyle=\sum_{k\in\mathbb{Z}}z_{k}e^{ik\psi t}.$

(21)

Recall that for a real (as opposed to complex-valued) periodic orbit, $z_{k}\in\mathbb{C}^{n+m}$ will satisfy the symmetry $z_{k}=\overline{z_{-k}}$. To substitute (21) into the differential equation in (20), we must examine how time delays transform under Fourier series. Observe

$$z(t+\mu)=\sum_{k\in\mathbb{Z}}z_{k}e^{ik\psi\mu}e^{ik\psi t},$$

which means that at the level of Fourier coefficients, a delay of $\mu$ corresponds to a complex rotation

$\displaystyle z_{k}\mapsto(\zeta_{\mu}(\psi)z)_{k}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,e^{ik\psi\mu}z_{k}.$

(22)

Note that this operator is linear on $C^{\mathbb{Z}}_{n+m}$ and bounded on $\ell_{\nu}^{1}(\mathbb{C}^{n+m})$. To densify our notation somewhat, we define $\zeta(\psi):\ell_{\nu}^{1}(\mathbb{C}^{n+m})\rightarrow\ell_{\nu}^{1}(\mathbb{C}^{n+m})^{J}$ by

$$\zeta(\psi)z=(\zeta_{\mu_{1}}(\psi)z,\dots,\zeta_{\mu_{J}}(\psi)z).$$

As for the derivative, we make use of the operator $(Kz)_{k}=kz_{k}$. Since $\tilde{f}$ is polynomial, subsituting (21) into the first equation of (20) will result in an equation of the form

$$\psi iKz=\mathbf{f}(\zeta(\psi)z,x,a,\alpha,\beta,\eta)$$

for a function $\mathbf{f}:(\mathbb{C}^{\mathbb{Z}}_{n+m})^{J}\times\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{2}\times\mathbb{R}^{m}\rightarrow(\mathbb{C}^{\mathbb{Z}}_{n+m})$ being a (formal) vector polynomial with respect to Fourier convolution, in the arguments $(\mathbb{C}^{\mathbb{Z}}_{n+m})^{J}$. Observe that we have abused notation and identified the function $z$ in (21) with its sequence of Fourier coefficients. As an example, the nonlinearity $z\mapsto z(t)^{2}z(t+\mu)$ is transformed in Fourier to the nonlinearity

$$z\mapsto(z*z)*(\zeta_{\mu}(\psi)z).$$