aa r X i v : . [ m a t h . D S ] J a n A proof of Jones’ conjecture
Jonathan Jaquette ∗† October 8, 2018
Abstract
In this paper, we prove that Wright’s equation y ′ ( t ) = − αy ( t − { y ( t ) } hasa unique slowly oscillating periodic solution for parameter values α ∈ ( π , . α > π . Furthermore, there are noisolas of periodic solutions to Wright’s equation; all periodic orbits arise from Hopfbifurcations. Key words.
Wright’s Equation · Jones’ Conjecture · Delay Differential EquationsComputer-Assisted Proofs · Branch and Bound · Krawczyk method
Contents ∗ Partially supported by NSF DMS 0915019, NSF DMS 1248071 † Department of Mathematics, Hill Center-Busch Campus, Rutgers, The State University of New Jersey,Piscataway, NJ, USA, 08854-8019. [email protected] Introduction
An often studied class of delay differential equations are negative feedback systems of theform: x ′ ( t ) = − αf ( x ( t − xf ( x ) > x = 0 and f ′ (0) >
0. One particularly well studied example of (1) iswhen f ( x ) = e x −
1, better known as Wright’s equation, which after making the change ofvariables y = e x − y ′ ( t ) = − α y ( t −
1) [1 + y ( t )] . (2)In [11], Jones proved that for α > π there exists at least one slowly oscillating periodicsolution (SOPS) . That is, a periodic solution y : R → R which is positive for at least oneunit of time (the delay time in Wright’s equation), negative for at least one unit of time,and then repeats. In this paper we prove there is a unique SOPS to (2) for α ∈ ( π , . Theorem 1.1 (Jones’ conjecture) . For every α > π there exists a unique slowly oscillatingperiodic solution to (2) . Amplitude ofperiodic orbits π /2 α π /2 9 π /2 Figure 1:
A bifurcation diagram for periodic so-lutions to Wright’s equation. There are no folds inthe principal branch of slowly oscillating periodic so-lutions (solid curve). While there may be folds inthe branches of rapidly oscillating periodic solutions(dotted curves), it is conjectured that this does notoccur. There are no isolas of periodic solutions (notdisplayed).
This work contributes a capstone tomany decades of mathematical work study-ing Wright’s equation. To briefly review,a principal branch of slowly oscillating pe-riodic orbits is born at α = π and con-tinues on for all α > π [22]. Moreover,Wright’s equation has supercritical Hopf bi-furcations at α = π +2 nπ for integers n ≥ n = 0, and rapidly oscillating pe-riodic orbits arising when n ≥ α ,the collection of periodic orbits forms a 2-dimensional manifold [23].A two-part geometric version of Jones’conjecture was proposed in [12]: (i) theprincipal branch of SOPS does not fold backon itself, and (ii) there are no other con-nected components ( isolas ) of SOPS. By[10, 12, 26, 29] the principal branch does nothave any folds α > π .In [10, 29] it is shown that there is aunique SOPS for α ≥ .
9. These proofs use that fact that if every SOPS is asymptoti-cally stable for some α > π , then there is a unique SOPS [30]. Using estimates describingSOPS for when α is large [20], Xie showed that there is a unique SOPS for all α ≥ .
67 [29].By using computer-assisted proofs to characterize SOPS to Wright’s equation [10], thismethod was extended to show there is a unique SOPS for α ∈ [1 . , . α close to the bifurcation value π the dynamics becomes center-like, andproving uniqueness through these stability arguments becomes infeasible. To overcome this2bstacle, we equate the problem of finding periodic orbits to (2) with a zero–finding problemin a space of Fourier coefficients. We then employ rigorous numerics to derive a computer-assisted proof that there is a unique SOPS to Wright’s equation for α ∈ ( π , . x ( t − Theorem 1.2.
Let x be a periodic solution to (1) at parameter α with period L and lapnumber N . Then there exists a SOPS x ( t ) = x ( rt ) to (1) at parameter α = rα where r := 1 − N − L . Thus, every periodic orbit is on a branch originating from one of the Hopf bifurcationsat α = π + 2 nπ . That is to say, there are no isolas of rapidly oscillating periodic solutions.However, this is not sufficient to show there are no folds in the branches of rapidly oscillatingperiodic solutions. The proofs for Theorem 1.1 and Theorem 1.2 are presented at the endof Section 6, and we discuss future directions in Section 7. In this paper we show that there is a unique slowly oscillating periodic orbit to (2) for all α ∈ ( π , . F defined in a space of Fourier coefficients(see Section 2.2). Since periodic solutions to (2) must have a high degree of smoothness, inparticular real analyticity [21, 28], their Fourier coefficients will decay very rapidly. That isto say, the functional we are interested in can be well approximated by a Galerkin projectiononto a finite number of Fourier modes.In finite dimensions, there are efficacious techniques for rigorously locating and enumer-ating the solutions to a system of nonlinear equations by way of interval arithmetic [8,18,19].We apply these techniques in infinite dimensions, specifically the branch and bound method,also referred to as a branch and prune method. That is, we first construct a bounded set X of Fourier coefficients which contains all the zeros of F (see Section 5). Then we partition X into a finite number of pieces { X n } which we refer to as cubes (see Definition 2.7). Foreach cube X n we are interested to know whether:(a) there exists a unique point ˆ x ∈ X n for which F (ˆ x ) = 0, or(b) there does not exist any points ˆ x ∈ X n for which F (ˆ x ) = 0.If we can show that ( a ) holds for one cube, and ( b ) holds for all the other cubes, then wewill have shown that F = 0 has a unique solution.This approach requires some additional preparation. Since periodic orbits to (2) form a2-manifold in phase space [23], the functional F we construct in Section 2.2 will not haveisolated zeros. The numerical techniques we employ are suited to finding isolated zeros, soit is necessary to reduce the dimension of the kernel by two. Along the principal branch α can be taken as one of the coordinate dimensions. We reduce this dimension by treating α as a parameter and performing our estimates uniformly in α . The other dimension can be3ttributed to time translation; if y ( t ) is a periodic orbit, then so is y ( t + τ ) for any τ ∈ R .We reduce this dimension by imposing a phase condition; we may assume without loss ofgenerality that the first Fourier coefficient is a positive real number (see Proposition 5.4).The central technique we use to determine whether ( a ) or ( b ) holds for a given cube is theKrawczyk method [8, 17–19]. For a function f ∈ C ( R n , R n ) the Krawczyk operator takesas input a rectangular set X ⊆ R n and produces as output a rectangular set K ( X ) ⊆ R n .This set K ( X ) has the properties that, ( i ) if K ( X ) ⊆ X , then there exists a unique pointˆ x ∈ X for which f (ˆ x ) = 0, and ( ii ) if ˆ x ∈ X and f (ˆ x ) = 0, then ˆ x ∈ K ( X ). Clearly ( i )implies ( a ), and if X ∩ K ( X ) = ∅ then ( b ) follows. Additionally, even if we can prove neither( a ) nor ( b ) our situation could still improve; we can replace X X ∩ K ( X ) without losingany solutions.Adjustments are needed to generalize the Krawczyk operator to infinite dimensionalsystems. In [7] a Krawczyk operator is defined in Hilbert space to study fixed points andperiod-2 orbits in an infinite dimensional map. In Section 2.1 we present a generalizationof the Krawczyk operator to Banach spaces. a Figure 2:
The main result of this paper is a collection of “cubes”in Fourier space which cover the Fourier coefficients of SOPS to(2). The first Fourier coefficient of this cover is plotted here withrespect to α . Inside each green cube there exists a unique SOPScorresponding to each α , essentially by Theorem 2.2. Inside eachblue cube the only SOPS that can exist are on the principalbranch, by [26]. To determine whether ( a ) or( b ) holds the Krawczyk operatorby itself is not always sufficient,and we combine several additionaltests to create a single pruning op-erator (see Section 4). One prob-lem is that y ≡ k F k is bounded away from zero on a cube X n , then ( b ) holds.Algorithm 6.1 follows the standard format of a global branch and bound method. Inshort, for a collection of cubes we successively prune each of its cubes. If ( a ) holds for agiven cube, then it is set aside and added to a list of solutions. If ( b ) holds for a given cube,then that cube is discarded. If the pruning operator significantly reduces the size of a cube,then the pruning operator is applied again. If none of these are the case, then the cubeis split in half, and both pieces are added back to the collection of cubes to inspect. Thisprocess repeats until all of the cubes have been removed or reduced to a sufficiently smallsize.The output of Algorithm 6.1 is three collections of cubes: A , B , and R (see Figure 2).In Theorem 6.2 we show that these sets have the properties that, ( i ) each cube in A has a4nique solution with respect to α , ( ii ) the cubes in B are near the Hopf bifurcation, withany solutions contained therein residing on the principal branch, and ( iii ) all solutions to F = 0 are contained in S A ∪ B ∪ R .Ideally R = ∅ , and this will often be the case if the zeros of F are simple and thealgorithm is allowed to run a sufficiently long time. However we are trying to verify not justsimple, isolated solutions, but a 1-parameter family of solutions. As such, sometimes whena cube is split in two this division will bisect the curve of solutions (see Figure 5). Whenthis occurs the algorithm will be forced to subdivide many cubes near where the solutioncurve was bisected, resulting in the variably sized cubes noticeable in Figure 2. To addressthis we recombine the cubes in R which have the same α values, then subsequently use theKrawczyk operator to show that ( a ) holds on the recombined cubes (see Algorithm 6.3). Inthis fashion, we prove Theorem 1.1. In numerical analysis there are many variations on the theme of Newton’s method: x n +1 x n − Df ( x n ) − f ( x n ) . As inverting a matrix is computationally expensive, one alternativemethod is to replace DF ( x n ) − with a fixed matrix A † ≈ Df ( x ) − . If f ( x ) ≈
0, then theNewton-Kantorovich theorem gives conditions for when the map T ( x ) = x − A † f ( x ) definesa contraction map in a neighborhood about x . The Krawczyk operator may be thoughtof as a way of bounding the image of T , itself being defined on rectangular sets X ⊆ R n and having the property that T ( X ) ⊆ K ( X, x ). Rectangular, in the sense that X can begiven as the product of intervals in the coordinate directions of R n . Here we generalize theKrawczyk operator to non-rectangular subsets of Banach spaces. Definition 2.1.
Let
Y, Z denote Banach spaces and let A † : Z → Y be an injective, boundedlinear operator. Fix a convex, closed and bounded set X ⊆ Y , a neighborhood U ⊇ X , anda Frechet differentiable function f : U → Z . Let ( I − A † Df ( X ))( X − ¯ x ) = conv [ x ,x ∈ X ( I − A † Df ( x ))( x − ¯ x ) , where conv denotes the closure of the convex hull. For a point ¯ x ∈ X we define the Krawczykoperator K ( X, ¯ x ) as: K ( X, ¯ x ) := ¯ x − A † f (¯ x ) + ( I − A † Df ( X ))( X − ¯ x ) ⊆ Y. (3)Typically ¯ x is taken to be the center of X , and A † is taken to be an approximate inverse of DF (¯ x ). If K ( X, ¯ x ) ⊆ X for a rectangular set X ⊆ R n , then there exists a unique ˆ x suchthat f (ˆ x ) = 0. In Theorem 2.2 we prove an analogous result. The existence of a fixed pointis achieved by the Schauder fixed point theorem. However to prove uniqueness, droppingthe rectangular condition causes problems even in finite dimensions; in Theorem 2.2 ( iv ) weprescribe a hypothesis sufficient for proving uniquessness in our level of generality. Theorem 2.2.
Suppose K is a Krawczyk operator as given in Definition 2.1 and T := x − A † f ( x ) .(i) If x ∈ X , then T ( x ) ∈ K ( X, ¯ x ) .(ii) If ˆ x ∈ X and f (ˆ x ) = 0 , then ˆ x ∈ K ( X, ¯ x ) . iii) If K ( X, ¯ x ) ⊆ X and X is compact, then there exists a point ˆ x ∈ X such that f (ˆ x ) = 0 .(iv) If K ( X, ¯ x ) ⊆ X and there exists ≤ λ < such that ( I − A † Df ( X ))( X − ¯ x ) ⊆ λ · ( X − ¯ x ) , then there exists a unique point ˆ x ∈ X such that f (ˆ x ) = 0 .Proof. (i) Fix a point x ∈ X and write h = x − ¯ x . By the mean-value theorem for Frechetdifferentiable functions [1], we have: T ( x ) = ¯ x − A † f (¯ x ) + Z DT (¯ x + th ) · h dt = ¯ x − A † f (¯ x ) + lim N →∞ N X i =1 1 N (cid:0) I − A † Df (¯ x + iN h ) (cid:1) · h ∈ ¯ x − A † f (¯ x ) + conv (cid:0)(cid:0) I − A † Df ( X ) (cid:1) · ( x − ¯ x ) (cid:1) ⊆ K ( X, ¯ x ) . (ii) If there is some ˆ x ∈ X such that f (ˆ x ) = 0, then ˆ x = T (ˆ x ) ∈ K ( X, ¯ x ).(iii) Since T ( X ) ⊆ K ( X, ¯ x ) by (i) and K ( X, ¯ x ) ⊆ X by assumption, therefore T ( X ) ⊆ X .As T is continuous and X is convex and compact, then by the Schauder fixed pointtheorem there exists some ˆ x ∈ X such that ˆ x = T (ˆ x ). Since A is injective, the zerosof f are in bijective correspondence with the fixed points of T , thereby f (ˆ x ) = 0.(iv) Inductively define: X = X , x = ¯ x , and X n +1 = T ( X n ), x n +1 = T ( x n ). Note thatas T ( X ) ⊆ X then X n +1 ⊆ X n for all n . We show that X n ⊆ x n + λ n ( X − x ). Thisis clearly true for n = 0. For n ≥ X n +1 ⊆ K ( X n , x n )= x n − A † f ( x n ) + ( I − A † Df ( X n )) · ( X n − x n ) ⊆ x n +1 + ( I − A † Df ( X )) · λ n ( X − x ) ⊆ x n +1 + λ n +1 ( X − x ) . Since λ n k X − x k can be made arbitrarily small and { x n } ∞ n = N ⊆ X N , it follows that { x n } is a Cauchy sequence. As X is complete, then lim x n = ˆ x and additionally T ∞ n =0 X n = ˆ x . Thereby ˆ x is the unique fixed point of T in X = X and the uniquezero of f in X . As in [12,26], we convert Wright’s equation into a functional equation on the space of Fouriercoefficients. For a continuous periodic function y : R → R with frequency ω >
0, we maywrite it as: y ( t ) = X k ∈ Z c k e iωkt (4)where c k ∈ C and P k ∈ Z | c k | < ∞ . By [26] it suffices to work with sequences { c k } ∞ k =1 tostudy periodic solutions to (2). This is because real-valued functions have Fourier coefficients6atisfying c − k = c ∗ k , and periodic solutions to (2) necessarily satisfy c = 0. Hence we definethe following Banach spaces: ℓ := ( { c k } ∞ k =1 : c k ∈ C and ∞ X k =1 | c k | < ∞ ) k c k ℓ =2 ∞ X k =1 | c k | (5)Ω s := (cid:26) { c k } ∞ k =1 : c k ∈ C and sup k ∈ N k s | c k | < ∞ (cid:27) k c k s = sup k ∈ N k s | c k | . (6)The smoother a function is the faster its Fourier coefficients will decay; if a function is s –times continuously differentiable, then its Fourier coefficients will be in Ω s . Since periodicsolutions to (2) are real analytic [21, 28], it follows that their Fourier coefficients will be inΩ s for all s ≥ y is a solution to Wright’s equation, then by substituting (4) into (2) we obtain: X k ∈ Z iωkc k e iωkt = − α X k ∈ Z c k e − iωk e iωkt ! X k ∈ Z c k e iωkt ! . (7)By matching the e iωkt terms, subtracting the RHS, and dividing through by α , we obtainthe following sequence of equations for k ∈ Z below:[ F ( α, ω, c )] k = (cid:0) i ωα k + e − iωk (cid:1) c k + X k ,k ∈ Z k + k = k e − iωk c k c k (8)= (cid:0) i ωα k + e − iωk (cid:1) c k + k − X j =1 e − iωj c j c k − j + ∞ X j =1 (cid:16) e − iω ( j + k ) + e iωj (cid:17) c ∗ j c j + k . (9)Dividing through by α ensures that the parameter dependence in F is solely concentratedin the linear part. In this manner y is a periodic solution with frequency ω to Wright’sequation at parameter α if and only if [ F ( α, ω, c )] k = 0 for all k ∈ Z [10, 26].To more succinctly express the functional F we introduce additional notation. For asequence c = { c k } ∞ k =1 we denote the projection onto the k -coefficient by [ c ] k := c k . Wedefine unnormalized basis elements e j ∈ ℓ , Ω s for j ∈ N by:[ e j ] k = ( k = j, k = j. We define the discrete convolution a ∗ b for a, b ∈ ℓ component-wise by:[ a ∗ b ] k := X | k | + | k | = k a k b k = k − X j =1 a j b k − j + ∞ X j =1 a ∗ j b k + j + a k + j b ∗ j , where a − k = a ∗ k and b − k = b ∗ k , and the sum is taken over k , k ∈ Z . The space ℓ is aBanach algebra, which is to say that k a ∗ b k ℓ ≤ k a k ℓ k b k ℓ for all a, b ∈ ℓ . While Ω s is not a Banach algebra per se , if s ≥ B ≥ k a ∗ b k s ≤ B k a k s k b k s for all a, b ∈ Ω s (see [12, 27]). Lastly, we define a linear operator K : Ω s → Ω s +1 and a continuous family of linear operators U ω : Ω s → Ω s − as below:[ K c ] k := c k /k, [ U ω c ] k := e − ikω c k . U ω is necessary for its continuity, as ∂∂ω U ω = − i K − U ω .We may extend U ω to act on bi-infinite sequences { c k } k ∈ Z using the same component-wisedefinition. Additionally, this extension is compatible with our definition of the discreteconvolution, as [ U ω c ] ∗ k = [ U ω c ] − k whenever c ∗ k = c − k . In Definition 2.3 we rewrite (8) inoperator notation and list several propositions, the proofs of which are left to the reader. Definition 2.3.
Define the function F : R × Ω s → Ω s − as: F ( α, ω, c ) := ( i ωα K − + U ω ) c + ( U ω c ) ∗ c. (10) Proposition 2.4 (Theorem 2.2 in [26]) . Let α, ω > . If c ∈ ℓ solves F ( α, ω, c ) = 0 , then y ( t ) , given by (4) with c = 0 and c − k = c ∗ k , is a periodic solution of (2) with period π/ω .Vice versa, if y ( t ) is a periodic solution of (2) with period π/ω , then its Fourier coefficientssatisfy c = 0 , c − k = c ∗ k , { c k } ∞ k =1 ∈ ℓ and solve F ( α, ω, { c k } ∞ k =1 ) = 0 . Proposition 2.5.
For each α > and s ≥ the function F : R × Ω s → Ω s − is Frechetdifferentiable, with partial derivatives given as: ∂∂ω F ( α, ω, c ) = i K − ( α − I − U ω ) c − i ( K − U ω c ) ∗ c (11) ∂∂c F ( α, ω, c ) · h = ( i ωα K − + U ω ) h + ( U ω c ) ∗ h + ( U ω h ) ∗ c, (12) where h ∈ Ω s . Proposition 2.6.
Define γ ( k, n ) := e − iω ( n + k ) + e iωn and γ ( k, n ) := e − iωn + e iω ( n − k ) .Writing c k = a k + ib k , the component-wise derivatives of F are given as: ∂∂ω [ F ( α, ω, c )] k = ik ( α − − e − iωk ) c k − i k − X j =1 je − iωj c j c k − j − i ∞ X j =1 (cid:16) ( j + k ) e − iω ( j + k ) − je iωj (cid:17) c ∗ j c j + k .∂∂a n [ F ( α, ω, c )] k = ( i ωα k + e − iωk ) + ( γ c n + k + γ c k − n if ≤ n < kγ c n + k + γ c ∗ n − k if k ≤ n. i ∂∂b n [ F ( α, ω, c )] k = ( i ωα k + e − iωk ) + ( − γ c n + k + γ c k − n if ≤ n < k − γ c n + k + γ c ∗ n − k if k ≤ n. By working in a space of rapidly decaying Fourier coefficients, we are able to closely approx-imate the value of F using a Galerkin projection. Since F : R × Ω s → Ω s − has distinctdomain and range, we need to define two sets of projection maps. We define projectionmaps π α , π ω : R × Ω s → R and π c : R × Ω s → Ω s on points x = (˜ α, ˜ ω, ˜ c ) ∈ R × Ω s as: π α ( x ) := ˜ α π ω ( x ) := ˜ ω π c ( x ) := ˜ c. (13)For a fixed integer M ∈ N , define the projection maps π M , π ∞ : Ω s → Ω s by: π M ( c ) := M X k =1 [ c ] k e k π ∞ ( c ) := c − π M ( c ) . (14)8efine the projection maps π ′ M , π ′∞ : R × Ω s → R × Ω s by: π ′ M ( c ) := ( π α ( x ) , π ω ( x ) , π M ◦ π c ( x )) π ′∞ ( c ) := (0 , , π ∞ ◦ π c ( x )) . (15)For any bounded set X ⊆ R × Ω s , define: | X | k := sup x ∈ X | [ π c ( x )] k | . We define for F its Galerkin projection and remainder F M , F ∞ : R × Ω s → Ω s − as follows: F M ( x ) := π M ◦ F ( π ′ M ( x )) , F ∞ ( x ) := F ( x ) − F M ( x ) . (16)By construction F = F M + F ∞ .To show that there is a unique SOPS to (2) we need to evaluate F not just on singlepoints but on voluminous subsets of its domain. The central subset of R × Ω s we considerin this paper are cubes which we define as follows: Definition 2.7.
For M ∈ N , s ≥ , C > define a cube X := X M × X ∞ ⊆ R × Ω s to beof the following form: X M := [ α, α ] × [ ω, ω ] × M Y k =1 [ A k , A k ] × [ B k , B k ] (17) X ∞ := { c k ∈ C : | c k | ≤ C /k s } ∞ k = M +1 . (18)To denote the union of a collection of cubes S := { X i ⊆ R × Ω s } we define S S := S X ∈S X ⊆ R × ˜Ω s .There are primarily two reasons we have chosen to consider cubical subsets of R × Ω s .Firstly, cubes are particularly easy to refine into smaller pieces. This is useful because tobegin using a branch and bound method, we need to obtain global bounds on the solutionspace, and then partition these bounds into smaller pieces. In practice, we reduce the sizeof a cube by either subdividing it along a lower dimension into two cubes, or replacing thecube by its intersection with the Krawczyk operator: X X ∩ K ( X, ¯ x ). In both thesecases the resulting object is again a cube. In this manner, we can use cubes to cover thesolutions to F = 0, and then refine the cover using successively smaller cubes.Secondly, cubes facilitate explicit computations of F M and analytical estimates of F ∞ .While formally F M is an infinite dimensional map, computationally, we may consider F M to be a map R × C M → C M . To calculate F M , we simply truncate the second sum in (9)at j = M − k . As the π ′ M projection of a cube is given as a finite product of intervals,it is well suited for using interval arithmetic [18] to bound the image of F M ( X ). On theother hand, bounding F ∞ requires significantly more analysis. Below is a simple, yet everrecurring estimate in our calculations: ∞ X k = M +1 k s ≤ Z ∞ M x s dx = 1( s − M s − , (19)where we take s >
1. For example, if a cube X ⊆ R × Ω s satisfies s >
1, then k π c x k ℓ ≤ P Mk =1 | X | k + C ( s − M s − for all x ∈ X . This specific bound on the ℓ norm is later used inAlgorithm 4.1 to check whether Lemmas 2.8 or 2.9 apply.9 emma 2.8 (Theorems E.1 and E.2 in [26]) . Let ω ≥ . , α ∈ (0 , , and define g ( α, ω ) := q(cid:0) − ωα (cid:1) + 2 ωα (1 − sin ω ) . (20) If F ( α, ω, c ) = 0 , then either c ≡ or g ( α, ω ) ≤ k c k ℓ . Lemma 2.9 (Theorem 4.10 [26]) . For each α ∈ ( π , π + 0 . there is a unique (up totime translation) periodic solution to Wright’s equation with Fourier coefficients satisfying k c k ℓ ≤ . and having frequency | ω − π | ≤ . . We note that while Lemma 2.8 is stated only for ω ≥ . α ∈ (0 , k c k ℓ as opposed to a bound on k y ′ k L as in the original paper. This allows us to use thestronger result derived in the proof of [26, Theorem 4.10], namely that the solution exists and is unique, as opposed to the exact result stated in [26, Theorem 4.10], which is thatthere is most one periodic solution.The remainder of this section is dedicated to proving Lemma 2.12, which estimates F ∞ ,its derivatives, and convolution products resulting from points inside of a cube. Theseestimates are used in Definition 3.2 to construct an outer approximation to the Krawczykoperator. The reader is encouraged to skip the proof of Lemma 2.12 on a first reading,which is best summarized as bounding various infinite sums by various finite sums and theestimate in (19). These bounds are presented in Definition 2.11, all of which are given as afinite number of operations, explicitly computable in terms of C and the π ′ M -projection ofa given cube. In Lemma 2.10 we define the constant γ M which is needed for the definitionof (26). Lemma 2.10 (Lemma 24 [27]) . Let s ≥ and let s ∗ be the largest integer such that s ∗ ≤ s and define: γ k := 2 (cid:20) kk − (cid:21) s + (cid:20) k − k + π − (cid:21) (cid:20) k + 12 (cid:21) s ∗ − . For k ≥ , we have that P k − k =1 k s k s ( k − k ) s ≤ γ k . If ≤ M ≤ k , then γ k ≤ γ M . Definition 2.11.
Fix a cube X with s > , define C := sup x ∈ X k π c x k s , and select a point ¯ x = (¯ α, ¯ ω, ¯ c ) ∈ X such that ¯ x = π ′ M (¯ x ) . Define H = X − ¯ x , and define ∆ ω ∈ R such that ∆ ω ≥ sup x ∈ H | π ω ( x ) − ¯ ω | .Define h, g iM , g iiM to be functions of the form g M : X g M ( X ) ∈ R M and define g i ∞ , g ii,a ∞ , g ii,b ∞ to be functions of the form g ∞ : X g ∞ ( X ) ∈ R as follows: h ( X )] k := 2 C ( s − M s − ( M + k + 1) s + 2 C M X j = M − k +1 | X | j ( j + k ) s (21)[ g iM ( X )] k := 2 C ∆ ω M X j = M − k +1 | X | j ( j + k ) ( s − + C ∆ ω ( s − M + k + 1) s M ( s − + C ∆ ω ( s − M + k + 1) ( s − M ( s − (22)[ g iiM ( X )] k := 4 C ( s − M + k + 1) s M s − + 2 C M X j = M − k +1 | H | j ( j + k ) s (23) g i ∞ ( X ) := max M +1 ≤ k ≤ M k s M X j = k − M | ¯ c j ¯ c k − j | (24) g ii,a ∞ ( X ) := max M +1 ≤ k ≤ M k s M X j = k − m | H | j | X | k − j + 2 C (2 s + 1)( s − M s − + C M X j =1 ( | X | j + | H | j ) (cid:18)(cid:18) M + j + 1 M + 1 (cid:19) s + 1 (cid:19) (25) g ii,b ∞ ( X ) := C γ M +1 C C (cid:18) s − M + 2)( s −
2) + ss − (cid:19) . (26) Lemma 2.12.
Fix a cube X with M ≥ , s > , a point ¯ x ∈ X such that ¯ x = π ′ M (¯ x ) , anddefine H = X − ¯ x . Then the following inequalities hold: sup x ∈ X | F ∞ ( x ) | k < [ h ( X )] k ≤ k ≤ M (27)sup x ∈ X,h ∈ H (cid:12)(cid:12) ∂∂ω F ∞ ( x ) · π ω ( h ) (cid:12)(cid:12) k ≤ [ g iM ( X )] k ≤ k ≤ M (28)sup x ∈ X,h ∈ H (cid:12)(cid:12) ∂∂c F ∞ ( x ) · π c ( h ) (cid:12)(cid:12) k ≤ [ g iiM ( X )] k ≤ k ≤ M (29) | F ∞ (¯ x ) | k ≤ k s g i ∞ ( X ) M + 1 ≤ k (30)sup x ∈ X,h ∈ H | π c ( h ) ∗ π c ( x ) | k ≤ k s g ii,a ∞ ( X ) M + 1 ≤ k (31)sup x ,x ∈ X (cid:12)(cid:12) ( K − π c ( x )) ∗ π c ( x ) (cid:12)(cid:12) k ≤ k s − g ii,b ∞ ( X ) M + 1 ≤ k. (32)Throughout, let us write X M = π ′ M ( X ), H M = π ′ M ( H ), and H ∞ = π ′∞ ( H ), noting alsothat H ∞ = π ′∞ ( X ). Proof of (27) . We show that | F ∞ ( x ) | k < [ h ( X )] k for 1 ≤ k ≤ M and all x ∈ X . Fix x = ( α, ω, c ) ∈ X , and write c M = π M ( c ) and c ∞ = π ∞ ( c ). We compute: π M ◦ F ∞ ( x ) = π M ◦ ( F ( x ) − F ( π ′ M x ))= π M ◦ (( U ω c ) ∗ c − ( U ω c M ) ∗ c M )= π M ◦ (( U ω c M ) ∗ c ∞ + ( U ω c ∞ ) ∗ c M + ( U ω c ∞ ) ∗ c ∞ )11ince | U ω c | k = | c | k , it follows that for 1 ≤ k ≤ M we compute the estimate below: | ( U ω c M ) ∗ c ∞ | k + | ( U ω c ∞ ) ∗ c M | k ≤ ∞ X j =1 | c ∗ M | j | c ∞ | k + j + | c M | k + j | c ∗∞ | j =2 M X j = M − k +1 | c ∗ M | j | c ∞ | j + k ≤ M X j = M − k +1 | X | j C ( j + k ) s . The last estimate uses the property that | c j | ≤ C /j s for j ≥ M + 1.We calculate ( U ω c ∞ ) ∗ c ∞ as below, again using | c j | ≤ C /j s for j ≥ M + 1. | ( U ω c ∞ ) ∗ c ∞ | k ≤ ∞ X j = M +1 | c ∗∞ | j | c ∞ | k + j + | c ∞ | j + k | c ∗∞ | j ≤ ∞ X j = M +1 C j s ( j + k ) s ≤ C ( s − M s − ( M + k + 1) s . Hence for 1 ≤ k ≤ M , it follows that: | F ∞ ( x ) | k ≤ C ( s − M s − ( M + k + 1) s + 2 C M X j = M − k +1 | X | j ( j + k ) s = [ h ( X )] k . Proof of (28) . We show that (cid:12)(cid:12) ∂∂ω F ∞ ( x ) · π ω ( h ) (cid:12)(cid:12) k ≤ [ g iM ( X )] k for 1 ≤ k ≤ M and all x ∈ X and h ∈ H . Select some x = ( α, ω, c ) ∈ X and write c M = π M ( c ) and c ∞ = π ∞ ( c ). From(11) we can calculate ∂∂ω F ∞ ( x ) as follows: ∂∂ω F ∞ ( x ) = − i ( K − U ω c ) ∗ c + iπ M ( K − U ω c M ) ∗ c M = − iπ ∞ (cid:0) K − U ω c M (cid:1) ∗ c M − i (cid:0) K − U ω c M (cid:1) ∗ c ∞ − i (cid:0) K − U ω c ∞ (cid:1) ( c M + c ∞ ) . Hence, for 1 ≤ k ≤ M we may calculate the following: (cid:12)(cid:12) ∂∂ω F ∞ ( x ) (cid:12)(cid:12) k ≤ sup c M ∈ X M ; c ∞ ,c ′∞ ∈ H ∞ (cid:12)(cid:12) ( K − c M ) ∗ c ∞ + ( K − c ∞ ) ∗ c M + ( K − c ∞ ) ∗ c ′∞ (cid:12)(cid:12) k . (33)For 1 ≤ k ≤ M and any c M ∈ X M , c ∞ ∈ H ∞ we can simplify the first two summands in(33) as follows:( K − c M ) ∗ k c ∞ = ∞ X j =1 [ K − c ∗ M ] j [ c ∞ ] k + j + [ K − c M ] k + j [ c ∗∞ ] j = ∞ X j = M +1 − k j [ c ∗ M ] j [ c ∞ ] k + j ( K − c ∞ ) ∗ k c M = ∞ X j =1 [ K − c ∗∞ ] j [ c M ] k + j + [ K − c ∞ ] k + j [ c ∗ M ] j = ∞ X j = M +1 − k ( k + j )[ c ∞ ] k + j [ c ∗ M ] j . K − c M ) ∗ k c ∞ + ( K − c ∞ ) ∗ k c M = M X j = M − k +1 (2 j + k )[ c ∞ ] j + k [ c ∗ M ] j (cid:12)(cid:12) ( K − c M ) ∗ c ∞ (cid:12)(cid:12) k + (cid:12)(cid:12) ( K − c ∞ ) ∗ c M (cid:12)(cid:12) k ≤ M X j = M − k +1 (2 j + k ) C ( j + k ) s | X | j ≤ C M X j = M − k +1 | X | j ( j + k ) s − . (34)Again, we used the estimate | c j | ≤ C /j s for j ≥ M + 1. We estimate the third summandin (33) for c ∞ , c ′∞ ∈ H ∞ as follows:( K − c ∞ ) ∗ k c ′∞ = ∞ X j = M +1 j [ c ∗∞ ] j [ c ′∞ ] k + j + ( j + k )[ c ∞ ] j + k [ c ′∞∗ ] j (cid:12)(cid:12) ( K − c ∞ ) ∗ c ′∞ (cid:12)(cid:12) k ≤ ∞ X j = M +1 C j ( s − ( j + k ) s + C j s ( j + k ) ( s − ≤ C ( s − M + k + 1) s M ( s − + C ( s − M + k + 1) ( s − M ( s − . (35)By combining the estimates from (34) and (35) into (33), and recalling our choice of ∆ ω inDefinition 2.11, then for 1 ≤ k ≤ M we obtain the following:sup x ∈ X,h ∈ H (cid:12)(cid:12) ∂∂ω F ∞ ( x ) · π ω ( h ) (cid:12)(cid:12) k ≤ C ∆ ω M X j = M − k +1 | X | j ( j + k ) ( s − + C ∆ ω ( s − M + k + 1) s M ( s − + C ∆ ω ( s − M + k + 1) ( s − M ( s − = [ g iM ( X )] k . Proof of (29) . We show that (cid:12)(cid:12) ∂∂c F ∞ ( x ) · π c ( h ) (cid:12)(cid:12) k ≤ [ g iiM ( X )] k for 1 ≤ k ≤ M and all x ∈ X and h ∈ H . Let ( α, ω, c ) ∈ X and h ∈ π c ( H ). From (12) we calculate ∂∂c ( F ( X ) − F M ( X )) · h below: ∂∂c ( F ( x ) − F ( π ′ M x )) · h = (( U ω h ) ∗ c + ( U ω c ) ∗ h ) − (( U ω h ) ∗ c M + ( U ω c M ) ∗ h )=( U ω h ) ∗ ( c − c M ) + ( U ω ( c − c M )) ∗ h. Since c − c M ∈ H ∞ , it follows that: | ∂∂c [ F ( x ) − F ( π ′ M x )] · h | k ≤ sup h ∈ H,h ′ ∈ H ∞ · | h ∗ h ′∞ | k . For h ∈ H and h ′ ∈ H ∞ and for 1 ≤ k ≤ M , we calculate h ∗ k h ′ below, using the property13hat [ h ′ ] j = 0 for j ≤ M . h ∗ k h ′ = ∞ X j =1 [ h ∗ ] j [ h ′ ] k + j + [ h ] k + j [ h ′∗ ] j = M X j = M − k +1 [ h ∗ ] j [ h ′ ] k + j + ∞ X j = M +1 [ h ∗ ] j [ h ′ ] k + j + [ h ] k + j [ h ′∗ ] j . By applying the estimates | h j | ≤ | H | j for j ≤ M , and | h | j , | h ′ | j ≤ C /j s for j ≥ M + 1, weobtain the following: (cid:12)(cid:12) ∂∂c F ∞ ( x ) · h (cid:12)(cid:12) k ≤ M X j = M − k +1 | H | j C ( j + k ) s + ∞ X j = M +1 C j s ( j + k ) s ≤ C M X j = M − k +1 | H | j ( j + k ) s + 4 C ( s − M + k + 1) s M s − = [ g iiM ( X )] k . Proof of (30) . We show that | F ∞ (¯ α, ¯ ω, ¯ c ) | k ≤ k s g i ∞ ( X ) for M + 1 ≤ k . Since π ′ M (¯ x ) = ¯ x and [¯ c ] k = 0 for k ≥ M + 1, it follows that:[ F ∞ (¯ α, ¯ ω, ¯ c )] k = ( k ≤ M P k − j =1 e − iωj ¯ c j ¯ c k − j otherwise. (36)As ¯ c j ¯ c k − j = 0 when either j > M or k − j > M , then it follows that: | F ∞ (¯ α, ¯ ω, ¯ c ) | k ≤ M X j = k − M | ¯ c j ¯ c k − j | . Noting that | F ∞ (¯ α, ¯ ω, ¯ c ) | k = 0 for k > M , we calculate: | F ∞ (¯ α, ¯ ω, ¯ c ) | k ≤ k − s max M +1 ≤ k ≤ M k s M X j = k − M | ¯ c j ¯ c k − j | = k − s g i ∞ ( X ) . Proof of (31) . We show that | h ∗ c | k ≤ k s g ii,a ∞ ( X ) for M + 1 ≤ k and all c ∈ π c ( X ) and h ∈ π c ( H ). Fix x = ( α, ω, c ) ∈ X and h ∈ π c ( H ), and write c M = π M ( c ) , c ∞ = π ∞ ( c ) , h M = π M ( h ), and h ∞ = π ∞ ( h ). We may expand h ∗ c as follows: h ∗ c = h M ∗ c M + h M ∗ c ∞ + c M ∗ h ∞ + h ∞ ∗ c ∞ . (37)14he composition h M ∗ c M only has non-zero components for M + 1 ≤ k ≤ M , thereby itis bounded by the computable value below: h M ∗ k c M ≤ k s max { k s · h M ∗ k c M : M + 1 ≤ k ≤ M }≤ k s max M +1 ≤ k ≤ M k s M X j = k − m | H | j | X | k − j . (38)We calculate c M ∗ h ∞ for k ≥ M + 1, noting that [ h ∞ ] k − j = 0 if k − j ≤ M , as below: c M ∗ k h ∞ = k − X j =1 [ c M ] j [ h ∞ ] k − j + ∞ X j =1 [ c ∗ M ] j [ h ∞ ] k + j + [ c M ] k + j [ h ∗∞ ] j = M X j = k − M − [ c M ] j [ h ∞ ] k − j + M X j =1 [ c ∗ M ] j [ h ∞ ] k + j Using the estimates | c j | ≤ | X | j for j ≤ M and | h j | ≤ C /j s for j ≥ M + 1, we calculate thefollowing: | c M ∗ h ∞ | k ≤ M X j = k − M − | X | j C ( k − j ) s + M X j =1 | X | j C ( k + j ) s ≤ C k s M X j = k − M − | X | j (cid:18) kk − j (cid:19) s + M X j =1 | X | j . (39)Note that kk − j is decreasing with k . To maximize the coefficient of | X | j in the first sum of(39), we choose the smallest k such that j ≤ k − M −
1. Hence, for each coefficient, wechoose k = M + j + 1 as an upper bound. We obtain the following: | c M ∗ h ∞ | k ≤ C k s M X j =1 | X | j (cid:18)(cid:18) M + j + 1 M + 1 (cid:19) s + 1 (cid:19) . (40)An analogous calculation produces a bound for | h M ∗ c ∞ | as given below: | h M ∗ c ∞ | k ≤ C k s M X j =1 | H | j (cid:18)(cid:18) M + j + 1 M + 1 (cid:19) s + 1 (cid:19) . (41)Lastly we estimate | h ∞ ∗ c ∞ | k . For h ∞ , c ∞ ∈ H ∞ and k ≥ M + 1 we calculate: h ∞ ∗ c ∞ = k − X j =1 [ h ∞ ] j [ c ∞ ] k − j + ∞ X j =1 [ h ∗∞ ] j [ c ∞ ] k + j + [ h ∞ ] k + j [ c ∗∞ ] j = k − M − X j = M +1 [ h ∞ ] j [ c ∞ ] k − j + ∞ X j = M +1 [ h ∗∞ ] j [ c ∞ ] k + j + [ h ∞ ] k + j [ c ∗∞ ] j . | h j | ≤ C /j s for M + 1 ≤ j we obtain: | h ∞ ∗ c ∞ | k ≤ k − M − X j = M +1 C j s ( k − j ) s + 2 ∞ X j = M +1 C j s ( k + j ) s ≤ C k − M − X j = M +1 j s ( k − j ) s + 2 k s C ( s − M s − . The remaining sum is only nonzero for k ≥ M + 1), and we bound it as follows: k − M − X j = M +1 j s ( k − j ) s = 1 k s k − M − X j = M +1 (cid:18) j + 1 k − j (cid:19) s ≤ k s k/ X j = M +1 (cid:18) j (cid:19) s ≤ s +1 k s ( s − (cid:18) M s − − k/ s − (cid:19) . This estimate is maximized in the k · k s norm by taking k → ∞ . Thereby, we obtain thefollowing estimate: | h ∞ ∗ c ∞ | k ≤ k s C (2 s + 1)( s − M s − . (42)By combining the results from (38 - 42) into (37), it follows that if M + 1 ≤ k , then | h ∗ c | k ≤ k s g ii,a ∞ ( X ). Proof of (32) . We show that (cid:12)(cid:12) ( K − π c ( x )) ∗ π c ( x ) (cid:12)(cid:12) k ≤ k s − g ii,b ∞ ( X ) for M + 1 ≤ k and all x , x ∈ X . For i = 1 , c i ∈ π c ( X ) and recall that C ≥ k c i k s by Definition 2.11.We can write ( K − c ) ∗ k c as below:( K − c ) ∗ k c = k − X j =1 j [ c ] j [ c ] k − j + ∞ X j = k +1 j [ c ∗ ] j [ c ] k + j + ( k + j )[ c ] k + j [ c ∗ ] j . (43)Using | c | j ≤ C /j s and | c | k + j ≤ C / ( k + j ) s for k ≥ M + 1, we obtain a bound on | ( K − c ) ∗ c | k as below: | ( K − c ) ∗ c | k ≤ k − X j =1 jC C j s ( k − j ) s + ∞ X j =1 C C j s − ( k + j ) s + ∞ X j =1 C C ( k + j ) s − j s ≤ C k − X j =1 j s − ( k − j ) s + C C ( k + 1) s (cid:18) s − (cid:19) + C C ( k + 1) s − (cid:18) s − (cid:19) . Since 5 ≤ M , thereby 6 ≤ M + 1 ≤ k and by Lemma 2.10 we can simplify the remainingsum as follows: k − X j =1 j s − ( k − j ) s = k k − X j =1 j s ( k − j ) s ≤ k γ k k s ≤ γ M +1 k s − . k ≥ M + 1, it follows that: | ( K − π c ( x )) ∗ π c ( x ) | k ≤ k s − (cid:18) C γ M +1 C C (cid:18) s − M + 2)( s −
2) + ss − (cid:19)(cid:19) = 1 k s − g ii,b ∞ ( X ) . When defining a Krawczyk operator K ( X, ¯ x ) for a function f : Y → Z one must choosea linear operator A † : Z → Y . The map A † is typically chosen to approximate Df (¯ x ) − .Even in finite dimensions it may be impossible to exactly calculate the inverse of a matrixusing floating point arithmetic. To denote a fixed but numerically approximate definition,we introduce the notation : ≈ . Since we set up our theorems in an a posteriori format,the question of whether our numerical approximation is sufficiently accurate is answered bywhether our computer-assisted proof is successful or not.As with any method relying on a contraction mapping argument, the Krawczyk operatoris only truly effective in locating the zeros of a function if they are isolated. Since the non-trivial zeros of F are not isolated, and in fact form a 2-manifold [23], we do not define aKrawczyk operator corresponding directly to F : R × Ω s → Ω s − . We must first reducethe dimensionality of its domain by two.We reduce one of the dimensions by imposing a phase condition; we may assume withoutloss of generality that the first Fourier coefficient is a positive real number (see Proposition5.4). To that end, we define a codimension − s ⊆ Ω s as follows:˜Ω s := { c ∈ Ω s : c = c ∗ } . To reduce the other dimension, we consider α as a parameter and perform our estimatesuniformly in α .For a cube X ⊆ R × ˜Ω s we define a Krawczyk operator to find the zeros of functions F α : R × ˜Ω s → Ω s − for all α ∈ π α ( X ) . To that end, we would like to define a map A † to be an approximate inverse of the derivative DF ¯ α (¯ ω, ¯ c ) ∈ L ( R × ˜Ω s , Ω s − ) for some(¯ α, ¯ ω, ¯ c ) ∈ X . We construct this approximate inverse by combining A † M , a 2 M × M realmatrix on the lower Fourier modes, with the operator − ( i ¯ α ¯ ω ) K π ′∞ on the higher Fouriermodes.As is ever the case, we may only explicitly perform a finite number of operations onfundamentally finite dimensional objects, and because of this we defined Galerkin projectionsin (14) and (15). To ensure the sum F = F M + F ∞ makes sense, the maps π M , π ′ M aredefined to be but finite rank maps onto a subspace of an infinite dimensional Banach space.To emphasize this finite dimensional subspace as a space in its own right, as well as the newdomain R × ˜Ω s , we define the following projection and inclusion maps:˜ π M : Ω s ։ R M , ˜ π ′ M : R × ˜Ω s ։ R M , ˜ i M : R M ֒ → Ω s , ˜ i ′ M : R M ֒ → R × ˜Ω s . ˜ π M ◦ ˜ i M = id R M , ˜ π ′ M ◦ ˜ i ′ M = id R M , ˜ i M ◦ ˜ π M = id Ω s , ˜ i ′ M ◦ ˜ π ′ M = id R × ˜Ω s . We define the linear operator A † below in Definition 3.1 as follows: We note that A † willbe injective if the 2 M × M matrix A † M has rank 2 M .17 efinition 3.1. Fix a cube X ⊆ R × ˜Ω s . For a point (¯ α, ¯ ω, ¯ c ) = ¯ x ∈ X such that ¯ x = π ′ M (¯ x ) , define the following linear operators: A M : ≈ ˜ π M ◦ DF ¯ α (¯ ω, ¯ c ) ◦ ˜ i ′ M A M ∈ L ( R M , R M ) A † M : ≈ A − M A † M ∈ L ( R M , R M ) A (¯ x, M ) := ˜ i M ◦ A M ◦ ˜ π ′ M + i ¯ ω ¯ α K − π ′∞ A (¯ x, M ) ∈ L ( R × ˜Ω s , Ω s − ) A † (¯ x, M ) := ˜ i ′ M ◦ A † M ◦ ˜ π M − i ¯ α ¯ ω K π ∞ A † (¯ x, M ) ∈ L (Ω s − , R × ˜Ω s ) . While a Krawczyk operator K ( X, ¯ x ) given as in Definition 2.1 is sufficient from a math-ematical perspective, from a computational perspective it leaves something to be desired.We address this deficiency in Definition 3.2 by defining an explicitly computable operator K ′ ( X, ¯ x ) as an outer approximation to K ( X, ¯ x ), which is to say that K ( X, ¯ x ) ⊆ K ′ ( X, ¯ x ).In Theorem 3.3 we prove this, and in Theorem 3.4 we give an analogue of Theorem 2.2.In practice, use interval arithmetic [18] to compute an outer approximations for thearithmetic combination of sets (e.g. A + B = S a ∈ A,b ∈ B a + b ). This allows us to boundthe image of functions over rectangular domains, which is to say domains given as theproduct of intervals. By employing outward rounding, interval arithmetic can be rigorouslyimplemented on a computer [24]. In every step an outer approximation is constructedas a rectangular domain, and the end result will too be an outer approximation. Whileobtaining a tight approximation is desirable, it is not required; as long as we have an outerapproximation, that is sufficient. Definition 3.2.
Fix a cube X ⊆ R × ˜Ω s as in Definition 2.7 with M ≥ , s > and C > . Fix some ¯ x = (¯ α, ¯ ω, ¯ c ) ∈ X such that ¯ x = π ′ M (¯ x ) and ∆ ω ≥ sup x ∈ X | π ω ( x ) − ¯ ω | .Fix A := A (¯ x, M ) and A † := A † (¯ x, M ) as in Definition 3.1. Define the following functions: g ii ∞ ( X ) := 2 ¯ α ¯ ω ( M + 1) g ii,a ∞ ( X ) + sup α ∈ π α ( X ) ∆ ω ¯ α ¯ ω (cid:0) ( α − + 1) C + g ii,b ∞ ( X ) (cid:1) + sup α ∈ π α ( X ) ,ω ∈ π ω ( X ) (cid:18) | − ¯ αα ω ¯ ω | + ¯ α ¯ ω ( M + 1) (cid:19) C (44) g M ( X ) := g iM ( X ) + g iiM ( X ) (45) g ∞ ( X ) := ¯ α/ ¯ ωM +1 g i ∞ ( X ) + g ii ∞ ( X ) . (46) Define K ′ ( X, ¯ x ) := K ′ M ( X, ¯ x ) × K ′∞ ( X, ¯ x ) by: K ′ M ( X, ¯ x ) := ¯ x − A † M F M (¯ x ) + ( I M − A † M A M ) · π ′ M ( X − ¯ x )+ A † M ( A M − DF M ( X ))( X − ¯ x ) ± A † M g M ( X ) (47) K ′∞ ( X, ¯ x ) := { c k ∈ C : | c k | < g ∞ ( X ) /k s } ∞ k = M +1 , (48) where F M (¯ x ) ⊆ R M is calculated to include the image of F M (¯ x ) for all α ∈ π α ( X ) , where DF M ( X ) ⊆ L ( R M , R M ) is calculated to include the image of ˜ π M ◦ DF α ( ω, c ) ◦ ˜ i ′ M for all ( α, ω, c ) ∈ X , and where ± A † M g M ( X ) ⊆ R M is calculated to be a set satisfying: [ | v | k ≤| g M ( X ) | k A † M · v ⊆ ± A † M g M ( X ) . heorem 3.3. Fix a cube X as in Definition 2.7 with M ≥ , s > and C > . Fix apoint ¯ x ∈ X such that ¯ x = π ′ M (¯ x ) , and fix A := A (¯ x, M ) , A † := A † (¯ x, M ) as in Definition3.1. Fix some α ∈ π α ( X ) , and for f ≡ F α : R × ˜Ω s → Ω s − let K be given as in Definition2.1. Then K ( X, ¯ x ) ⊆ K ′ ( X, ¯ x ) . Proof.
Let H := X − ¯ x . We begin by proving that π ′ M ( K ( X, ¯ x )) ⊆ π ′ M ( K ′ ( X, ¯ x )), firstshowing that: π ′ M ◦ ( I − A † DF ( X )) · H ⊆ K ′ M ( X, ¯ x ) − (cid:16) ¯ x − A † M F M (¯ x ) (cid:17) . (49)Fix some x ∈ X and h = ( h ω , h c ) ∈ H . We start by adding and subtracting A † A , rewritingthe LHS of (49) as follows: π ′ M ( I − A † DF ( x )) · h =( I M − A † M A M ) · π ′ M ( h ) + π ′ M A † ( A − DF ( x )) · h =( I M − A † M A M ) · π ′ M ( h )+ A † M ( A M − DF M ( x )) · π ′ M ( h ) + A † M π M DF ∞ ( x ) · π ′ M ( h ) . By (28) and (29) it follows that | π M DF ∞ ( x ) · h | k ≤ [ g iM ( X ) + g iiM ( X )] k . Thereby, it followsthat: A † M π M DF ∞ ( x ) · h ⊆ ±| A † M | · g M ( X ) for all x ∈ X and h ∈ H . Hence from the defi-nition of K ′ ( X, ¯ x ) given in (47), then (49) follows. From (36) we have that π M F ∞ (¯ x ) = 0,hence π ′ M (¯ x − A † F (¯ x )) = ¯ x − A † M F M (¯ x ). It then follows that π M ◦ K ( X, ¯ x ) ⊆ K ′ M ( X, ¯ x ).We now prove that π ′∞ ( K ( X, ¯ x )) ⊆ π ′∞ ( K ′ ( X, ¯ x )), first showing that: (cid:13)(cid:13) π ′∞ ◦ ( I − A † DF ( X )) · ( X − ¯ x ) (cid:13)(cid:13) s ≤ g ii ∞ ( X ) . (50)Fix some x = ( α, ω, c ) ∈ X and h = ( h ω , h c ) ∈ H . We start by adding and subtracting A † A , rewriting the LHS of (50) as follows: π ′∞ ( I − A † DF ( x )) · h = π ′∞ ( I − A † A ) · h + π ′∞ A † ( A − DF ( x )) · h = π ′∞ ◦ A † ( A − DF ( x )) · h = π ′∞ ◦ A † (cid:0) A − ∂∂c DF ( x ) (cid:1) · h c − π ′∞ ◦ A † ∂∂ω DF ( x ) · h ω . We calculate − π ∞ A † ∂∂ω F ( x ) · h ω writing ∂∂ω F ( x ) as in (11) below: − π ∞ ◦ A † ∂∂ω F ( X )) · h ω = − iπ ∞ ¯ α ¯ ω K (cid:0) i K − ( α − I − U ω ) c − i ( K − U ω c ) ∗ c (cid:1) · h ω = h ω ¯ α ¯ ω π ∞ (cid:0) ( α − I − U ω ) c − K ( K − U ω c ) ∗ c (cid:1) . Using | c | j ≤ C /j s and (32) we obtain for k ≥ M + 1 that: (cid:12)(cid:12) π ∞ ◦ A † ∂∂ω F ( x )) · ∆ ω (cid:12)(cid:12) k ≤ ∆ ω ¯ α ¯ ω (cid:18) ( α − + 1) C k s + 1 k g ii,b ∞ ( X ) k s − (cid:19)(cid:13)(cid:13) π ∞ ◦ A † ∂∂ω F ( x )) · ∆ ω (cid:13)(cid:13) s ≤ ∆ ω ¯ α ¯ ω (cid:0) ( α − + 1) C + g ii,b ∞ ( X ) (cid:1) . (51)For ( α, ω, c ) ∈ X we calculate π ∞ A † ( A − ∂∂c F ) · h c below: π ∞ ◦ A † ( A − ∂∂c F ( x )) · h c = − i ¯ α ¯ ω K (cid:0)(cid:0) i ¯ ω ¯ α K − − ( i ωα K − + U ω ) (cid:1) h c − ( U ω h c ) ∗ c − ( U ω c ) ∗ h c (cid:1) = π ∞ (cid:0) (1 − ¯ αα ω ¯ ω ) I + i ¯ α ¯ ω K U ω (cid:1) h c − π ∞ i ¯ α ¯ ω K (( U ω c ) ∗ h c + ( U ω h c ) ∗ c ) . (cid:13)(cid:13) π ∞ ◦ A † ( A − ∂∂c F ( x )) · h c (cid:13)(cid:13) s ≤ (cid:18) | − ¯ αα ω ¯ ω | + ¯ α ¯ ω ( M + 1) (cid:19) C + 2 ¯ α ¯ ω ( M + 1) g ii,a ∞ ( X ) . (52)By combining (51) and (52) and taking a supremum over α and ω , we obtain the definitionof g ii ∞ in (44), whereby (50) follows.To show that π ∞ K ( X, ¯ x ) ⊆ K ′∞ ( X, ¯ x ) note that from (30) it follows that: k π ∞ (¯ x − A † F (¯ x )) k s = k − i ¯ α ¯ ω K π ∞ F (¯ x ) k s ≤ ¯ α/ ¯ ωM + 1 g i ∞ ( X ) . Expanding out π ∞ K ( X, ¯ x ), it follows that: k π ∞ K ( X, ¯ x ) k s ≤k π ∞ (¯ x − A † F (¯ x )) k s + k π ∞ ( I − ADF ( X )) · ( X − ¯ x ) k s ≤ ¯ α/ ¯ ωM + 1 g i ∞ ( X ) + g ii ∞ ( X ) = g ∞ ( X ) . Thus π ∞ K ( X, ¯ x ) ⊆ K ′∞ ( X, ¯ x ). Thus, we have proved both that π ′ M ( K ′ ( X, ¯ x )) ⊆ π ′ M ( K ( X, ¯ x ))and π ′∞ ( K ′ ( X, ¯ x )) ⊆ π ′∞ ( K ( X, ¯ x )). Hence it follows that K ( X, ¯ x ) ⊆ K ′ ( X, ¯ x ). Theorem 3.4.
Fix a cube X as in Definition 2.7 with M ≥ , s > and C > . Fix apoint ¯ x ∈ X such that ¯ x = π ′ M (¯ x ) . Let K ( X, ¯ x ) and K ′ ( X, ¯ x ) be given as in Definition 2.1and 3.2 respectively. If K ′ ( X, ¯ x ) ⊆ X , and moreover g ∞ ( X ) < C and: ˜ π ′ M (cid:16) K ′ M ( X, ¯ x ) + A † M F M (¯ x ) (cid:17) ⊆ int (˜ π ′ M ( X )) , then for all α ∈ π α ( X ) there exists a unique point ˆ x α = ( α, ˆ ω α , ˆ c α ) ∈ X such that F (ˆ x α ) = 0 .Proof. Fix α ∈ π α ( X ). By Theorem 2.2, in order to show that there exists a unique solutionto F α = 0, it suffices to show that there is some 0 ≤ λ < I − A † DF ( X ))( X − ¯ x ) ⊆ λ ( X − ¯ x ) . We find a λ M which works for the π ′ M -projection and a λ ∞ which works for the π ′∞ -projection. Since K ( X, ¯ x ) ⊆ K ′ ( X, ¯ x ) by Theorem 3.3 and ˜ π ′ M (cid:16) K ′ M ( X, ¯ x ) + A † M F M (¯ x ) (cid:17) ⊆ int (˜ π ′ M ( X )), it follows from the definition of K ( X, ¯ x ) in (3) that:˜ π ′ M (cid:0) K ( X, ¯ x ) + A † F (¯ x ) (cid:1) ⊆ int (˜ π ′ M ( X ))˜ π ′ M (cid:0) ( I − A † DF ( X ))( X − ¯ x ) (cid:1) ⊆ int (˜ π ′ M ( X − ¯ x )) (53)Since ˜ π ′ M (cid:0) ( I − A † DF ( X ))( X − ¯ x ) (cid:1) is compactly contained inside of ˜ π ′ M ( X − ¯ x ) ⊆ R M ,there is some positive distance separating the LHS of (53) away from the boundary of˜ π ′ M ( X − ¯ x ). It follows that there must exist some 0 ≤ λ M < π ′ M (cid:0) ( I − A † DF ( X ))( X − ¯ x ) (cid:1) ⊆ λ M · ˜ π ′ M ( X − ¯ x ).Since K ′∞ ( X, ¯ x ) ⊆ π ′∞ X it follows that g ∞ ( X ) ≤ C , and by our additional assumptionthis is in fact a strict inequality. If we define λ ∞ := g ii ∞ ( X ) /C <
1, then by (50) it followsthat: π ∞ ( I − A † DF ( X )) · ( X − ¯ x ) ≤ λ ∞ π ∞ ( X − ¯ x ) .
20f we define λ := max { λ M , λ ∞ } < I − A † DF ( X )) · ( X − ¯ x ) ≤ λ ( X − ¯ x ) . By Theorem 2.2 there exists a unique point ˆ x α = ( α, ˆ ω α , ˆ c α ) ∈ X such that F α (ˆ ω α , ˆ c α ) = 0.Moreover, this is true for all α ∈ π α ( X ). For a given cube, we want to know if it contains any solutions to F = 0. We try to determinethis by combining several different tests into one pruning operator described in Algorithm4.1. It is called a pruning operator because even if we cannot determine whether a cubecontains a solution, we may still be able to reduce the size of the cube without losing anysolutions.We describe the tests performed in Algorithm 4.1. Most simply, if we can prove that | F ( X ) | k > ≤ k ≤ M , then F has no zeros in X . From Lemma 2.8, we know thatif a cube has a small k · k ℓ norm then it cannot contain any nontrivial zeros. Furthermore, ifa cube is contained in the neighborhood of the Hopf bifurcation explicitly given by Lemma2.9, then the only solutions that can exist therein are on the principal branch. If none ofthose situations apply, then we calculate the outer approximation of the Krawczyk operatorgiven in Definition 3.2. If the hypothesis of Theorem 3.4 is satisfied, then there exists aunique solution. Alternatively, if X ∩ K ( X, ¯ x ) = ∅ , then there do not exist any solutions in X . If none of these other situations apply, then we replace X by X ∩ K ( X, ¯ x ). Algorithm4.1 arranges these steps in order of ease of computation. Algorithm 4.1 (Prune) . Take as input a cube X with M ≥ and s > . The output is apair { f lag, X ′ } where f lag ∈ Z and X ′ ⊆ X is a cube.1. Compute δ := 2 P Mk =1 | X | k + C ( s − M s − .2. If for all ( α, ω, · ) ∈ X we have α ∈ (0 , , ω ≥ . , and δ < g ( α, ω ) for g defined in (20) , then return { , ∅} .3. If for all ( α, ω, · ) ∈ X we have | α − π | ≤ . , | ω − π | ≤ . and δ < . , thenreturn { , X } .4. If inf x ∈ X | F M ( x ) | k > h k ( X ) for h k defined in (21) and some ≤ k ≤ M , then return { , ∅} .5. Fix some ¯ x ∈ X such that ¯ x = π ′ M (¯ x ) and π ′ M (¯ x ) is approximately the center of π ′ M ( X ) . Construct K ′ ( X, ¯ x ) as in Definition 3.2.6. If K ′ ( X, ¯ x ) ⊆ X , g ∞ ( X ) < C , and ˜ π M (cid:16) K ′ M ( X, ¯ x ) + A † M F M (¯ x ) (cid:17) ⊆ int (˜ π M ( X )) ,then return { , X } .7. If X ∩ K ′ ( X, ¯ x ) = ∅ , then return { , ∅} .8. Else return { , X ∩ K ′ ( X, ¯ x ) } . Theorem 4.2.
Let { f lag, X ′ } denote the output of Algorithm 4.1 with input a cube X . i) If f lag = 1 , then F ( x ) = 0 for all nontrivial x ∈ X .(ii) If f lag = 2 , then the only solutions to F = 0 in X are on the principal branch.(iii) If f lag = 3 , then for all α ∈ π α ( X ) there is a unique ˆ ω α ∈ π ω ( X ) and ˆ c α ∈ π c ( X ) such that F ( α, ˆ ω α , ˆ c α ) = 0 .(iv) If there are any points ˆ x ∈ X for which F (ˆ x ) = 0 , then ˆ x ∈ X ′ .Proof. To prove ( i ) we must check the output from Steps 2, 4, and 7. To prove ( ii ) we mustcheck Step 3. To prove ( iii ) we must check Step 6. The proof of ( iv ) follows from ( i ), ( ii ),( iii ), and Step 8. We organize the proof into the steps of the algorithm.1. It follows from (19) that k c k ℓ < δ for all c ∈ π c ( X ).2. Since α ∈ (0 ,
2] and ω ≥ .
1, Lemma 2.8 applies. If k c k ℓ < δ < g ( α, ω ), then byLemma 2.8 the only solutions to F ( α, ω, c ) = 0 are trivial, which is to say c = 0.3. If Step 3 returns f lag = 2, then by Lemma 2.9 there is at most one SOPS c ∈ X withfrequency ω , and it lies on the branch of SOPS originating from the Hopf bifurcationat α = π .4. Suppose that inf x ∈ X | F M ( x ) | k > h k ( X ) for some 1 ≤ k ≤ M . Since sup x ∈ X | F ∞ ( x ) | k
0, and so X cannot contain any zeros of F .5. Note that K ( X, ¯ x ) ⊆ K ′ ( X, ¯ x ) by Theorem 3.4.6. If Step 6 returns f lag = 3, then the hypothesis of Theorem 3.4 is satisfied. Hence for all α ∈ π α ( X ) there is a unique ˆ ω α ∈ π ω ( X ) and ˆ c α ∈ π c ( X ) such that F ( α, ˆ ω α , ˆ c α ) = 0.7. By Theorem 2.2 all solutions in X are contained in K ( X, ¯ x ). Hence, all of the zerosof F in X are contained in X ∩ K ( X, ¯ x ) ⊆ X ∩ K ′ ( X, ¯ x ).If X ∩ K ′ ( X, ¯ x ) = ∅ then X ∩ K ( X, ¯ x ) = ∅ , whereby there cannot be any solutions in X .8. As proved in Step 7, all solutions in X are contained in X ∩ K ′ ( X, ¯ x ). The goal of this section is to construct a bounded region in R × Ω s which contains all ofthe nontrivial zeros of F . This is ultimately achieved in Algorithm 5.7, which is discussedin Section 5.2, along with other estimates pertaining specifically to Wright’s equation.In Section 5.1, we discuss generic algorithms used to construct bounds in Fourier space.Algorithm 5.1 converts pointwise bounds on a periodic function and its derivatives into acube containing its Fourier coefficients. Algorithm 5.3 modifies a cube so that after a timetranslation, any periodic function contained therein will satisfy the phase condition c = c ∗ .22 .1 Converting Pointwise Bounds into Fourier Bounds To translate pointwise bounds on a periodic function into bounds on its Fourier coefficientswe use the unnormalized L inner product, which we define for g, h ∈ L ([0 , π/ω ] , C ) as: h g, h i := Z π/ω g ( t ) h ( t ) ∗ dt. (54)For a function y given as in (4), its Fourier coefficients may be calculated as c k = π/ω (cid:10) y ( t ) , e iωkt (cid:11) .By applying (54) to a priori estimates on y we are able to derive bounds on its Fourier co-efficients. For example, in [28] it is shown that − < y ( t ) < e α − e α ≥ | c k | ≤ π/ω ( e α −
1) for all k ∈ Z .With more detailed estimates on y we can produce tighter bounds on its Fourier co-efficients. In [2, 10] such estimates are numerically derived in a rigorous fashion. One ofthe results from this analysis is a pair of bounding functions which provide upper and lowerbounds on SOPS to (2) at a given parameter value. Formally, a bounding function is definedto be an interval valued function χ ( t ) = [ ℓ ( t ) , u ( t )] where ℓ, u : R → R .These functions ℓ, u are constructed in [2, 10] using rigorous numerics, and in particularinterval arithmetic. As a matter of computational convenience, these functions are definedas piecewise constant functions which change value only finitely many times (see Figure3). For functions of this form, calculating a supremum over a bounded domain is reducedto finding the maximum of a finite set, and calculating an integral is reduced into a finitesum. For elementary functions such as sin or cos, interval arithmetic packages have beendeveloped which allow us to rigorously bound their image over arbitrary domains [24].23: Bounds for y k A k, B k, − . , . − . , − . − . , . − . , . − . , . − . , . y ′ k A k, B k, − . , . − . , − . − . , . − . , . − . , . − . , . y ′′ k A k, B k, − . , . − . , − . − . , . − . , . − . , . − . , . y ′′′ k A k, B k, − . , . − . , . − . , . − . , . − . , . − . , . Depicted in the figures are functions ℓ s , u s : R → R which bound a periodic function y and itsderivatives y ( s ) . Depicted in the tables are the values for A k,s and B k,s produced by Algorithm 5.1 whichbound the Fourier coefficients c k = a k + ib k of y . Algorithm 5.1 describes a method for obtaining rigorous bounds on the Fourier coeffi-cients of a periodic function y . This algorithm applies the inner product h· , ·i to bounds notjust on the function y but on its derivatives as well. Examples of these bounds are given inFigure 3, where we note that by the third Fourier coefficient, the tightest estimate is givenby the third derivative. We will use y ( s ) denotes the s th derivative of a function y , whereaswe will use Y s to denote a bounding function of index s , which bounds the derivative y ( s ) .We have stated Algorithm 5.1 so that it does not estimate the zeroth Fourier coefficient,as periodic solutions to (2) necessarily have a trivial zeroth Fourier coefficient. The algorithmcould be modified in the obvious way to bound the zeroth Fourier coefficient of a functionas well. Algorithm 5.1.
Take as input projection dimension M ∈ N , period bounds [ L, L ] , and acollection of interval-valued functions: { Y s ( t ) = [ ℓ s ( t ) , u s ( t )] : ℓ s , u s : R → R } Ss =0 . The output is an ( α -parameterless) cube X ⊆ R × Ω S .1. Define I ω := [2 π/L, π/L ] . . For ≤ k ≤ M and ≤ s ≤ S define δ c , δ s ∈ R + so that: δ c ≥ sup ω ∈ I ω ,y s ∈ Y s Z LL | cos( ωkt ) y s ( t ) | dt, δ s ≥ sup ω ∈ I ω ,y s ∈ Y s Z LL | sin( ωkt ) y s ( t ) | dt, and define a + k,s , a − k,s , b + k,s , b − k,s ∈ R + so that: a + k,s ≥ δ c + sup ω ∈ I ω ,y s ∈ Y s Z L cos( ωkt ) y s ( t ) dta − k,s ≤ − δ c + inf ω ∈ I ω ,y s ∈ Y s Z L cos( ωkt ) y s ( t ) dtb + k,s ≥ δ s + sup ω ∈ I ω ,y s ∈ Y s Z L sin( ωkt ) y s ( t ) dtb − k,s ≤ − δ s + inf ω ∈ I ω ,y s ∈ Y s Z L sin( ωkt ) y s ( t ) dt.
3. For ≤ k ≤ M and ≤ s ≤ S define: A ′ k,s := 12 πk s " inf ω ∈ I ω a − k,s ω s − , sup ω ∈ I ω a + k,s ω s − , B ′ k,s := 12 πk s " inf ω ∈ I ω b − k,s ω s − , sup ω ∈ I ω b + k,s ω s − . (55) Define the intervals A k,s and B k,s as follows: A k,s := A ′ k,s if s ≡ − B ′ k,s if s ≡ − A ′ k,s if s ≡ B ′ k,s if s ≡ , B k,s := − B ′ k,s if s ≡ − A ′ k,s if s ≡ B ′ k,s if s ≡ A ′ k,s if s ≡ .
4. For ≤ k ≤ M define: A k := \ ≤ s ≤ S A k,s , B k := \ ≤ s ≤ S B k,s .
5. For each ≤ k ≤ M , define ¯ a k := mid ( A k,S ) , ¯ b k := mid ( B k,S ) , ¯ c k = ¯ a k + i ¯ b k , and ¯ c − k = ¯ c ∗ k . Define y SM ( t, ω ) as in (56) , and define C > so that (57) holds. y SM ( t, ω ) := M X k = − M ¯ c k ( iωk ) S e iωkt (56) C ≥ sup ω ∈ I ω ,y S ∈ Y S πω S − Z L (cid:12)(cid:12) y S ( t ) − y SM ( t, ω ) (cid:12)(cid:12) dt. (57)
6. Define a cube X := X M × X ∞ ⊆ R × Ω S by: X M := I ω × M Y k =1 A k × B k X ∞ := (cid:8) c k ∈ C : | c k | ≤ C /k S (cid:9) ∞ k = M +1 . roposition 5.2. Let the cube X be the output of Algorithm 5.1 with input M ∈ N , [ L, L ] ⊆ R and bounding functions { Y s } Ss =0 . Fix a function ˆ y with period L and con-tinuous derivatives ˆ y ( s ) for ≤ s ≤ S . If L ∈ [ L, L ] and ˆ y ( s ) ( t ) ∈ Y s ( t ) for all ≤ s ≤ S and t ∈ [0 , L ] , then the frequency and Fourier coefficients of ˆ y satisfy ( ω, { c k } ∞ k =1 ) ∈ X .Proof. We organize the proof into the steps of the algorithm.1. If the period of ˆ y is L ∈ [ L, L ] then it will have frequency ˆ ω = 2 π/L and ˆ ω ∈ [2 π/L, π/L ].2. Let us define a k,s := D cos(ˆ ωkt ) , ˆ y ( s ) ( t ) E , b k,s := D sin(ˆ ωkt ) , ˆ y ( s ) ( t ) E . We show that a k,s ∈ [ a − k,s , a + k,s ]. Since L ∈ [ L, L ] it follows that: D cos(ˆ ωkt ) , ˆ y ( s ) ( t ) E = Z L cos(ˆ ωkt )ˆ y ( s ) ( t ) dt = Z L cos(ˆ ωkt )ˆ y ( s ) ( t ) dt + Z LL cos(ˆ ωkt )ˆ y ( s ) ( t ) dt. (58)To estimate the rightmost summand in (58) we calculate: (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Z LL cos(ˆ ωkt )ˆ y ( s ) ( t ) dt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Z LL (cid:12)(cid:12)(cid:12) cos(ˆ ωkt )ˆ y ( s ) ( t ) (cid:12)(cid:12)(cid:12) dt ≤ sup ω ∈ I ω ,y s ∈ Y s Z LL | cos( ωkt ) y s ( t ) | dt ≤ δ c . We obtain a bound on a k,s by appropriately taking an infimum and a supremum in(58) as follows:inf ω ∈ I ω ,y s ∈ Y s Z L cos( ωkt ) y s ( t ) dt − δ c ≤ a k,s ≤ sup ω ∈ I ω ,y s ∈ Y s Z L cos( ωkt ) y s ( t ) dt + δ c . Hence a k,s ∈ [ a − k,s , a + k,s ], and by analogy b k,s ∈ [ b − k,s , b + k,s ].3. Let c k = a k + ib k denote the Fourier coefficients of ˆ y . We show that a k ∈ A k,s and b k ∈ B k,s . Firstly, we calculate the derivative ˆ y ( s ) as follows:ˆ y ( s ) ( t ) = X k ∈ Z c k ( i ˆ ωk ) s e i ˆ ωkt . We can express the Fourier coefficients of ˆ y in terms of the Fourier coefficients of itsderivatives ˆ y ( s ) ; below, we calculate c k in terms of a k,s and b k,s as follows: Z π/ ˆ ω c k ( i ˆ ωk ) s e i ˆ ωkt · e − i ˆ ωkt dt = D ˆ y ( s ) ( t ) , e i ˆ ωkt E (59)2 π ˆ ω c k ( i ˆ ωk ) s = D ˆ y ( s ) ( t ) , cos(ˆ ωkt ) E − i D ˆ y ( s ) ( t ) , sin(ˆ ωkt ) E i s a k + i s +1 b k = a k,s − i b k,s π ˆ ω s − k s . A ′ k,s and B ′ k,s in (55) it follows that: a k,s π ˆ ω s − k s ∈ A ′ k,s , b k,s π ˆ ω s − k s ∈ B ′ k,s . By matching the real and imaginary parts, which only depend on s (mod 4), we obtainthat a k ∈ A k,s and b k ∈ B k,s .4. Since a k ∈ A k,s and b k ∈ B k,s for all k and 0 ≤ s ≤ S , it follows that: a k ∈ \ ≤ s ≤ S A k,s , b k ∈ \ ≤ s ≤ S B k,s .
5. We calculate c k for k ≥ M + 1 starting from (59) and using the fact that the functions e i ˆ ωkt are L –orthogonal: c k ( i ˆ ωk ) S = 12 π/ ˆ ω D e i ˆ ωkt , ˆ y ( S ) ( t ) E = 12 π/ ˆ ω * e i ˆ ωkt , ˆ y ( S ) ( t ) − M X j = − M ¯ c j ( i ˆ ωj ) S e i ˆ ωjt + = 12 π/ ˆ ω D e i ˆ ωkt , ˆ y ( S ) ( t ) − y SM ( t, ˆ ω ) E . By taking absolute values, and the suprema over ω ∈ I ω and y S ∈ Y S we obtain thefollowing. (cid:12)(cid:12) c k ( i ˆ ωk ) S (cid:12)(cid:12) ≤ π/ ˆ ω Z L (cid:12)(cid:12) e − i ˆ ωkt (cid:12)(cid:12) (cid:12)(cid:12)(cid:12) ˆ y ( S ) ( t ) − y SM ( t, ˆ ω ) (cid:12)(cid:12)(cid:12) dt | c k | k S ≤ sup ω ∈ I ω ,y S ∈ Y S πω S − Z L (cid:12)(cid:12) y S ( t ) − y SM ( t, ω ) (cid:12)(cid:12) dt ≤ C . Hence | c k | ≤ C /k S for all k ≥ M + 1.6. In Step 1 we showed that ˆ ω ∈ I ω . In Steps 2-4 we showed that c k ∈ [ X ] k for 1 ≤ k ≤ M ,and in Step 5 we showed that | c k | ≤ C /k S for k ≥ M + 1. Algorithm 5.3.
Take as input an ( α -parameterless) cube X ⊆ R × Ω s . The output is an( α -parameterless) cube X ′ ⊆ R × ˜Ω s .
1. For [ X ] = A × B , with A = [ A , A ] and B = [ B , B ], define an interval Θ ⊆ R so that: Θ ⊇ S a ∈ A ,b ∈ B tan − ( b /a ) if A > S a ∈ A ,b ∈ B tan − ( b /a ) + π if A < S a ∈ A ,b ∈ B − tan − ( a /b ) + π if B > S a ∈ A ,b ∈ B − tan − ( a /b ) − π if B < − π, π ] otherwise.27. Rotate every Fourier coefficient’s phase by − Θ k . That is, define: A ′ := " inf a ∈ A ,b ∈ B q a + b , sup a ∈ A ,b ∈ B q a + b , B ′ := [0 , , and for 2 ≤ k ≤ M define intervals A ′ k , B ′ k ⊆ R such that: A ′ k ⊇ [ θ ∈ Θ ,a k ∈ A k ,b k ∈ B k cos( θk ) a k + sin( θk ) b k B ′ k ⊇ [ θ ∈ Θ ,a k ∈ A k ,b k ∈ B k − sin( θk ) a k + cos( θk ) b k .
3. Define a cube X ′ := X ′ M × X ′∞ ⊆ R × Ω S by X ′ M := I ω × M Y k =1 A ′ k × B ′ k X ′∞ := (cid:8) c k ∈ C : | c k | ≤ C /k S (cid:9) ∞ k = M +1 . Proposition 5.4.
For an input cube X , let X ′ denote the output of Algorithm 5.3. Supposethat y : R → R is a periodic function given as in (4) with frequency and Fourier coefficientssatisfying ( ω, { c k } ∞ k =1 ) ∈ X . Then there exists some τ ∈ R such that the Fourier coefficients c ′ of y ( t + τ ) satisfy ( ω, { c ′ k } ∞ k =1 ) ∈ X ′ . Furthermore c ′ is a real non-negative number.Proof. We organize the proof into the steps of the algorithm.1. Write the first Fourier coefficient of y as c = a + ib . We may write c = re iθ where r = p a + b and if c = 0, then θ is unique up to an integer multiple of 2 π . By therules for arctan we can calculate: θ = tan − ( b /a ) if a > − ( b /a ) + π if a < − tan − ( a /b ) + π if b > − tan − ( a /b ) − π if b < . Since a ∈ A and b ∈ B , it follows that θ ∈ Θ.2. For any τ we can calculate the Fourier series of y ( t + τ ) as follows: y ( t + τ ) = X k ∈ Z c k e iωk ( t + τ ) = X k ∈ Z c k e iωkτ e iωkt . If we choose τ = − θ/ω , then c ′ = c e iωτ = p a + b is a real, non-negative numberand moreover c ′ ∈ [ X ′ ] .3. The Fourier coefficients of y ( t + τ ) are given by c ′ k = e − ikθ c k , hence ( ω, { c ′ k } ∞ k =1 ) ∈ X ′ .28 .2 Bounds for Wright’s Equation The culmination of this subsection is Algorithm 5.7 which, for a given range of parameters,constructs a collection of cubes covering the solution space to F α = 0. This algorithm beginswith pointwise bounds on SOPS to (2). To obtain these pointwise bounds, we use the resultsfrom [10]. One of the results [10] achieves is, for a given range of parameters I α , it producesa collection of bounding functions X , such that if there is a SOPS to the exponential versionof Wrights equation at parameter α ∈ I α , then it will be bounded by one of the boundingfunctions in X . Recall that solutions to the exponential version of Wrights equation solve(1) where f ( x ) = e x −
1, and can be transformed into the quadratic version of Wright’sequation (2) using the change of variable y = e x − a priori estimates, some of which are iteratively constructed, and require a selectionof parameters i , j ∈ N . These are used to construct numerical bounding functions havingtime resolution n T ime ∈ N . A pruning operator is defined on these bounding functions, andthe spacing between the zeros of a SOPS, and the parameter N P eriod ∈ N defines how manytimes this pruning operator is applied in this initial construction of the bounding functions.Then a branch and prune algorithm is executed, with a stopping criterion defined by theparameters ǫ , ǫ ∈ R . We formally state the results of this algorithm below: Theorem 5.5 (See Theorem 5.2 in [10]) . Fix some I α = [ α min , α max ] such that α min ≥ π .Suppose that x : R → R is periodic with period L , and is a SOPS to (1) at parameter α ∈ I α with f ( x ) = e x − . Furthermore, assume without loss of generality that x (0) = 0 and x ′ (0) > .If L and X denote the output of [10, Algorithm 5.1] ran with input I α , then there existssome [ L i , L i ] ∈ L and χ i ∈ X for which L ∈ [ L i , L i ] and x ( t ) ∈ χ i ( t ) for all t . In [10] the authors applied this algorithm to prove there is a unique SOPS for α ∈ [1 . , . α = π . To remedy this, wemodify the pruning operator in [10] with the addition of the following Proposition 5.6. Thisallows for a new way to potentially conclude that a given bounding function cannot containany SOPS. Proposition 5.6. If y is a nontrivial periodic solution to (2) at parameter α ∈ (0 , andfrequency ω ≥ . , then: sup | y ( t ) | > − + q √ ωπα g ( α, ω ) . Proof.
Define M := sup | y ( t ) | . From [26, Lemma 4.1] we know that if F ( α, ω, c ) = 0, then: k c k ℓ ≤ πω √ k y ′ k ∞ ≤ πω √ αM (1 + M ) . From Lemma 2.8, the only solutions satisfying k c k ℓ < g ( α, ω ) are trivial. Hence ( α, ω, c )would only be a trivial solution at best if the following inequality is satisfied: k c k ℓ ≤ πω √ αM (1 + M ) < g ( α, ω ) . Solving the quadratic equation M + M − ω √ πα g ( α, ω ) < y ′′ ( t ) = − α [ y ′ ( t −
1) [1 + y ( t )] + y ( t − y ′ ( t )] y ′′′ ( t ) = − α [ y ′′ ( t − y ( t )] + 2 y ′ ( t − y ′ ( t ) + y ( t − y ′′ ( t )] . Note that we can always express the derivative y ( s ) ( t ) in terms of y ( r ) ( t ) and y ( r ) ( t − ≤ r ≤ s −
1. That is, we can inductively define functions f s : R s → R such thatfor all t we have: y ( s ) ( t ) = f s (cid:16) y ( t ) , y ( t − , y ′ ( t ) , y ′ ( t − , . . . , y ( s − ( t ) , y ( s − ( t − (cid:17) . (60)If we start with a bounding function for y , then by appropriately adding and multiplyingthe bounding functions for y ( r ) , taking wider bounds whenever necessary, we can obtainbounding functions for any derivative of y (see for example Figure 3).Algorithm 5.7 proceeds by first constructing bounding functions for y and its derivatives,and then applying Algorithm 5.1 to obtain a cube containing its Fourier coefficients. Thenit applies Algorithm 5.3 to impose the phase condition that c = c ∗ . In this manner we a a a a Figure 4:
Depicted above is the output of Algorithm 5.7 projected onto the ( ω, a ) plane. From left toright, the input I α was taken to be [ π , . . , . . , . . , . C increases with α , a , and period length 2 π/ω . obtain a collection of cubes which contains all of the Fourier coefficients to SOPS to (2). Wethen apply Algorithm 4.1 to each cube, discarding it if possible. This allows us to discardbetween 5% and 60% of cubes (see N ′ grid in Table 1).One problem however, is that the Fourier projection of two distinct bounding functionsoften overlap considerably. To address this we combine overlapping cubes together. Whilewe could combine all of our cubes into one big cube, this would not be efficient. Instead, wedivide our cover along a grid in the ω × a plane (see Figure 4). Algorithm 5.7.
Fix an interval of I α ⊆ [ α min , α max ] , integers M, S ∈ N and a subdivisionnumber N ∈ N , and the computational parameters for [10, Algorithm 5.1]. The output is a(finite) collection of cubes S = { X i ⊆ R × ˜Ω s } .1. Let X , L be the output of [10, Algorithm 5.1] with input I α and appropriate computa-tional parameters. . Use the change of variables y = e x − to define a collection of functions: Y := n Y i ( t ) = [ e ℓ i ( t ) − , e u i ( t ) −
1] : χ i = [ ℓ i ( t ) , u i ( t )] ∈ X o .
3. Inductively define Y s for ≤ s ≤ S so that corresponding to each Y i ∈ Y there existsa Y si = [ Y si , Y si ] ∈ Y s such that for f s defined in (60) we have: Y si ( t ) ≤ inf { y r } s − r =0 ∈{ Y ri } s − r =0 f s (cid:0) y ( t ) , y ( t − , . . . , y s − ( t ) , y s − ( t − (cid:1) Y si ( t ) ≥ sup { y r } s − r =0 ∈{ Y ri } s − r =0 f s (cid:0) y ( t ) , y ( t − , . . . , y s − ( t ) , y s − ( t − (cid:1) .
4. Define S ′ := { X ′ i ⊆ R × Ω s } to be the collective output of Algorithm 5.1 run with M ∈ N , and each of the sets L i ∈ L and { Y si } Ss =0 ∈ {Y s } Ss =0 as input.5. Define S ′′ := { X ′′ i ⊆ R × ˜Ω s } to be the collective output of Algorithm 5.3 run witheach of the sets X ′ i ∈ S ′ as input.6. Define S ′′′ by taking the product of I α with the cubes in S ′′ . That is, define S ′′′ := { I α × X ′′ i ⊆ R × ˜Ω s : X ′′ i ∈ S ′′ } .7. For each X ∈ S ′′′ , let { f lag, X ′ } denote the output of Algorithm 4.1 with input X . If f lag = 1 , then remove X from S ′′′ . Otherwise replace X by X ′ .8. Subdivide the ω × a space covered by S ′′′ into an N × N grid. That is, define an indexset B := { , , . . . , N } × { , , . . . , N } and define intervals I ω , I a ⊆ R so that: I ω ⊇ [ X ∈S ′′′ π ω ( X ) , I a ⊇ [ X ∈S ′′′ π a ( X ) . Subdivide I ω and I a into N subintervals of equal width, { I ωi } Ni =1 and { I a i } Ni =1 , sothat I ω = S Ni =1 I ωi and I a = S Ni =1 I a i .9. For each β = ( β , β ) ∈ B , take the union of cubes in S ′′′ whose ( ω, a ) –projectionintersects I ωβ × I a β . That is, define: ˜ X β := { ( α, ω, c ) ∈ R × ˜Ω s : ω ∈ I ωβ , [ c ] ∈ I a β } , and define X β to be a cube such that: X β ⊇ [ X ∈S ′′′ X ∩ ˜ X β .
10. Define S := { X β : β ∈ B } . Theorem 5.8.
Fix an interval I α = [ α min , α max ] such that α min ≥ π , and let S denotethe output of Algorithm 5.7. If a function y as given in (4) is a SOPS to Wright’s equationat α ∈ I α , then there exists a time translation so that its Fourier coefficients are in S S .Proof. Every SOPS y to the quadratic version of Wright’s equation given in (2) correspondsto a SOPS x to the exponential version of Wright’s equation given in (1) with f ( x ) = e x − x : R → R to the exponential version of Wright’s equation with period L . Weorganize the proof into the steps of the algorithm.31. By Theorem 5.5 there exists an interval L i ∈ L and a bounding function χ i ∈ X andsuch that L ∈ L i and x ( t ) ∈ χ i ( t ) for all t ∈ R .2. The change of variables between the exponential and quadratic versions of Wright’sequation is given by y = e x −
1. Hence for the interval L i ∈ L and the boundingfunction Y i ∈ Y , it follows that L ∈ L i and y ( t ) ∈ Y i ( t ) for all t ∈ R .3. Since y ∈ Y i it follows that its derivatives satisfy y ( s ) ∈ Y si for all 0 ≤ s ≤ S .4. Let ω and c denote the frequency and Fourier coefficients of y respectively. If X ′ i isthe output of Algorithm 5.1 with input M ∈ N , L i and { Y si } Ss =0 , then by Proposition5.2 it follows that ( ω, { c k } ∞ k =1 ) ∈ X ′ i .5. Let X ′′ i denote the output of Algorithm 5.3 with input X ′ i . By Theorem 5.4, thereexists a τ ∈ R such that the Fourier coefficients c ′ of y ( t + τ ) satisfy ( ω, { c ′ k } ∞ k =1 ) ∈ X ′′ i .6. We have shown that if y is a SOPS to (2) at parameter α having frequency ω , thenup to a time translation ( α, ω, c ) ∈ S S ′′′ . By Proposition 2.4 the SOPS to (2) atparameter α ∈ I α correspond to the non-trivial zeros of F in S S ′′′ . Hence, if there isa solution F (ˆ x ) = 0 for some x ∈ R × ˜Ω s with π α (ˆ x ) ∈ I α , then ˆ x ∈ S S ′′′ .7. Let { f lag, X (4) i } denote the output of Algorithm 4.1 with input X ′′′ i ∈ S ′′′ . By The-orem 4.2 we can replace each X ′′′ ∈ S ′′′ with X (4) i , and it will still be the case that S S ′′′ contains all of the solutions to F = 0. In particular, if f lag = 1 then X (4) i = ∅ and we may remove X ′′′ i in this case.8. If ( α, ω, c ) ∈ S S ′′′ and a = [ c ] , then by construction ω ∈ I ω and a ∈ I a . As I ω × I a = S ( β ,β ) ∈ B I ωβ × I a β , then there is some ( β , β ) ∈ B such that ( ω, a ) ∈ I ωβ × I a β .9. As S X ∈S ′′′ X ⊆ S β ∈ B ˜ X β , then it follows that S X ∈S ′′′ X ⊆ S β ∈ B X β . That is to say S S ′′′ ⊆ S S .10. Hence, S S contains the Fourier coefficients of any possible SOPS. After Algorithm 5.7 has constructed a collection of cubes S covering the solution spaceto F = 0, we run a branch and prune algorithm. This algorithm iteratively inspects theelements in X ∈ S and then constructs three new lists of cubes: A , B and R . To summarize,first we compute the output P rune ( X ) = { f lag, X ′ } from Algorithm 4.1. If f lag = 1, thenthere are no solutions in X , and we can remove X from S . If f lag = 2, then the cube isin the neighborhood of the Hopf bifurcation, and we add X ′ to B . If f lag = 3, then for all α ∈ π α ( X ) there exists a unique solution to F α = 0 in X ′ , and we add X ′ to A . If X ′ is toosmall, then we add it to R . If the Krawczyk operator appears to be effective at reducing thesize of the cube, then the pruning operation is performed again. Otherwise X ′ is subdividedalong some lower dimension and the resulting pieces are added back to S .32he most obvious difference between our algorithm and the classical algorithm is thatwe are working in infinite dimensions. While we store 2 M + 1 real valued coordinates in agiven cube, as in [5, 7] the subdivision is only performed along a subset of these dimensions.Choosing which dimension to subdivide along can greatly affect the efficiency of a branch andbound algorithm, and there are heuristic methods for optimizing this choice [4]. Howeversince we are finding all the zeros along a 1-parameter family of solutions, these branchingmethods are not entirely applicable. To determine which dimension to subdivide we selectthe dimension with the largest weighted diameter. That is, for a collection of weights { λ i } di =0 we define: w ( X, i ) := ( λ i · diam ( π α ( X )) if i = 0 ,λ i · diam (cid:0) [˜ π ′ M ( X )] i (cid:1) otherwise. (61) Algorithm 6.1 (Branch & Bound) . Take as input a collection of cubes S = { X i ⊆ R × ˜Ω s } with M ≥ and s > , and as computational parameters: a halting criteria ǫ > , acontinue-pruning criteria δ ≥ , a maximum subdivision dimension ≤ d ≤ M and a setof weights { λ i } di =0 . The output is three lists of cubes: A , B and R .1. If S is empty, terminate the algorithm.2. Select an element X ∈ S and remove X from S .3. Define { f lag, X ′ } = P rune ( X ) to be the output of Algorithm 4.1 with input X .4. If f lag = 1 , then reject X and GOTO Step 1.5. If f lag = 2 , then add X ′ to B and GOTO Step 1.6. If f lag = 3 , then add X ′ to A and GOTO Step 1.7. If max ≤ i ≤ d w ( X ′ , i ) < ǫ , then add X ′ to R and GOTO Step 1.8. Define m = ⌊ d/ ⌋ . If (1 + δ ) < vol (˜ π ′ m ( X )) vol (˜ π ′ m ( X ′ )) , then define X := X ′ and GOTO Step 3.9. Subdivide X ′ into two pieces, X ′ and X ′ , along a dimension which maximizes w ( X ′ , i ) ,and so that X ′ = X ′ ∪ X ′ . Add the two new cubes to S and GOTO Step 1. Theorem 6.2.
Let S = { X i ⊆ R × ˜Ω s } with M ≥ and s > . Let A , B and R be theoutput of Algorithm 6.1 run with input S and various computational parameters.(i) If F (ˆ x ) = 0 for some ˆ x ∈ S S , then ˆ x ∈ S A ∪ B ∪ R .(ii) For each X ∈ A and α ∈ π α ( X ) , there is a unique ˆ x = ( α, ˆ ω α , ˆ c α ) ∈ X such that F (ˆ x ) = 0 .(iii) For each X ∈ B , if there is a solution ˆ x ∈ X to F = 0 , then ˆ x is on the principalbranch.Proof. We prove the claims of the theorem.(i) Suppose there is some solution ˆ x ∈ X for some X ∈ S . We show that ˆ x ∈ S S∪A∪B∪R at every step of the algorithm. If we replace X by X ′ as in Step 3, then ˆ x ∈ X ′ byTheorem 4.2. In Step 4, if f lag = 1 then in fact X ′ = ∅ , so X could not have containedany solutions in the first place. In Steps 5, 6 and 7, the cube X ′ is added to one of A ,33 or R . Hence, as ˆ x ∈ X ′ then ˆ x ∈ S S ∪ A ∪ B ∪ R . If in Step 8 we decide to prunethe cube X ′ again, then we may repeat the argument made for Steps 3-7. In Step 9we divide X ′ into two new cubes X ′ and X ′ for which X ′ = X ′ ∪ X ′ . Hence ˆ x will becontained in at least one of X ′ or X ′ , and both cubes are added to S , so we cannotlose the solution in Step 9.Thus we have shown that ˆ x ∈ S S ∪ A ∪ B ∪ R at every step. Since the algorithm canonly stop when S = ∅ , it follows that every solution ˆ x initially contained in S S willeventually be contained in S A ∪ B ∪ R .(ii) The only way a cube X ′ can be added to A is in Step 6. That is, for some cube X ∈ S the output of Algorithm 4.1 returned { , X ′ } . Thus, it follows from Theorem 4.2 thatfor all α ∈ π α ( X ) there is a unique ˆ x = ( α, ˆ ω α , ˆ c α ) ∈ X such that F (ˆ x ) = 0.(iii) The only way a cube X ′ can be added to B is in Step 5. That is, for some cube X ∈ S the output of Algorithm 4.1 returned { , X ′ } . Thus, it follows from Lemma 2.9 thatthe only solutions to F = 0 in X ′ are those on the principal branch.If a cube has no zeros inside of it yet there is a solution close to its boundary, thenproving that the cube does not contain any solutions can be very difficult, resulting inan excessive number of subdivisions. This phenomenon is common to branch and boundalgorithms and is referred to as the cluster effect [25]. As we wish to enumerate not justisolated solutions but a 1-parameter family of solutions, the difficulty of the cluster effect ismultiplied. Furthermore, we cannot expect that the boundary of a cube will almost nevercontain a solution. In particular, when we subdivide a cube we may also bisect the curveof solutions, and further subdivisions will not remedy this problem (see Figure 5). As such,we should not expect that R 6 = ∅ .To address this issue we apply Algorithm 6.3 to the output of Algorithm 6.1. In Step 1we recombine cubes in R which overlap in the α dimension. In Step 2 we split the cubes in R along the α -dimension to make them easier to prune, which we do in Step 3. Ideally byStep 4 all of the cubes have been removed from R , having been added to either A or B .Even if R = ∅ at this point, it is not immediately clear that the only solutions are onthe principal branch. For two distinct cubes X , X ∈ A , if there is some α such that α ∈ π α ( X ) and α ∈ π α ( X ), then there could very well be two distinct solutions at theparameter α . In fact, since we subdivide along the α –dimension it is to be expected thata cube will share an α –value with one or two other cubes. In Steps 6-9 of Algorithm 6.3 wecheck to make sure that when two cubes have α –values in common, then there is a uniquesolution associated to each α ∈ π α ( X ) ∩ π α ( X ). Algorithm 6.3.
Take as input sets A , B , R produced by Algorithm 6.1 and a computationalparameter n ∈ N . The output is a pair of intervals I A α , I B α and either success or failure.1. Combine the elements in R whose α -components overlap in more than just a point.That is, for all X, Y ∈ R , if diam ( π α ( X ) ∩ π α ( Y )) > , then replace X and Y in theset R with a new cube Z containing X ∪ Y .2. Subdivide each X ∈ R along the α -dimension.3. For all X ∈ R calculate { f lag, X ′ } = P rune ( n ) ( X ) , the output of Algorithm 4.1iterated at most n times with initial input X . If f lag = 1 , then remove X from R . If .57 1.575 1.58 1.585 1.59 1.595 1.600.050.10.150.2 a a Figure 5:
An example output of Algorithm 6.1. The cubes in A are in green, the cubes in B are in blue,and the cubes in R are in red. f lag = 2 , then remove X from R and add X ′ to B . If f lag = 3 , then remove X from R and add X ′ to A .4. If R 6 = ∅ then return FAILURE.5. Define I A α = S X ∈A π α ( X ) and I B α = S X ∈B π α ( X ) .6. Construct a cover I ′B of the parts of cubes in A which intersect with S B . That is,define I B = { X ∈ A : π α ( X ) ∩ I B α } . Then define I ′B by, for each X ∈ I B , taking the α -component of X and setting it equal to π α ( X ) ∩ I B α and adding the modified cube to I ′B .7. For all X ∈ I ′B calculate { f lag, X ′ } = P rune ( n ) ( X ) , the output of Algorithm 4.1iterated n –times with initial input X . If f lag = 2 then return FAILURE.8. Construct a cover I ′A of the parts of cubes in A which intersect with another cube in A . That is, define I A = { ( X, Y ) ∈ A × A : X = Y, π α ( X ) ∩ π α ( Y ) = ∅} . Then define I ′A by, for each ( X, Y ) ∈ I A , defining a new cube Z which contains X ∪ Y , replacingthe α -component of Z by π α ( X ) ∩ π α ( Y ) , and adding Z to I ′A .9. For all Z ∈ I ′A calculate { f lag, Z ′ } = P rune ( n ) ( Z ) , the output of Algorithm 4.1iterated n –times with initial input Z . If f lag = 3 then return FAILURE.10. If the algorithm reaches this point, return SUCCESS. Theorem 6.4.
Let A , B , R denote the output of Algorithm 6.1 run with input S = { X i ⊆ R × ˜Ω s } where M ≥ and s > . Suppose having received input A , B , R and n ∈ N ,Algorithm 6.3 returns SUCCESS and intervals I A α and I B α .(i) If α ∈ I A α \ I B α , then there is a unique solution ˆ x α = ( α, ˆ ω α , ˆ c α ) ∈ S S such that F α (ˆ ω α , ˆ c α ) = 0 .(ii) If α ∈ I B α , then the only solutions to F α = 0 in S S are on the principal branch.Proof. We describe the first 4 steps of the algorithm and then prove the theorem.35. Let R denote the initial input to the algorithm and R ′ denote the resulting set pro-duced by Step 1. By its construction, it follows that S R ⊆ S R ′ .2. If we subdivide the cubes in R ′ , then it is still true that S R ⊆ S R ′ .3. As described in the proof of Theorem 6.2, if f lag = 1 , , X , add X ′ to B and add X ′ to A . Appropriate, that is, in thesense that the conclusion of Theorem 6.1 will hold for these modified sets A , B and R .4. If we cannot show that every region of phase-space lies in either A or B then we areunable to prove the theorem. Otherwise, every solution to F = 0 in S S is containedin S A ∪ B .We prove claim ( i ). If α ∈ I A α \ I B α there is a solution ˆ x α to F α = 0 in S A . Suppose thereexists a second distinct solution ˆ x ′ α to F α = 0. Since each cube X ∈ A contains a uniquesolution for all α ∈ π α ( X ), there would exist distinct cubes X, Y ∈ A such that ˆ x α ∈ X and ˆ x ′ α ∈ Y . It follows then that there exists some Z ∈ I ′ α such that ˆ x α , ˆ x ′ α ∈ Z . Since itis determined by Step 9 that f lag = 3 in the output of P rune ( n ) ( Z ), therefore by Theorem4.2 there exists a unique solution to F = 0 in Z . Thereby ˆ x α = ˆ x ′ α , and if α ∈ I A α \ I B α , thenthere is a unique solution ˆ x α = ( α, ˆ ω α , ˆ c α ) ∈ S S such that F α (ˆ ω α , ˆ c α ) = 0.We prove claim ( ii ). Suppose there exists some ˆ x α such that α ∈ I B α and F α (ˆ ω, ˆ c ) = 0.Since the algorithm passed through Step 4, it follows that ˆ x α ∈ S A ∪ B . If ˆ x α ∈ S B ,then ˆ x α is on the principal branch by Theorem 6.2. If ˆ x α ∈ S A , then there exists a cube X ∈ I ′B such that ˆ x α ∈ X . If the Algorithm 6.3 is successful, then when Algorithm 4.1is run n –times with initial input X it will produce f lag = 2. Hence by Theorem 4.2 thissolution ˆ x α ∈ S A must be on the principal branch. Proof of Theorem 1.1.
We implemented the algorithms discussed in this paper using MAT-LAB version R2017b (see [9] for the code). The calculations were performed on IntelXeon E5-2670 and Intel Xeon E5-2680 processors, and used INTLAB for the interval arith-metic [24]. A summary of the algorithms’ runtime is given in Table 1.For the intervals I α taking the values (containing at least) [ π , . . , . . , . . , . i = 2, j = 20, n T ime = 32, N P eriod = 10, N P rune = 4, ǫ = 0 .
05 and ǫ = 0 .
05. We then ran Algorithm5.7 using computational parameters M = 10 and S = 3, and N = 15 producing outputs S I α (see Figure 4). By Theorem 5.8, if y is a SOPS at parameter α ∈ I α given as in (4), then( α, ω, c ) ∈ S S I α . By Proposition 2.4 the SOPS to (2) at parameters α ∈ I α are in bijectivecorrespondence with the nontrivial zeros of F inside S S I α .On each of the collections of cubes S I α we ran Algorithm 6.1, using the following com-putational parameters: For the stopping criterion we used ǫ = 0 . α ∈ [ π , .
6] and ǫ = 0 .
01 otherwise. For the continue-pruning criterion, in every case we used δ = 0 . d = 6, corresponding to thevariables α, ω, a ∈ R and c , c ∈ C . For the set of weights, in each case we used λ = 8(corresponding to α ) and λ i = 1 otherwise.The output of Algorithm 6.1 are sets A I α , B I α , R I α . On each of these resulting outputswe ran Algorithm 6.3 using n = 5, and in each case it was successful, producing sets I A α and I B α . When I α = [ π , .
6] then I B α = [ π , π + 0 . I A α = [ π + 0 . , . I A α = I α . By Theorem 6.2, this shows that for all α ∈ [ π +0 . , .
9] there exists36 unique solution to F α = 0 in S S , and if α ∈ [ π , π + 0 . α = π on or offthe principal branch, and there are no folds in the principal branch for α ∈ ( π , π + 0 . α ∈ ( π , .
9] there exists a unique solution to (2). By [10] and [29] there exists aunique SOPS to (2) for α ∈ [1 . , .
0] and α ≥ .
67 respectively. Hence there exists a uniqueSOPS to (2) for all α > π . I α N bf N ′ grid N grid T bf T grid T ∗ bb T verify [ π , .
6] 1604 614 181 602.5 5.2 2 . ∗ . , .
7] 985 861 165 461.6 5.6 4 . ∗ . , .
8] 604 566 143 335.1 3.6 10 . ∗ . , .
9] 292 277 97 135.6 1.9 67 . ∗ Computational benchmarks from the computer-assisted proof of Theorem 1.1. N bf – the numberof bounding functions output by [10, Algorithm 5.1]. N ′ grid – the number of cubes in S ′′′ after Step7 in Algorithm 5.7. N grid – the number of cubes output by Algorithm 5.7. T bf – the run time (min.)of [10, Algorithm 5.1]. T grid – the run time (min.) of Algorithm 5.7. T ∗ bb – the run time (min.) ofAlgorithm 6.1 parallelized using 20 workers. T verify – the run time (min.) of Algorithm 6.3. Proof of Theorem 1.2.
By [14] every global solution to (1) has a positive, integer valued lapnumber V ( x, t ). For non-zero x the lap number will be an odd integer, defined by fixing thesmallest possible σ ≥ t such that x ( σ ) = 0 and defining: V ( x, t ) = ( the x ( s ) in ( σ − , σ ]; or1 if no σ exists . Let us fix x as a periodic solution to (1) with period L . For any t ∈ R the lap number V ( x , t ) remains constant, and we can define N := V ( x , t ). If N = 1 then x must be aSOPS. If N ≥ n := N − and r := 1 − nL . By [14], it follows that2 /N < L < / ( N − < r < N − . Defining x ( t ) := x ( rt ) and α = rα wecalculate the derivative of x ( t ) as: x ′ ( t ) = − α f ( x ( rt − x ( rt −
1) = x ( rt − nL ) = x ( r ( t − x ( t − . Hence it follows that x ′ ( t ) = − α f ( x ( t − V ( x ) ≥ x is a rescaling of a periodic solution x with period length L = L /r >
2. Hence x is arescaling of a SOPS. One pertinent question that remains concerns the period length of SOPS to Wright’s equa-tion.
Conjecture 7.1.
The period length of SOPS to (2) increases monotonically in α . The rigorous numerics in [10,12] strongly suggests this to be true when α ≤
6, and when α ≥ . L satisfies | L − α − e α | < . α − by [20]. It is known that theperiod length increases monotonically when α ∈ ( π , π + 6 . × − ] by [26]. HoweverConjecture 7.1 is unresolved for α > π + 6 . × − .Another question, proposed in [2], is the generalized Wright’s conjecture.37 onjecture 7.2. For every α > the set U ( α ) , the closure of the forward extension by thesemiflow of a local unstable manifold at zero, is the global attractor for (2) . This is known to be true for α ≤ π by [2,26,28] and is unresolved for α > π . Conjecture7.2 can be reduced to a question about the number of rapidly oscillating periodic solutions,and moreover Conjecture 7.1 implies Conjecture 7.2. To wit, by the Poincar´e-Bendixsontheorem for monotone feedback systems [15], the ω -limit set of any initial data to (2) iseither 0 or a periodic orbit. The lap number organizes the attractor into Morse sets S N by [14], and by [6] there is always a connecting orbit from the unstable manifold of theorigin to the Morse set S N . Hence, to prove Conjecture 7.2, it would suffice to show thateach Morse set consists of exactly one periodic orbit.By Theorem 1.2 there are no isolas of periodic orbits, so multiple rapidly oscillatingperiodic solutions can only arise if there is a fold in one of the branches of rapidly oscillatingperiodic solutions. If Conjecture 7.1 holds, then such a fold can be ruled out using rescalingequation in Theorem 1.2. In particular, if there are two SOPS at parameters α < α withperiod lengths L , L and the equality α = α (1 + nL ) = α (1 + nL ) holds, then therewill be two distinct rapidly oscillating periodic solutions at parameter α . This equalitycannot hold if L < L whenever α < α . Thereby Conjecture 7.1 implies Conjecture 7.2.There are still further questions about Wright’s equation. In [16] the authors show asemi-conjugacy of Wright’s equation, and negative feedback systems more generally, onto afamily of finite dimensional ODEs. Outside the dynamics described by this semi-conjugacy,are there any other interesting dynamics in (2)? Furthermore, do the stable and unstablemanifolds of the periodic orbits in (2) intersect transversely?There are many future directions for the rigorous numerics of infinite dimensional dy-namical systems. Perhaps one of the most striking features of Figure 2 and Figure 5 is thenon-uniform size of cubes. This seems to be a result of applying the branch and boundmethod to a 1-parameter family of solutions instead of a collection of isolated solutions.One approach would be to first validate a neighborhood around the branch of solutions ( ´ala [12]) and then use a branch and bound method to ensure that there are no solutions out-side of this neighborhood. In this paper, we used a collection of weights { λ } di =0 to mitigatethis problem. When using all equal weights ( λ i = 1 for all i ), the vast majority of cubesoutput by Algorithm 6.1 ended up in R . Having a better heuristic for deciding along whichdimension to branch would be very useful, particularly so if it does away with the a priori need to select a maximal subdivision dimension d as a computational parameter.Integral to the success of our algorithm (allowing it to finish in finite time) are theestimates derived in [10] which bound all of the slowly oscillating periodic solutions toWright’s equation. Since most initial conditions are attracted to the single SOPS in Wright’sequation, it was sufficient for the methods in [10] to be relatively simple. Future work couldbe done toward bounding all periodic orbits when there are multiple (unstable) solutions,or when the dimension is higher, as well as bounding all periodic solutions to ODEs andPDEs.Another question, explored in [13], is “what the best Banach space to work in?” In thispaper we consider the space Ω s of Fourier coefficients with algebraic decay. In Algorithm5.7, the estimates for obtaining a priori estimates on the Fourier coefficients of SOPS alwaysimprove in absolute terms by using larger value of S . However, the value of C will increasewhen using a larger S . It would likely be beneficial to initially run Algorithm 5.7 with alarge S , and then convert these bounds into a smaller S so that C will shrink as well.38owever, for other applications and other infinite dimensional problems, the question ofwhat is the optimal Banach space remains. Acknowledgments
The author thanks Konstantin Mischaikow for many insightful discussions, as well as JohnMallet-Paret and Roger Nussbaum for discussions on future work.The author acknowledges the Office of Advanced Research Computing (OARC) at Rutgers,The State University of New Jersey for providing access to the Amarel cluster and associatedresearch computing resources that have contributed to the results reported here.
References [1] A. Ambrosetti and G. Prodi.
A primer of nonlinear analysis . Number 34. CambridgeUniversity Press, 1995.[2] B. B´anhelyi, T. Csendes, T. Krisztin, and A. Neumaier. Global attractivity of thezero solution for Wright’s equation.
SIAM Journal on Applied Dynamical Systems ,13(1):537–563, 2014.[3] S.-N. Chow and J. Mallet-Paret. Integral averaging and bifurcation.
Journal of Differ-ential Equations , 26(1):112–159, 1977.[4] T. Csendes and D. Ratz. Subdivision direction selection in interval methods for globaloptimization.
SIAM Journal on Numerical Analysis , 34(3):922–938, 1997.[5] S. Day and W. D. Kalies. Rigorous computation of the global dynamics of integrod-ifference equations with smooth nonlinearities.
SIAM Journal on Numerical Analysis ,51(6):2957–2983, 2013.[6] B. Fiedler and J. Mallet-Paret. Connections between Morse sets for delay differentialequations.
J. reine angew. Math , 397:23–41, 1989.[7] Z. Galias and P. Zgliczy´nski. Infinite dimensional Krawczyk operator for finding pe-riodic orbits of discrete dynamical systems.
International Journal of Bifurcation andChaos , 17(12):4261–4272, 2007.[8] E. Hansen and G. W. Walster.
Global optimization using interval analysis: revised andexpanded , volume 264. CRC Press, 2003.[9] J. Jaquette. MATLAB code available at: .[10] J. Jaquette, J.-P. Lessard, and K. Mischaikow. Stability and uniqueness of slowlyoscillating periodic solutions to Wright’s equation.
Journal of Differential Equations ,263(11):7263–7286, 2017.[11] G. S. Jones. The existence of periodic solutions of f ( x ) = − αf ( x − { f ( x ) } . Journal of Mathematical Analysis and Applications , 5(3):435–450, 1962.3912] J.-P. Lessard. Recent advances about the uniqueness of the slowly oscillating periodicsolutions of Wright’s equation.
Journal of Differential Equations , 248(5):992–1016,2010.[13] J.-P. Lessard and J. D. Mireles James. Computer assisted fourier analysis in sequencespaces of varying regularity.
SIAM Journal on Mathematical Analysis , 49(1):530–561,2017.[14] J. Mallet-Paret. Morse decompositions for delay-differential equations.
Journal ofdifferential equations , 72(2):270–315, 1988.[15] J. Mallet-Paret and G. R. Sell. The Poincar´e–Bendixson theorem for monotone cyclicfeedback systems with delay.
Journal of differential equations , 125(2):441–489, 1996.[16] C. McCord and K. Mischaikow. On the global dynamics of attractors for scalar delayequations.
Journal of the American Mathematical Society , 9(4):1095–1133, 1996.[17] R. E. Moore. A test for existence of solutions to nonlinear systems.
SIAM Journal onNumerical Analysis , 14(4):611–615, 1977.[18] R. E. Moore, R. B. Kearfott, and M. J. Cloud.
Introduction to interval analysis . SIAM,2009.[19] A. Neumaier.
Interval methods for systems of equations , volume 37. Cambridge uni-versity press, 1990.[20] R. Nussbaum. Asymptotic analysis of some functional-differential equations. In C. L.Bednarek, A. R., editor,
Dynamical systems, II , pages 277–301, 1982.[21] R. D. Nussbaum. Periodic solutions of analytic functional differential equations areanalytic.
Michigan Math. J. , 20:249–255, 1973.[22] R. D. Nussbaum. A global bifurcation theorem with applications to functional differ-ential equations.
Journal of Functional Analysis , 19(4):319–338, 1975.[23] B. T. Regala.
Periodic solutions and stable manifolds of generic delay differentialequations . PhD thesis, Brown University, 1989.[24] S. M. Rump. INTLAB–INTerval LABoratory. In
Developments in reliable computing ,pages 77–104. Springer, 1999.[25] H. Schichl and A. Neumaier. Exclusion regions for systems of equations.
SIAM journalon numerical analysis , 42(1):383–408, 2004.[26] J. B. van den Berg and J. Jaquette. A proof of Wright’s conjecture. arXiv preprintarXiv:1704.00029 , 2017.[27] J. B. Van Den Berg and J.-P. Lessard. Chaotic braided solutions via rigorous numerics:Chaos in the Swift–Hohenberg equation.
SIAM Journal on Applied Dynamical Systems ,7(3):988–1031, 2008.[28] E. M. Wright. A non-linear difference-differential equation.
J. reine angew. Math ,194(1-4):66–87, 1955. 4029] X. Xie.
Uniqueness and stability of slowly oscillating periodic solutions of differentialdelay equations . PhD thesis, Rutgers University, 1991.[30] X. Xie. Uniqueness and stability of slowly oscillating periodic solutions of delay equa-tions with unbounded nonlinearity.