Validated numerics for period-tupling and touch-and-go bifurcations of symmetric periodic orbits in reversible systems
VValidated numerics for period-tupling and touch-and-gobifurcations of symmetric periodic orbits in reversible systems
Irmina Walawska ∗ , Daniel Wilczak ∗∗ Faculty of Mathematics and Computer Science, Jagiellonian University, Łojasiewicza 6, 30-348 Krak´ow, Poland.
Abstract
We propose a general framework for computer-assisted verification of the presence of symme-try breaking, period-tupling and touch-and-go bifurcations of symmetric periodic orbits for re-versible maps. The framework is then adopted to Poincar´e maps in reversible autonomous Hamil-tonian systems.In order to justify the applicability of the method, we study bifurcations of halo orbits in theCircular Restricted Three Body Problem. We give a computer-assisted proof of the existenceof wide branches of halo orbits bifurcating from L , , -Lyapunov families and for wide range ofmass parameter. For two physically relevant mass parameters we prove, that halo orbits undergomultiple period doubling, quadrupling and third-order touch-and-go bifurcations. Keywords: bifurcations of periodic orbits, reversible systems, validated numerics, halo orbits.
1. Introduction
In the past 40 years there were proposed very e ffi cient methods for numerical continuationof periodic orbits for general ODEs, or periodic orbits satisfying certain symmetries with specialfocus on applications to Hamiltonian systems [19, 37, 22]. They are implemented in very e ffi cientpackages, such as AUTO [18], MATCONT [17] and CONTENT [27].Most of the methods are similar in the spirit. The family of periodic orbits satisfies certainfinite-dimensional implicit equations, which is solved by a Newton-like scheme. The equationsusually involve period (return time) of the orbit, space variables and / or value of the Hamiltonian.The Newton-like iteration requires computation of monodromy matrix, which is not an issue ∗ The work of the first author was supported by the Polish National Science Center under grant UMO-2016 / / N / ST6 / ∗∗ Corresponding author. This research is partially supported by the Polish National Science Center under MaestroGrant No. 2014 / / A / ST1 / / / B / ST1 / Email addresses:
[email protected] (Irmina Walawska),
[email protected] (Daniel Wilczak)
Preprint submitted to Communications in Nonlinear Science and Numerical Simulation March 5, 2019 a r X i v : . [ m a t h . D S ] M a r owadays, where many e ffi cient tools for integration of ODEs along with variational equationsare available [4, 23, 10].Bifurcations of periodic orbits can be detected by looking at changes of determinant of mon-odromy matrix [54], inspecting so-called stability parameter in low dimensional systems [21, 22],analysis of normal form after Lyapunov–Schmidt reduction [16] or solving an equation specificfor the type of bifurcation [38]. Literature on the topic is really wide, and the list of methods andreferences mentioned above is clearly incomplete.The present paper is complementary to the above results and methods. We focus on reversiblemaps and reversible Hamiltonian systems. The primary result of the paper is a general frameworkfor rigorous computer-assisted verification, that a branch of symmetric periodic orbits undergoesperiod-tupling and / or touch-and-go bifurcation. Finding an approximate numerical candidate forbifurcation point is just a preliminary step of the validation algorithm and for this purpose we canuse any of the methods mentioned above [16, 22, 21, 38, 54]. Then, checking several inequali-ties on the (Poincar´e) map under consideration and its derivatives on an explicit neighbourhood(input to the algorithm) of the approximate bifurcation point we can prove, that it contains a bi-furcation point and we can conclude about type of bifurcation. As an output of the algorithmwe obtain guaranteed bounds on both the bifurcation point and the parameter of the system, atwhich bifurcation occurs.In order to justify applicability of the framework we apply it to the Circular RestrictedThree Body Problem (CR3BP) [44]. We give a computer-assisted proof of the existence ofwide branches of so-called halo orbits bifurcating from the families of planar Lyapunov orbitsaround L , , libration points. We also prove, that for some physically relevant mass parametersof the system, these branches undergo period-tupling and touch-and-go bifurcations. These rig-orous results are justification of some numerical observations from previous articles, in particular[21, 19]. We also would like to emphasize, that we have found a new phenomenon regarding bi-furcations of halo orbits near L libration point, when the mass parameter of the system tends tozero. We have observed, that the energy at the symmetry breaking bifurcation, which creates haloorbit, as a function of the mass parameter is not monotone, when the relative mass tends to zero— see Remark 19. Although not validated, this observation has been made using nonrigorousnumerics with very high accuracy.Rigorous numerical investigation of periodic orbits to ODEs became quite standard [2, 9,11, 24, 25, 41, 45, 48, 49, 50]. To the best of our knowledge, there are very few results re-garding computer-assisted verification of bifurcations of periodic orbits of ODEs, and the fieldremains widely open. Validation, that a family of periodic orbits undergoes a bifurcation usuallyinvolves computation of rigorous bounds on higher order derivatives of the trajectories with re-spect to initial condition (except some special cases when the system admits additional structure).The algorithm capable to do that appeared just in 2011 [53] (implemented as a part of publiclyavailable CAPD library [10]) and to the best of our knowledge there are no other publicly avail-able algorithms for rigorous integration of higher order variational equations. The C r -Lohneralgorithm from [53] has been already applied to study period-dubling bifurcations of periodic or-bits [51] in the R¨ossler system [43], homoclinic tangencies of periodic orbits in a time-periodicforced-damped pendulum equation [52] and non-local cocoon bifurcations [28] in the Michelsonsystem [32]. Very recently [6], another approach to compute periodic orbits for ODEs withoutintegration of the system has been proposed. Periodic solutions are approximated via piecewiseChebyshev polynomials and then their existence is validated by analysis of certain nonlinear op-erator on a Banach space of Chebyshev coe ffi cients. The e ffi ciency of the method is illustrated onthe example of Equilateral Circular Restricted Four Body Problem. Similar approach has been2uccessfully applied for validated computation of bifurcations of equilibria for ODEs and steadystates for PDEs — see for example [8, 30, 31]. A di ff erent, geometric method for validation ofbifurcations of steady states in the Kuramoto-Sivashinsky PDE is also proposed in [55].Our algorithm for proving the existence of period-tupling and touch-and-go bifurcation issimilar in the spirit to that proposed in [51, 55], but exploits the presence of a reversing symmetryof the system. After fixing an appropriate coordinate system in a neighbourhood of an apparentbifurcation points, we perform validated Lyapunov-Schmidt reduction. In this way, the analysisof bifurcation is transformed to zero-finding problem of a bivariate scalar-valued function, called bifurcation function . Then, checking some inequalities on derivatives of the bifurcation functionwe can prove, that its set of zeroes is the union of two smooth curves which intersect at an uniquepoint — the bifurcation point. Finally, some non-degeneracy conditions (inequalities on higherorder derivatives of the bifurcation function) let us to conclude about the type of bifurcation. Inthis paper we restrict to symmetry breaking, period-tupling and touch-and-go bifurcations.The article is organized as follows. In Section 2 we introduce notation and main definitionsused in the paper. Theoretical results, which constitute a basis of the computational frameworkfor general reversible (Poincar´e) maps are presented in Section 3 and then adopted to autonomousreversible Hamiltonian systems in Section 4. Finally, in Section 6 we give formal statements oftheorems regarding continuation and bifurcations of halo orbits in the CR3BP.In the Appendix we present two auxiliary yet new results, which are included mainly forself-consistency of the paper. Frequently we have to show, that a solution to an implicit equationis defined over an explicit domain. This step appears for example in continuation of periodicorbits or in validated Lyaunov-Schmidt reduction. Although there are available methods for thispurpose (see for example [7, 20, 53, 9]), in Appendix A we provide an improvement, whichtakes advantage from higher order derivatives and flattening the implicit function by a smooth,well chosen substitution. In Appendix B we propose own and short algorithm for finding anapproximate bifurcation point. It takes advantage from the Lyapunov–Schmidt reduction andcomputation of higher order derivatives of the (Poincar´e) map under consideration. Using it, wecould easily localize approximate bifurcation points in the CR3BP with the accuracy 10 − .
2. Preliminaries
For a map f : D ⊂ X → X , a predicate C : X → { true , false } and a set U ⊂ D we introducethe following notation Fix( f , U , C ) = { x ∈ U : ( f ( x ) = x ) ∧ C ( x ) } . Thus, Fix( f , U , C ) is the set of fixed points of f in U satisfying constraint C . Although thereis a natural correspondence between subsets of X and univariate predicates defined on X , wewill separate U and C to emphasize rather rare property C in a larger (usually open) set U . Wewill also write Fix( f , U ) if C ( x ) ≡ true . For a set M ⊂ R × X and ν ∈ R , we define its slice by M ν = { x ∈ X : ( ν, x ) ∈ M } . For a map f : M ⊂ R × X → X and fixed ν ∈ R we define f ν : M ν → X by f ν ( x ) = f ( ν, x ). Definition 1.
A homeomorphism R : X → X is called a reversing symmetry for f : M ⊂ X → Xif ( R ◦ f )( M ) ⊂ M and for x ∈ M there holds ( R − ◦ f ◦ R ◦ f )( x ) = x .
3t is easy to see that any reversing symmetry R for f is also a reversing symmetry for all iterations f n defined on their proper domains. Definition 2.
A homeomorphism S : X → X is called a symmetry for f : M ⊂ X → X if S ( M ) ⊂ M and for x ∈ M there holds S ( f ( x )) = f ( S ( x )) . For a map S : X → X we define a predicate C S : X → { true , false } by C S ( x ) = ( x ∈ Fix( S , X )) . (1)We will use this notation to select points satisfying certain symmetries. Our primary object of interest is a C -smooth function f : M → X defined on an open set M ⊂ R × X , where X is a smooth manifold. In this article we make the following standingassumption: C1: for ν ∈ R , if M ν (cid:44) ∅ then the function f ν : M ν → X is a di ff eomorphism onto image.The above assumption fits the applications we have in mind, i.e. f ν will be a one-parameterfamily of Poincar´e maps. Thus, the domain of each map f ν may vary with the parameter ν .Bifurcations are usually defined by their normal forms [13]. The following f k ν ( x ) = x ( ν − x ) , (2) f k ν ( x ) = x ( ν − x ) (3)describe period-tupling and transcritical bifurcations, respectively. When an eigenvalue of thederivative of a reversible planar map at a symmetric fixed point crosses 1 : k resonance, k ≥ f k ν ( x ) − x = Definition 3.
Let C : X → { true , false } be a predicate and let k be a positive integer. Wesay that f ν : M ν → X has period k –tupling bifurcation at ( ν ∗ , x ∗ ) ∈ M, if there exists V = [ ν , ν ] × U ⊂ M, such that ( ν ∗ , x ∗ ) ∈ int V and there are continuous and smooth in the interiorof their domains functionsx f p : [ ν , ν ] −→ int U , x b , x b : [ ν ∗ , ν ] −→ int U , such that x b ( ν ∗ ) = x b ( ν ∗ ) = x f p ( ν ∗ ) = x ∗ and the following conditions are satisfied. ix ( R ) ν < ν * Fix ( R ) ν = ν * Fix ( R ) ν > ν * Fix ( R ) ν < ν * Fix ( R ) ν = ν * Fix ( R ) ν > ν * Figure 1: Upper row: period quadrupling bifurcation of a reversible planar map. After bifurcation a period-4 orbit iscreated, which intersects Fix( R ) at exactly two points. Lower row: third-order touch-and-go bifurcation of a reversiblemap. x fp ( ν ) x b ( ν ) x b ( ν ) V ν * ν ν ν x * x x fp ( ν ) x b ( ν ) V ν * ν ν ν x * x Figure 2: Geometry of (left) period-tupling and (right) touch-and-go bifurcations. Periodic points:
Fix( f p ν , U , C ) = { x f p ( ν ) } , ν ∈ [ ν , ν ∗ ] , ≤ p ≤ k , Fix( f p ν , U , C ) = { x f p ( ν ) } , ν ∈ [ ν ∗ , ν ] , ≤ p < k , Fix( f k ν , U , C ) = { x f p ( ν ) , x b ( ν ) , x b ( ν ) } , ν ∈ [ ν ∗ , ν ] , f k ν , U , C ) = , ν ∈ ( ν ∗ , ν ] . If k is even, then for ν ∈ ( ν ∗ , ν ] there holdsf k / ν ( x b ( ν )) = x b ( ν ) , f k / ν ( x b ( ν )) = x b ( ν ) . Definition 4.
Let C : X → { true , false } be a predicate and let k be a positive integer. Wesay that f ν : M ν → X has k th –order touch-and-go bifurcation at ( ν ∗ , x ∗ ) ∈ M, if there existsV = [ ν , ν ] × U ⊂ M, such that ( ν ∗ , x ∗ ) ∈ int V and there are continuous and smooth in theinterior of their domains functions x f p , x b : [ ν , ν ] −→ int U , satisfying the following conditions. Periodic points:
Fix( f p ν , U , C ) = { x f p ( ν ) } , ν ∈ [ ν , ν ] , ≤ p < k , Fix( f k ν , U , C ) = { x f p ( ν ) , x b ( ν ) } , ν ∈ [ ν , ν ] . The curves x f p and x b intersect at exactly one point x ∗ = x f p ( ν ∗ ) = x b ( ν ∗ ) . The above two definitions are slightly more general than the normal forms (2)–(3). Indeed, it iseasy to see that the functions f k ν ( x ) = x ( ν − x )( ν n + x m ) , f k ν ( x ) = x ( ν − x )( ν n + x m )also satisfy geometric conditions from Definition 3 and Definition 4 for any n , m ∈ N .
3. Validation of bifurcations of symmetric periodic orbits: reversible maps
Let f : M → X be a family of R –reversible maps satisfying our standing assumption C1 . Inthis section we propose a general framework for computer-assisted verification that a branch of R -symmetric period-2 points for f ν undergoes symmetry breaking, period-tupling or touch-and-go bifurcation. Let V ⊂ M be an explicit neighbourhood of an approximate bifurcation point(both are input to the algorithm). The algorithm is split into two steps.In the first step we validate, that the set of points ( ν, x ) ∈ V , such that x ∈ Fix( R ) is eitherperiod-2 k point or period-2 point for f ν is a union of two regular curves intersecting at exactlyone point. This is a common step for period-tupling and touch-and-go bifurcations and it will bedescribed in Section 3.1. In the case of symmetry breaking bifurcation, we need an additionalconstraint, which allows to isolate the primary curve of period-2 points and the branching-o ff curve of period-2 points. In Section 3.3 we will show, that such a problem can be solved, if thesystem has an additional symmetry, which commutes with R .The second step is specific for each type of bifurcation. We provide checkable by means ofrigorous numerics conditions, which guarantee, that the two curves intersect in a manner givenby Definition 3 and Definition 4.Finally, in Section 3.4 we show that some additional non-degeneracy conditions lead to stan-dard unfolding of the computed bifurcations to their normal forms (2) and (3).6 .1. Bifurcation as two intersecting curves We are mostly interested in studying bifurcations of symmetric periodic orbits of R –reversibleODEs. In continuous-time R –reversible dynamical systems, a trajectory is R –symmetric, if it iseither an equilibrium belonging to Fix( R ) or it intersects twice the set Fix( R ) — see [29]. Inorder to obtain isolated symmetric periodic orbits we assume, that dim X = n and Fix( R ) is an n -dimensional submanifold.Let us fix k ≥ ν, ˆ x ) is a good numerical approximation of a bifurcationpoint, i.e. one of the eigenvalues of D f ˆ ν ( ˆ x ) is close to 1 : k resonance. Since all types ofbifurcations we are studying are local phenomena, we assume, that there is an open interval J and an open set U ⊂ X , such that (ˆ ν, ˆ x ) ∈ J × U and for ν ∈ J the function f k ν is defined on U .We also assume that there are local coordinates in U in whichFix( R , U ) = { ( p , q ) ∈ U : p , q ∈ R n , q = } . We will use the same symbols ( p , q ) as a coordinate system near f ˆ ν ( ˆ x ). Solutions to π q ( f k ν ( p , = R –symmetric periodic points of f ν with principal period not larger than 2 k . The solution setto (4) near (ˆ ν, ˆ x ) is expected to be a union of two curves, one of which corresponds to the fixedpoints of f ν and the second corresponding to period-2 k points of f ν . Remark 5.
There are rigorous numerical methods for validation, that a solution set to an im-plicit equation forms a regular curve over an explicit domain — see for example [9, 7, 51]. In thispaper we propose slightly improved version of the Interval Newton Operator [1, 34, 39]. Sinceit is only an auxiliary step of the main algorithm for validation of bifurcations, we postpone thisimprovement to Appendix A.
The case k = k >
1, we may apply the method described in Appendix A to solve(4) on an explicit range of parameter values ν ∈ J . In what follows we assume, that there is asmooth curve J (cid:51) ν → x f p ( ν ) = ( p ( ν ) , p ( ν ) , ∈ R × R n − × R n , (5)such that for ν ∈ J and i = , . . . , k − f i ν , U , C R ) = { x f p ( ν ) } . (6)In most cases, the function x f p cannot be computed exactly, but using rigorous numerics we canprove, that it exists and we can find bounds for x f p ( ν ) and its derivatives. Although it is notrequired, it is desirable for further numerical computation to choose the coordinate system inFix( R , U ) and the decomposition p = ( p , p ) so that p (cid:48) (ˆ ν ) ≈ R -symmetric, period-2 k points, which solves(4) goes as in [51]. At first, we perform the Lyapunov-Schmidt reduction [13]. We split q = ( q , q ) ∈ R × R n − and using the method from Appendix A we solve for a function p = p ( ν, p )satisfying implicit equation π q ( f k ν ( p , p , = . (7)Assuming that p = p ( ν, p ) is locally unique solution to (7) in J × U , we can define so-called bifurcation function G k ( ν, p ) = π q ( f k ν ( p , p ( ν, p ) , . (8)7n this way we reduced the problem of finding zeros of (4) to a problem of finding zeros of abivariate scalar-valued function (8). By the construction, the function G k vanishes at ( ν, p ( ν )),for ν ∈ J — see (5). Therefore, we can factorize it in the following way G k ( ν, p ) = ( p − p ( ν )) g k ( ν, p ) , (9)where g k ( ν, p ) = (cid:90) ∂ G k ∂ p ( ν, t ( p − p ( ν )) + p ( ν )) dt . (10)The function g k is called the reduced bifurcation function . Recall, that in most cases the function ν → p ( ν ) is unknown and we have only a rigorous bound on it and its derivatives. Similarly as p , the function g k cannot be computed exactly. However, it possible to bound values and partialderivatives of g k using integral representation (10).In the second step, we apply the method from Appendix A to prove, that the set of zeroes of g k can be parametrized as a smooth curve p → ν ( p ). To sum up, in addition to C1 , we makethe following standing assumptions, which are all checkable be means of rigorous numericalmethods: C2: P is a closed interval and P is closed set, such that p ( ν ) ∈ int P for ν ∈ J and thesolution set to (7) in J × P × P is a graph of smooth function p : J × P → int P ; C3: there exists a smooth function P (cid:51) p → ν ( p ) ∈ J , such that { ( ν, p ) ∈ J × P : g k ( ν, p ) = } = { ( ν ( p ) , p ) : p ∈ P } , where g k is defined by (10); C4: there holds 0 (cid:60) ∂ G k ∂ν∂ p ( J × P ) + ∂ G k ∂ p ( J × P ) p (cid:48) ( J ) . We will finish this section by showing, that assumptions
C1–C3 imply, that the two curves, whichsolve (4) intersect at some point, while C4 implies that this intersection is unique. Lemma 6.
Let s : A → B and t : B → A be continuous, where A is an open interval, B isa closed interval and s ( A ) ⊂ int B. Then there exists ( a ∗ , b ∗ ) ∈ int ( A × B ) , such that ( a ∗ , s ( a ∗ )) = ( t ( b ∗ ) , b ∗ ) . Proof:
Since s ( A ) ⊂ int B , the set ( A × B ) \ { ( a , s ( a )) : a ∈ A } has two connected componentscontaining the points ( t (min B ) , min B ) and ( t (max B ) , max B ), respectively. Since t is continuous,there is b ∗ ∈ int B such that ( t ( b ∗ ) , b ∗ ) = ( a ∗ , s ( a ∗ )) for some a ∗ ∈ A . (cid:3) Lemma 7.
Under assumptions
C1–C4 , there is a unique point ( ν ∗ , p ∗ ) ∈ J × P , such that ν ∗ = ν ( p ∗ ) and p ∗ = p ( ν ∗ ) . Proof:
From
C2–C3 and Lemma 6 the curves
J (cid:51) ν → p ( ν ) ∈ P and P (cid:51) p → ν ( p ) ∈ J intersect at some point ( ν ∗ , p ∗ ). Assume that (ˆ ν, p (ˆ ν )) is another intersection point for someˆ ν ∈ J . For τ ∈ [0 ,
1] we define w ( τ ) = ν ∗ + τ (ˆ ν − ν ∗ ) and u ( τ ) = g ( w ( τ ) , p ( w ( τ )) . C3 we have u (0) = u (1) = = u (1) − u (0) = S (ˆ ν − ν ∗ ), where S = (cid:90) ∂ g ∂ν ( u ( τ )) + ∂ g ∂ p ( u ( τ )) p (cid:48) ( w ( τ )) d τ. (11)We will argue, that S ∈ ∂ G k ∂ν∂ p ( J × P ) + ∂ G k ∂ p ( J × P ) p (cid:48) ( J ) and thus by C4 there must be ˆ ν = ν ∗ .For ( ν, p ) ∈ J × P we have ∂ g ∂ p ( ν, p ) = (cid:90) ∂ G ∂ p ( ν, t ( p − p ( ν )) + p ( ν )) tdt and ∂ g ∂ν ( ν, p ) = (cid:90) ∂ G ∂ν∂ p ( ν, t ( p − p ( ν )) + p ( ν )) dt + (cid:90) ∂ G ∂ p ( ν, t ( p − p ( ν )) + p ( ν )) p (cid:48) ( ν )(1 − t ) dt . In particular, for fixed τ ∈ [0 ,
1] and ( ν, p ) = u ( τ ), both second partial derivatives of G areevaluated at a point which does not depend on t . Therefore ∂ g ∂ p ( u ( τ )) = (cid:90) ∂ G ∂ p ( u ( τ )) tdt = ∂ G ∂ p ( u ( τ )) . (12)Similarly ∂ g ∂ν ( u ( τ )) = (cid:90) ∂ G ∂ν∂ p ( u ( τ )) + ∂ G ∂ p ( u ( τ )) p (cid:48) ( w ( τ ))(1 − t ) dt = ∂ G ∂ν∂ p ( u ( τ )) + ∂ G ∂ p ( u ( τ )) p (cid:48) ( w ( τ )) . (13)Using (11)–(13) we obtain S = (cid:90) ∂ G ∂ν∂ p ( u ( τ )) + ∂ G ∂ p ( u ( τ )) p (cid:48) ( w ( τ )) d τ ∈ ∂ G k ∂ν∂ p ( J × P ) + ∂ G k ∂ p ( J × P ) p (cid:48) ( J ) . (cid:3) Assumptions
C1–C4 are easily checkable by means of rigorous numerics. The curve of R -symmetric period-2 k points is defined on an explicit domain P , which makes it possible forfurther continuation of this branch by the standard methods [9, 7, 20, 51] (see also Appendix A).On the other hand, Lemma 7 guarantees, that there is an unique bifurcation point in the domainunder consideration. The type of bifurcation, however, is unknown.In what follows, we derive conditions, which along with C1–C4 guarantee, that period tu-pling or touch-and-go bifurcation occurs in the sense of Definition 3 and Definition 4, respec-tively. 9 heorem 8.
Under assumptions
C1–C4 , if k > is odd and (cid:60) ∂ G k ∂ p ( J × P ) (14) then f ν has k th -order touch-and-go bifurcation at some point ( ν ∗ , x ∗ ) ∈ J × P × P × { } . Proof:
From Lemma 7 there is a unique intersection point ( ν ∗ , p ∗ ) of two curves of zeroes of G k in J × P . From (14) it follows that ∂ g k ∂ p ( ν ∗ , p ∗ ) = (cid:90) ∂ G k ∂ p ( ν ∗ , p ∗ ) tdt = ∂ G k ∂ p ( ν ∗ , p ∗ ) (cid:44) . (15)Di ff erentiating the identity g k ( ν ( p ) , p ) ≡ ∂ g k ∂ν ( ν ∗ , p ∗ ) ν (cid:48) ( p ∗ ) + ∂ g k ∂ p ( ν ∗ , p ∗ ) = . From the above and (15) we conclude, that ν (cid:48) ( p ∗ ) (cid:44)
0. Therefore, there is an interval [ ν , ν ] ⊂J , such that ν ∗ = ν ( p ∗ ) ∈ ( ν , ν ) and the inverse function (cid:101) x b = ν − is defined on it. Since p ∗ ∈ int P , the interval [ ν , ν ] can be chosen so that (cid:101) x b ([ ν , ν ]) ⊂ int P .Define U ε = P × P × ( − ε, ε ) n . We assumed, that Fix( R ) is a submanifold given locally by q =
0. Therefore, there exists su ffi ciently small ε >
0, such that for ν ∈ [ ν , ν ] we haveFix( f k ν , U ε , C R ) = Fix( f k ν , P × P × { } , C R ) = { x f p ( ν ) , x b ( ν ) } , where x b ( ν ) = ( (cid:101) x b ( ν ) , p ( ν, (cid:101) x b ( ν )) , ∈ U ε . From (6) it follows that Fix( f i ν , U ε , C R ) = { x f p ( ν ) } for i = , . . . , k −
1. Thus, all requirements from Definition 4 are satisfied on the set [ ν , ν ] × U ε and for the point x ∗ = ( p ∗ , p ( ν ∗ , p ∗ ) , (cid:3) The next theorem provides a framework for computer-assisted verification of the presence ofperiod k -tupling bifurcation. Theorem 9.
Assume
C1–C4 . If k ≥ is even and ∂ G k ∂ p ( J × P ) · (cid:32) ∂ G k ∂ν∂ p ( J × P ) (cid:33) − ⊂ R − , (16) then f ν has period k-tupling bifurcation at some point ( ν ∗ , x ∗ ) ∈ J × P × P × { } . Proof:
From Lemma 7 we know, that the curves
J (cid:51) ν → p ( ν ) ∈ P and P (cid:51) p → ν ( p ) ∈ J intersect at exactly one point ( ν ∗ , p ∗ ). We will prove, that p ∗ is a proper local minimum of p → ν ( p ). First we will show, that ν (cid:48) ( p ∗ ) = x ∗ = ( p ∗ , p ( ν ∗ , p ∗ ) , { p n } n ∈ N , such that p n (cid:44) p ∗ andlim n →∞ p n = p ∗ . Define ν n = ν ( p n ) and x n = ( p n , p ( ν n , p n ) , p n (cid:44) p ∗ and the inter-section point is unique it follows, that x n is period-2 k point. Since p ∗ ∈ int P , k is even and10 ∗ is a fixed point for f ν ∗ we conclude, that for n large enough the point y n = f k ν n ( x n ) satisfies π p y n ∈ int P and ν ( π p ( y n )) = ν n , lim n →∞ π p ( y n ) = π p ( x ∗ ) = p ∗ . The point y n is an R -symmetric, 2 k -periodic point for f ν n , so it solves (4). Therefore y n ∈ P × P ×{ } . Since x n (cid:44) y n , assumption C2 implies, that also p n (cid:44) π p y n . By the Rolle’s Theorem there isa point in the interval joining p n and π p y n at which ν (cid:48) vanishes. In consequence, arbitrary closeto p ∗ there is a point at which ν (cid:48) is zero, which by the smoothness of ν (cid:48) ( p ) implies ν (cid:48) ( p ∗ ) = ν (cid:48)(cid:48) ( p ∗ ) >
0. Di ff erentiating g k ( ν ( p ) , p ) ≡ ∂ g k ∂ p ( ν ∗ , p ∗ ) = − ∂ g k ∂ν ( ν ∗ , p ∗ ) ν (cid:48) ( p ∗ ) = . On the other hand 0 = ∂ g k ∂ p ( ν ∗ , p ∗ ) = (cid:90) ∂ G k ∂ p ( ν ∗ , p ∗ ) tdt = ∂ G k ∂ p ( ν ∗ , p ∗ ) . Therefore ∂ g k ∂ν ( ν ∗ , p ∗ ) = (cid:90) ∂ G k ∂ p ∂ν (cid:0) ν ∗ , p ∗ (cid:1) + ∂ G k ∂ p (cid:0) ν ∗ , p ∗ (cid:1) (1 − t ) p (cid:48) ( ν ∗ ) dt = ∂ G k ∂ p ∂ν (cid:0) ν ∗ , p ∗ (cid:1) . (17)We also have ∂ g k ∂ p ( ν ∗ , p ∗ ) = (cid:90) ∂ G k ∂ p (cid:0) ν ∗ , p ∗ (cid:1) t dt = ∂ G k ∂ p ( ν ∗ , p ∗ ) . (18)Di ff erentiating twice g k ( ν ( p ) , p ) ≡ p , we obtain ∂ g k ∂ν ( ν ∗ , p ∗ ) ν (cid:48)(cid:48) ( p ∗ ) + ∂ g k ∂ν ( ν ∗ , p ∗ ) (cid:0) ν (cid:48) ( p ∗ ) (cid:1) + ∂ g k ∂ p ∂ν ( ν ∗ , p ∗ ) ν (cid:48) ( p ∗ ) + ∂ g k ∂ p ( ν ∗ , p ∗ ) ≡ . (19)Using ν (cid:48) ( p ∗ ) = ν (cid:48)(cid:48) ( p ∗ ) = − ∂ g k ∂ p ( ν ∗ , p ∗ ) (cid:32) ∂ g k ∂ν ( ν ∗ , p ∗ ) (cid:33) − > . This proves, that the function p → ν ( p ) has a proper local minimum at p ∗ .Now, we will construct the functions x b , x b as required by Definition 3. Since ν (cid:48)(cid:48) ( p ∗ ) > P in such a way, that p ∗ ∈ int P , ν (min P ) = ν (max P ) and ν is strictly convexon P . Given that p ( ν ) is smooth (so it has bounded slope near ν ∗ ) and ν (cid:48) ( p ∗ ) = ν ∈ [ ν ∗ , ν (min P )] there holds p ( ν ) ∈ int P . Shrinking J from below, if necessary, wemay also assume that p ( ν ) ∈ int P for ν < ν ∗ and ν ∈ J .Let us fix ν ∈ ( ν ∗ , ν (min P )). The function ν ( p ) is monotone on [min P , p ∗ ] and [ p ∗ , max P ].Therefore, for (cid:101) ν ∈ [ ν ∗ , ν ] { p : ν ( p ) = (cid:101) ν } = { (cid:101) x b ( (cid:101) ν ) , (cid:101) x b ( (cid:101) ν ) } (cid:101) x b i can be chosen to be continuous and smooth in ( ν ∗ , ν ). Define U ε = P × P × ( − ε, ε ) n . For some su ffi ciently small ε > x b , x b : [ ν ∗ , ν ] → int U ε given by x b i ( ν ) = (cid:0)(cid:101) x b i ( ν ) , p ( ν, (cid:101) x b i ( ν )) , (cid:1) . satisfy conditions from Definition 3. From (6) it follows thatFix( f i ν , U ε , C R ) = { x f p ( ν ) } for i = , . . . , k −
1, which completes the proof. (cid:3)
A special type of bifurcation covered by Definition 3 is when k =
1. In the case of reversiblemaps, we are looking for the set of fixed points of f ν , which near the bifurcation point is expectedto be a union of two intersecting curves. Thus, the principal branch of fixed points of f ν cannot beisolated by applying Interval Newton Operator (Appendix A) to (4). Therefore, we need an extraconstraint in order to isolate two curves that intersect at the bifurcation point. A possible scenariofor this kind of bifurcation is breaking some symmetry S , which provides desired additionalconstraint C S , as defined by (1).Let f : M ⊂ J × X → X be a family of R –reversible and S –symmetric maps. As in Sec-tion 3.1, we focus on computation of the set of R –symmetric fixed points for f ν .Let (ˆ ν, ˆ x ) ∈ M be an apparent bifurcation point and let U ⊂ X be an open set, such that ˆ x ∈ U and J × U ⊂ M . We assume, that there are local coordinates in U , such thatFix( R , U ) = { ( p , q ) ∈ U : p , q ∈ R n , q = } . We impose that Fix( R , U ) and Fix( S , U ) intersect transversally. In the above settings, we expectthat for ν ∈ J the set of double-symmetric fixed points of f ν in U is a single pointFix( f ν , U , C R ∧ C S ) = { x f p ( ν ) } . Thus, using the method described in Appendix A, it is possible to isolate the curve
J (cid:51) ν → x f p ( ν ) = ( p ( ν ) , p ( ν ) , ∈ R × R n − × R n (20)from the branching-o ff curves x b , ( ν ) ∈ Fix( R , U ), which intersect x f p ( ν ) at exactly one pointlocated in Fix( R , U ) ∩ Fix( S , U ).The remaining steps in validation of the existence of a symmetry breaking bifurcation go asdescribed in Section 3.1 and Section 3.2. We perform the Lyapunov-Schmidt reduction (7) andwe define the bifurcation functions G and g by (8) and (10), respectively.In Theorem 9 we assumed that the period of branching-o ff orbits is even. This lead us toa conclusion, that after half of iterations, the points su ffi ciently close to the bifurcation pointcome back to its vicinity. This allowed us to prove, that the reduced bifurcation function g k hasa local minimum at the bifurcation point, which was the crucial step in the proof of the presenceof period-tupling bifurcation.In the case of the symmetry breaking bifurcation we cannot use this argument, because k = R and S yields to the same conclusion.12 heorem 10. Let S and R be a symmetry and a reversing symmetry for f ν , respectively. Assumethat C1–C4 are satisfied for k = and with x f p given by (20). If the symmetries R and Scommute and (16) is satisfied, then f ν has symmetry breaking bifurcation at some point ( ν ∗ , x ∗ ) ∈J × P × P × { } . Proof:
From Lemma 7 there is an unique ( ν ∗ , p ∗ ) ∈ int J × P , such that p ( ν ∗ ) = p ∗ and ν ( p ∗ ) = ν ∗ . Let us set x ∗ : = ( p ( ν ∗ ) , p ( ν ∗ ) , ∈ Fix( S ) ∩ Fix( R ).We will show, that ν (cid:48) ( p ∗ ) =
0. Take a sequence (cid:110) p n (cid:111) n ∈ N of points from P , which convergesto p ∗ and such that p n (cid:44) p ∗ for n ∈ N . Let us define ν n = ν ( p n ), x n = ( p n , p ( ν n , p n ) , ∈ Fix( R )and y n = S ( x n ). By the commutativity of R and S we have R ( y n ) = R ( S ( x n )) = S ( R ( x n )) = S ( x n ) = y n , hence y n ∈ Fix( R ). Similarly f ν n ( y n ) ∈ Fix( R ), because R ( f ν n ( y n )) = R ( f ν n ( S ( x n ))) = R ( S ( f ν n ( x n ))) = S ( R ( f ν n ( x n ))) = S ( f ν n ( x n )) = f ν n ( y n ) . This shows, that y n solves (4) with ν = ν n and ν ( π p ( y n )) = ν n . Sincelim n →∞ π p ( y n ) = π p ( S ( x ∗ )) = p ∗ (21)we conclude, that y n ∈ P × P × { } for n large enough.Since p n (cid:44) p ∗ and the intersection point is unique, it follows that x n (cid:60) Fix( S ), so x n (cid:44) y n .This and C2 imply that also p n (cid:44) π p y n . By the Rolle’s Theorem there is a point in the intervaljoining p n and π p y n at which ν (cid:48) vanishes. Given that both p n and π p ( y n ) (see (21)) converge to p ∗ we conclude, that arbitrarily close to p ∗ there is a point at which ν (cid:48) is zero and thus ν (cid:48) ( p ∗ ) = ν (cid:48)(cid:48) ( p ∗ ) >
0, which proves that p ∗ is a proper minimum of ν ( p ). This makes it possible to define twobranches x b and x b , as required in Definition 3. (cid:3) Theorem 8 and Theorem 9 provide frameworks for computer-assisted verification, that thetwo types of bifurcations occurs in the sense of geometric conditions given by Definition 4 andDefinition 3, respectively. In this section we will show, that the same geometric assumptions leadto standard unfolding of these bifurcations.
Theorem 11.
Under assumptions of Theorem 9 (
C1–C4 and (16)), there is a smooth ν -dependentsubstitution p = p ( ν, w ) , which brings the bifurcation function (8) to the normal formG k ( ν, w ) = α w ( ν − β w ) + O ( | w | ) , (22) where α = α ( ν ) does not vanish at ν = and β = β ( ν ) is positive at ν = . Proof:
Let ν ∗ be the parameter value at which bifurcation occurs. Without loosing generality wemay assume, that ν ∗ =
0. From the proof of Theorem 9 we have that ν (cid:48) ( p ∗ ) =
0. Di ff erentiationof the identity g k ( ν ( p ) , p ) ≡ ∂ g k ∂ p ( ν ∗ , p ∗ ) = − ∂ g k ∂ν ( ν ∗ , p ∗ ) ν (cid:48) ( p ∗ ) = . (23)13irst, a ν -dependent substitution p = z + p ( ν ) brings G k to the form G k ( ν, z ) = zg k ( ν, z ) = zg k ( ν, z + p ( ν )) . (24)Put z ∗ =
0. Di ff erentiation of the product G k ( ν, z ) = zg k ( ν, z ) gives ∂ G k ∂ z ( ν ∗ , z ∗ ) = g k ( ν ∗ , z ∗ ) + z ∗ ∂ g k ∂ z ( ν ∗ , z ∗ ) = , ∂ G k ∂ z ( ν ∗ , z ∗ ) = ∂ g k ∂ z ( ν ∗ , z ∗ ) + z ∗ ∂ g k ∂ z ( ν ∗ , z ∗ ) = , (25)because g k ( ν ∗ , z ∗ ) = ∂ g k ∂ z ( ν ∗ , z ∗ ) = ∂ g k ∂ p ( ν ∗ , p ∗ ) =
0. The non-degeneracycondition (16) and C4 imply that ∂ G k ∂ z ( ν ∗ , z ∗ ) = ∂ G k ∂ p ( ν ∗ , p ∗ ) (cid:44) , ∂ G k ∂ z ∂ν ( ν ∗ , z ∗ ) = ∂ G k ∂ p ∂ν ( ν ∗ , p ∗ ) + ∂ G k ∂ p ( ν ∗ , p ∗ ) p (cid:48) ( ν ∗ ) (cid:44) . (26)Using (24)–(26) we can write G k in the form G k ( ν, z ) = c ν z + c ν z + c z + O ( | z | ) , where c i , i = , , ν -dependent coe ffi cients with c (0) (cid:44) c (0) (cid:44)
0. Consider a ν -dependent substitution of the form z = w − h ( ν ) w . Then G k ( ν, w ) = c ν w + ( − c h + c ) ν w + ( c − c ν h ) w + O ( | w | ) . In order to annihilate the term ν w we should take h = c / c , which is well defined near ν ∗ = c (0) (cid:44)
0. With such choice of h , we have G k ( ν, w ) = c ν w + d w + O ( | w | ) , where d = c − c ν/ c and d (0) = c (0) (cid:44)
0. Put β = − d / c . From (16) and (26) we have c (0) c (0) <
0. Therefore β (0) = − c (0) / c (0) >
0. Setting α = c we conclude, that G k takesthe required form (22). (cid:3) Theorem 12.
Under assumptions of Theorem 8 (
C1–C4 and (14)), there is a smooth ν -dependentsubstitution p = p ( ν, w ) which brings the bifurcation function (8) to the normal formG k ( ν, w ) = α w ( ν − β w ) + O ( | w | ) , (27) where α = α ( ν ) , β = β ( ν ) and they are non-zero at ν = . Proof:
Let ν ∗ be the parameter value at which bifurcation occurs. Without loosing generality wemay assume that ν ∗ =
0. Reasoning as in the proof of Theorem 11, substitution z = p − p ( ν )brings G k to the form G k ( ν, z ) = zg k ( ν, z ). Put z ∗ =
0. From
C1–C4 and (14) we have that ∂ G k ∂ z ( ν ∗ , z ∗ ) = g k ( ν ∗ , z ∗ ) + z ∗ ∂ g k ∂ z ( ν ∗ , z ∗ ) = , ∂ G k ∂ z ( ν ∗ , z ∗ ) = ∂ G k ∂ p ( ν ∗ , p ∗ ) (cid:44) , ∂ G k ∂ z ∂ν ( ν ∗ , z ∗ ) = ∂ G k ∂ p ∂ν ( ν ∗ , p ∗ ) + ∂ G k ∂ p ( ν ∗ , p ∗ ) p (cid:48) ( ν ∗ ) (cid:44) . Hence, in these coordinates we have G k ( ν, z ) = c ν z + c z + O ( | z | ) , with c (0) (cid:44) c (0) (cid:44)
0. Setting α = c and β = c / c we obtain the required normal form(27). (cid:3) . Validation of bifurcations of symmetric periodic orbits: reversible Hamiltonian systems In this section we show, how to adopt the general framework described in Section 3 to au-tonomous Hamiltonian systems. We consider an R –reversible Hamiltonian system˙ x = J ∇ H ( x ) , (28)where J is the standard symplectic matrix and H : R n + → R is C -smooth. First, we choose a C -smooth Poincar´e section Π ⊂ R n + . If Fix( R ) ⊂ Π , by [48, Lemma 3.3] the Poincar´e map P : Π → Π is R | Π –reversible, too. In what follows we will study bifurcations of R -symmetricfixed points for P .In autonomous Hamiltonian systems, a natural choice of the bifurcation parameter is thevalue of H , because it is a constant of motion. On the other hand, from the point of view ofrigorous numerics, it is much easier and e ffi cient to work in the phase space coordinates andparametrize periodic orbits for P by one of them.This dissonance between numerical e ffi ciency and formal description of a bifurcation issolved in the following way. We expect that two families of periodic points for P intersect ata bifurcation point. This geometric condition can be checked in the phase-space coordinates.Additional conditions on H , which will be given in this section, will guarantee, that H can beused (locally) as a bifurcation parameter.Let us give a brief overview of the construction we are going to perform. Period-tuplingand touch-and-go bifurcations are local phenomena. Thus, we can choose local coordinates( p , p , q , q ) ∈ R × R n × R × R n near an apparent bifurcation point, such that1. Π = { ( p , p , q , q ) ∈ R × R n × R × R n : q = } and2. Fix( R ) = { ( p , p , q , q ) ∈ R × R n × R × R n : q = , q = } ⊂ Π .In order to simplify further notation we will use ( p , p , q ) coordinates in Π and we will alwaysskip q = P and H . We also assume, that the vector field (28) is transverseto Π near an apparent bifurcation point, so that the Poincar´e map is well defined and smooth.Define Π h = { ( p , p , q ) ∈ Π : H ( p , p , q ) = h } .Assuming that near an apparent bifurcation point there holds ∂ H ∂ p ( p , p , q ) (cid:44) Π h (cid:51) ( p , p , q ) → ( p , q ) is a local di ff eomorphism. Thisallows us to parametrize Π h locally by ( p ( h , p , q ) , p , q ) for h belonging to some interval H . Ourgoal is to specify conditions on P and H , which will guarantee, that the map f ( h , p , q ) = π ( p , q ) P ( p ( h , p , q ) , p , q ) (29)is well defined on some domain and it has period-tupling and / or touch-and-go bifurcations in thesense of Definition 3 and Definition 4, respectively.Let us fix k >
1. Following Section 3, we split p = ( p , p ) and q = ( q , q ). We assume,that the set of R –symmetric fixed points of P near an apparent bifurcation point forms a regularcurve, which is parametrized by the coordinate p P (cid:51) p → x Hf p ( p ) : = ( p , p ( p ) , p ( p ) , ∈ R × R × R n − × R n (30)15nd defined on an explicit, open interval P . We also assume, that for i = , . . . , k − P i , P × P × P × { } )) = { x Hf p ( p ) : p ∈ P } . (31)First we perform the Lyapunov-Schmidt reduction. We assume that there is a set P × P × P and a smooth function p H : P × P → P , such that for ( p , p , p ) ∈ P × P × P π q ( P k ( p , p , p , = ⇐⇒ p = p H ( p , p ) . (32)Using this implicit function we can define the bifurcation function by G Hk ( p , p ) = π q ( P k ( p , p , p H ( p , p ) , G Hk ( p , p ) = ( p − p ( p )) g Hk ( p , p ), where the reduced bifurcation functionreads g Hk ( p , p ) = (cid:90) ∂ G Hk ∂ p ( p , p ( p ) + t ( p − p ( p ))) dt . (34)Now we can specify the standing assumptions in the context of Hamiltonian systems: HC2: P is a closed interval, and P is a closed set, such that p ( p ) ∈ int P for p ∈ P , where p ( p ) is given by (30) and there is a smooth function p H : P × P → int P solving (32); HC3: there exists a smooth function P (cid:51) p → p ( p ) ∈ P such that (cid:110) ( p , p ) ∈ P × P : g Hk ( p , p ) = (cid:111) = { ( p ( p ) , p ) : p ∈ P } , where g Hk is defined by (34); HC4: there holds 0 (cid:60) ∂ G Hk ∂ p ∂ p ( P × P ) + ∂ G Hk ∂ p ( P × P ) p (cid:48) ( P ) . Assumptions
HC2–HC4 and Lemma 7 imply, that the set of R –symmetric fixed points for P k in P × P × P × { } is the union of two regular curves, which intersect at unique point x ∗ = ( p ∗ , p ∗ , p H ( p ∗ , p ∗ ) , p H is defined by (32). These parametric curves are defined onexplicit range p ∈ P and p ∈ P , respectively, which makes it possible for further continuationof these branches by the method described in Appendix A. On the other hand, we have noinformation about the type of bifurcation at the intersection point.Generically we expect, that the value of H can be used as the bifurcation parameter. It mayhappen, however, that H is not monotone (or even constant) along one or both curves of periodicpoints. In the remaining part of this section we derive (easily checkable by means of rigorousnumerics) conditions, which guarantee, that1. H can be used as the parameter near x ∗ and2. the mapping (29) undergoes one of the bifurcations defined by Definition 3 and / or Defini-tion 4.For further use we define four functions s ( p , p , q ) = ( H ( p , p , q ) , p , q ) , x Hb ( p ) = (cid:16) p ( p ) , p , p H ( p ( p ) , p ) , (cid:17) , h f p ( p ) = H ( x Hf p ( p )) and h b ( p ) = H ( x Hb ( p )) . emma 13. If P ( x ∗ ) = x ∗ and ∂ H ∂ p ( x ∗ ) (cid:44) , (35) then there is an open interval H and an open set U, such that s ( x ∗ ) ∈ H × U, s − | H× U is well defined di ff eomorphism onto image and P k ◦ s − is defined on H × U. Proof:
From (35) it follows, that det ( Ds ( x ∗ )) = ∂ H ∂ p ( x ∗ ) (cid:44) s is a local di ff eomorphismnear x ∗ . The inverse function takes the form s − ( h , p , q ) = ( p ( h , p , q ) , p , q ) . Since x ∗ is a fixed point of P , we can also find open sets H and U , such that s ( x ∗ ) ∈ H × U andsuch that P k ◦ s − is defined on H × U . (cid:3) Theorem 14.
Assume
HC2–HC4 and let x ∗ = ( p ∗ , p ∗ , p H ( p ∗ , p ∗ ) , be the unique intersectionpoint from Lemma 7. If k > is odd, ∂ H ∂ p ( x ∗ ) (cid:44) , h (cid:48) f p ( p ∗ ) (cid:44) and h (cid:48) b ( p ∗ ) (cid:44) , then f h defined by(29) has k th -order touch-and-go bifurcation at s ( x ∗ ) . Proof:
Take
H × U from Lemma 13. From (30)–(31) the point x Hb ( p ), p ∈ P is of principalperiod 2 k for P , except the intersection p ∗ . Since h (cid:48) f p ( p ∗ ) (cid:44) h (cid:48) b ( p ∗ ) (cid:44)
0, we can choose[ h , h ] ⊂ H , h ∗ ∈ ( h , h ), such that h − f p and h − b are defined on [ h , h ] and both functions x f p = π ( p , q ) ◦ x Hf p ◦ h − f p and x b = π ( p , q ) ◦ x Hb ◦ h − b have range in U . By the construction they satisfy conditions from Definition 4 for the function f defined by (29). (cid:3) Theorem 15.
Assume
HC2–HC4 and let x ∗ = ( p ∗ , p ∗ , p H ( p ∗ , p ∗ ) , be the unique intersectionpoint from Lemma 7. If k ≥ is even, ∂ H ∂ p ( x ∗ ) (cid:44) , h (cid:48) f p ( p ∗ ) > and h (cid:48)(cid:48) b ( p ∗ ) > , then f h definedby (29) has period k-tupling bifurcation at s ( x ∗ ) . Proof:
Take
H × U from Lemma 13. We will show that h (cid:48) b ( p ∗ ) =
0. Let us fix a sequence p n → p ∗ , such that p ∗ (cid:44) p n and x Hb ( p n ) is defined for all n ∈ N . Let us set x n = x Hb ( p n ) and y n = P k ( x n ). We also have y n (cid:44) x n , because by (31) the principal period of x n is 2 k . Thisand HC3 imply, that p n (cid:44) π p y n . Since k is even and x ∗ is a fixed point for P , we have thatlim n →∞ y n = x ∗ , y n (cid:44) x n and H ( y n ) = H ( x n ) = h b ( p n ). Hence, in the interval joining p n and π p y n there is a point at which h (cid:48) b vanishes and in consequence h (cid:48) b ( p ∗ ) =
0. The assumption h (cid:48)(cid:48) b ( p ∗ ) > p ∗ is a proper local minimum of h b .The remaining part of the construction goes as in Theorem 9. The function h b is convex near p ∗ ∈ int P which allows us to define two continuous branches h b i : [ h ∗ , h ] → int P , i = , , of h − b , for some h ∈ H . The function h f p is invertible near p ∗ , because we assumed that h (cid:48) f p ( p ∗ ) (cid:44)
0. Restricting the domain, if necessary, we may assume, that the inverse is defined on[ h , h ] with h ∗ ∈ ( h , h ). Define x f p = π ( p , q ) ◦ x Hf p ◦ h − f p and x b i = π ( p , q ) ◦ x Hb ◦ h b i , i = , . h , h ] if necessary, we may assume, that the range of the abovethree functions is in U . Thus f h defined by (29) has period k -tupling bifurcation at s ( x ∗ ). (cid:3) We end this section by a version of Theorem 10 for Hamiltonian systems. Proceeding as inSection 3.3 we assume, that (28) admits a symmetry S and a reversing symmetry R . We alsoassume that the function defined by (30) satisfiesFix( P , P × P × P × { } , C R ∧ C S ) = (cid:110) x Hf p ( p ) : p ∈ P (cid:111) . Then, we define G H and g H by (33) and (34), respectively. We state the following result withouta proof, as it is similar to that of Theorem 10 and Theorem 15. Theorem 16.
Assume
HC2–HC4 with k = and let x ∗ = ( p ∗ , p ∗ , p H ( p ∗ , p ∗ ) , be the uniqueintersection point from Lemma 7. If S and R commute, ∂ H ∂ p ( x ∗ ) (cid:44) , h (cid:48) f p ( p ∗ ) > and h (cid:48)(cid:48) b ( p ∗ ) > ,then f h defined by (29) has symmetry breaking bifurcation at s ( x ∗ ) .
5. Bifurcations of odd periodic solutions in the Falkner-Skan equation
The Falkner-Skan equation [35] is a third order ODE given by f (cid:48)(cid:48)(cid:48) + f f (cid:48)(cid:48) + β (cid:104) − ( f (cid:48) ) (cid:105) = . (36)It is low-dimensional and without singularities. Therefore it is relatively easy for rigorous nu-merical investigation. Although for some physical reasons solutions of certain BVP for (36) arerelevant, we use the system to test the methodology introduced in Section 3 and thus we will fo-cus on periodic solutions. Using this example we will also provide the reader with some detailsregarding validation of the presence of period-tupling and touch-and-go bifurcations. A higher-dimensional hamiltonian system (Circular Restricted Three Body Problem), which is also morecomputationally demanding, will be studied in Section 6.Our aim is to prove that some family of odd periodic solutions of (36) parametrized by β undergoes period-doubling, third order touch-and-go and period-quadrupling bifurcations.The equation (36) can be rewritten as a system of first order equations x (cid:48) = y , y (cid:48) = z , z (cid:48) = β ( y − − xz (37)and in the sequel we will work with (37). The system (37) is reversible with respect to R ( x , y , z ) = ( − x , y , − z ). For all β > , ± ,
0) are R -symmetric equilibria. Numerical simu-lation shows, that there is a family of R -symmetric periodic orbits u β ( t ) = ( x β ( t ) , y β ( t ) , z β ( t ))parametrized by β >
1. These periodic orbits intersect the symmetry line Fix( R ) = { (0 , y ,
0) : y ∈ R } at exactly two points, which approach (0 , ± , β → + — seeFig. 3. Thus, the period of these orbits goes to infinity when β → + . It is also observed, thatmax t ∈ R (cid:107) x β ( t ) (cid:107) goes to infinity when β → + .For small β >
1, these periodic orbits are hyperbolic. For larger β they become elliptic andcrossing strong 1 : k resonances, k = , , β k , respectively,where ˆ β = . , ˆ β = . , ˆ β = . . (38)18 log β y - - - - - - log β y Figure 3: The y coordinate of two intersection points of Fix( R ) with an approximate R -symmetric periodic solution of(37) as a function of β ∈ (1 ,
100 000].
Approximate initial conditions for resonant R -symmetric periodic orbits are (0 , ˆ y k , y = . , ˆ y = . , ˆ y = . . (39)and their shapes are shown in Fig. 4. It seems, the family u β continues to exists for all β > β → ∞ . - - x - - - - y - - t - - x Figure 4: R -symmetric periodic orbits corresponding to 1 : 2 (solid), 1 : 3 (dashed) and 1 : 4 (dotted) resonances. Let us define a Poincar´e section
Π = { ( x , y , z ) : z = } . We will use ( x , y ) coordinates to describe points in Π . By P β : Π → Π we denote the associatedPoincar´e map for the system (37) with the parameter value β . The map P β is reversible withrespect to involution R ( x , y ) = ( − x , y ). We will also use the notation P ( β, x , y ) = P β ( x , y ).The aim of this section is to give a computer-assisted proof of the following theorem.19 heorem 17. There is a smooth family u β = ( x β , y β ) of R-symmetric period-two points for P β parametrized by β ∈ J = [ ,
100 000] . This family undergoes period-doubling, third ordertouch-and-go and period-quadrupling bifurcations at some points ( β ∗ k , , y ∗ k ) , k = , , , respec-tively, with β ∗ k ∈ J k : = ˆ β k + | − , · − , y ∗ k ∈ Y k : = ˆ y k + [ − , · − , where approximate bifurcations points are listed in (38)–(39). Proof:
The existence of a smooth branch of R -symmetric period-two points for P β has beenproved by means of the method described in Appendix A. We used an adaptive cover of theparameter range J = (cid:83) j = J j , where the diameters of intervals J j are smaller (approximately2 · − ) if J j is close to and quite large (above 141) in the second end of the parameter range.Then, for each subinterval J j we check the assumptions of the parametrized Interval NewtonMethod (Lemma 25). If succeed, we prove that the segments y ( J j ) glue into a smooth curve bychecking the the bounds on y ( J j ) resulted from the Interval Newton operator overlap, when theirdomains do.We will give more details regarding validation of bifurcations. The bifurcation function is G k ( β, y ) = π x P k β (0 , y )for k = , ,
4. Hence, the Lyapunov-Schmidt reduction is not needed. In order to apply thegeneral framework introduced in Section 3 we need first to check conditions
C2–C4 . We havethe following bounds on G . G ( J k × { ˆ y k } ) ∂ G ∂ y ( J k × Y k ) − G ( J k × { ˆ y k } ) · ∂ G ∂ y ( J k × Y k ) − k = − . , . · − . ,
30] [ − . , . · − k = − . , . · − . ,
4] [ − . , . · − k = − . , . · − . ,
9] [ − . , . · − We see that in each case ˆ y k − G ( J k × { ˆ y k } ) · ∂ G ∂ y ( J k × Y k ) − ⊂ Y k , which proves that, there is abranch ( β, y ( β )) of zeroes of G parametrized by β ∈ J k , for k = , , C3 we need bounds on the reduced bifurcation function g k – see (10). Wehave g k ( { ˆ β k } × Y k ) ∂ g k ∂β ( J k × Y k ) − g k ( { ˆ β k } × Y k ) · ∂ g k ∂β ( J k × Y k ) − k = − , · − . , − [2 . , . · − k = − . , . · − . ,
31] [ − . , . · − k = − . , . · − . , − . , . · − Note, that in order to obtain tiny bounds on g k ( { ˆ β k }× Y k ) we used high precision interval arithmetic[36]. Again, in each case we have ˆ β k − g k ( { ˆ β k } × Y k ) · ∂ g k ∂β ( J k × Y k ) − ⊂ J k , which proves that for k = , , g k has branch of zeroes ( β ( y ) , y ) parametrized by y ∈ Y k and thus C3 issatisfied.From the following estimates ∂ G k ∂β∂ y ( J k × Y k ) ∂ G k ∂ y ( J k × Y k ) y (cid:48) ( J k ) k = . ,
9] [ − , · − . , · − k = . , − [27 . , . . , · − k = . , − . , . . , · − (40)20t follows that 0 (cid:60) ∂ G k ∂β∂ y ( J k × Y k ) + ∂ G k ∂ y ( J k × Y k ) y (cid:48) ( J k )for k = , , C4 is also satisfied.There remains to check conditions specific for each type of bifurcation. In (40) we havealready computed bound on ∂ G ∂ y ( J × Y ), which proves that the assumptions of Theorem 8are satisfied. Thus, the proof of the existence of third order touch-and-go bifurcation for P in J k × { } × Y k is completed.For period-doubling and period-quadrupling bifurcations we have to check non-degeneracycondition (16). From (40) we already have, that for k = , ∂ G k ∂β∂ y ( J k × Y k ) >
0. Thus,it su ffi ces to check that ∂ G k ∂ y ( J k × Y k ) is non-zero. We have the following bounds ∂ G ∂ y ( J × Y ) ⊂ [1372 , ∂ G ∂ y ( J × Y ) ⊂ − [1860 , . Eventually, we have to check that the principal periods of bifurcating orbits. The cases k = , k = g ( J × Y ) ⊂ . , . This completes the proof. (cid:3)
6. Halo orbits in the Circular Restricted Three Body Problem
In this section we apply the general framework described in Section 4 to the
Circular Re-stricted Three Body Problem . First, we will give a short overview of the CR3BP and we listsome of its relevant properties. Then, we will give a computer-assisted proof, that the wellknown families of halo orbits undergo various types of bifurcations.
The CR3BP is a mathematical model, that describes the motion of a small body with negli-gible mass under the gravitational influence of two point like big bodies, called primaries, whichrotate around their common mass centre on a circle.Denote by µ the relative mass of the smaller primary. In a rotating coordinate system centredat the common mass centre of two big primaries, the dynamics of the small particle is governedby the following system of second-order di ff erential equations [26, 44]¨ x − y = ∂ Ω µ ( x , y , z ) ∂ x , ¨ y + x = ∂ Ω µ ( x , y , z ) ∂ y , ¨ z = ∂ Ω µ ( x , y , z ) ∂ z , (41)where Ω µ ( x , y , z ) =
12 ( x + y ) + − µ (cid:112) ( x + µ ) + y + z + µ (cid:112) ( x − + µ ) + y + z . C µ ( x , y , z , ˙ x , ˙ y , ˙ z ) = Ω µ ( x , y , z ) − ( ˙ x + ˙ y + ˙ z ) . The hyperplane { ( x , y , z = , ˙ x , ˙ y , ˙ z = } is invariant under the local flow induced by (41) andthe corresponding four-dimensional Hamiltonian system is called the Planar Circular RestrictedThree Body Problem (PCR3BP).
The system possesses two main symmetries S : ( x ( t ) , y ( t ) , z ( t )) −→ ( x ( t ) , y ( t ) , − z ( t )) , R : ( x ( t ) , y ( t ) , z ( t )) −→ ( x ( − t ) , − y ( − t ) , z ( − t )) , (42)It is important to note that R and S commute. This property is required for our method ofvalidation of the existence of symmetry breaking bifurcations — see Theorem 18. Let us define the following Poincar´e section
Π = { ( x , y , z , ˙ x , ˙ y , ˙ z ) ∈ R : y = } and denote by P µ : Π → Π the associated Poincar´e map. We will skip the dependency on µ andwrite P , if it will be clear from the context. Since the y variable is fixed and equal to zero on thesection, we will use ( x , z , ˙ x , ˙ y , ˙ z ) coordinates to describe points in Π . With some abuse of notationon R , we will denote by the same letter the reversing symmetry of the system restricted to pointson Π , i.e. R ( x , z , ˙ x , ˙ y , ˙ z ) = ( x , z , − ˙ x , ˙ y , − ˙ z ) . Since Fix( R ) = { ( x , z , ˙ x = , ˙ y , ˙ z =
0) : x , z , ˙ y ∈ R } ⊂ Π by [48, Lemma 3.3] the mapping P is reversible, too. Thus, the frameworks for computer-assistedverification of various types of bifurcations introduced in Section 4 can be applied to P . The CR3BP possesses five equilibrium points, called the libration points . All of them arelocated in the { z = , ˙ z = } invariant hyperplane and thus they are equilibrium points for thePCR3BP, as well. Three of libration points, commonly denoted by L , L and L , are collinearand are located on the x -axis. They are of saddle-centre type for the PCR3BP. It is well known[44], that for all µ ∈ (0 ,
1) there exists a family of R –symmetric periodic orbits, called planarLyapunov orbits (PLO), that surround these libration points — see also Figure 5. In [11, 12]the existence of Lyapunov orbits around L and L libration points for selected mass parametershas been proved in an explicit domain. Computer-assisted methods have been used to prove[2, 45, 49, 50], that for the mass parameter µ = . L and L for certain energy level, and there is countableinfinity of connecting orbits between them in both directions.For the full system (CR3BP) the libration points L , , are of saddle-centre-centre type, forall µ . In additional direction z there exists a second family of vertical Lyapunov orbits (VLO),22hich are double symmetric both with respect to symmetry S and reversing symmetry R definedby (42). These orbits intersect twice the x -axis. Therefore, an object (a spacecraft) located nearbyone of those orbits will be periodically collinear with the two main primaries. Thus eclipses orshadows are unavoidable for trajectories approaching planar or vertical Lyapunov orbits.There is a numerical evidence [19, 21, 22], that a branch of out-of-plane R -symmetric or-bits, called halo orbits , bifurcates from the Lyapunov family. In opposite to planar and verticalLyapunov orbits, they never cross the x -axis, except at the bifurcation point. For large verticalamplitudes z the halo orbits become more and more aligned to the ( y , z ) plane allowing contin-uous observation of both primaries without eclipses. Parts of L -Lyapunov and L -halo familiesare shown in Figure 5 for the relative mass corresponding to the Sun-Jupiter system. Figure 5: A branch of planar Lyapunov orbits near L libration point and bifurcating halo orbits for the mass parameter µ S J = . · − corresponding to the Sun-Jupiter system — see also (43). The halo orbits in the CR3BP are R -symmetric, out-of-plane periodic orbits, which bifurcatefrom the planar Lyapunov family — see Figure 5. Although, there were extensive numericalstudy of these orbits and their bifurcations (just to mention few papers [22, 21, 19]), to the bestof our knowledge, they were never proved to exist.The best theoretical result in this direction has been done in [14, 15]. The authors considera certain normal form, which approximates the CR3BP. They prove, that in the normal form sys-tem there is a branch of R -symmetric halo orbits bifurcating from R -symmetric planar Lyapunovorbits. Moreover, they give an explicit expression for the bifurcation point. This result is validfor all L , , points and for all mass parameters.The results of this section are complementary to those from [14]. We give a computer-assistedproof of the existence of halo orbits for the original CR3BP system for all libration points L , , and for some (not full) range of µ . We will also study continuation and bifurcations of L , -halofamilies. 23 heorem 18. There are continuous functionsh i : [ µ ∗ , µ ∗ ] × Z (cid:51) ( µ, z ) → ( x i ( µ, z ) , , z , , ˙ y i ( µ, z ) , ∈ Fix( R ) with [ µ ∗ , µ ∗ ] = [9 . · − , . and Z = [ − , · ∆ z , where ∆ z = − , such that h i ( µ, is a point of symmetry breaking bifurcation of out-of-plane family of R–symmetricperiodic orbits from the planar family of Lyapunov orbits near L i libration point and for µ ∈ [ µ ∗ , µ ∗ ] and z ∈ Z the point h i ( µ, z ) is an initial condition for an R–symmetricperiodic out-of-plane (halo) orbit. Proof:
We applied the framework from Section 4 to the family of Poincar´e maps P µ . For fixed i = , , i to simplify notation). The rangeof parameters [ µ ∗ , µ ∗ ] is initially subdivided into N smaller overlapping subintervals [ µ ∗ , µ ∗ ] = (cid:83) Nj = µ j , where µ j = [ µ j , µ j ]. For each j = , . . . , N we execute in parallel the following (inde-pendent) tasks.1. We find an approximate bifurcation point ˆ u j = ( ˆ x , , , ˆ˙ y ,
0) for the mass parameter ˆ µ j = ( µ j + µ j ). For this purpose we use the scheme described in Appendix B — see alsoRemark 26.2. The planar double-symmetric Lyapunov orbits can be easily isolated by restriction to theplanar system. Using Lyapunov-Schmidt reduction and the method from Appendix A wevalidate the existence of smooth, two-parameter families of periodic orbits u f p : µ j × X → X × { } × { } × ˙ Y × { } , h i : µ j × Z → X × Z × { } × ˙ Y × { } , which correspond to Lyapunov orbits and halo orbits, respectively. The diameters of X and˙ Y were hand-optimized to speed-up computation.3. The set W = X ×{ }×{ }× ˙ Y ×{ } is a bound for the bifurcation point for each µ ∈ µ j . Thenwe check the assumptions of Theorem 16, i.e. ddx C µ j ( W ) (cid:44) ddx C µ j ( u f p ( µ j , X )) (cid:44) d dz C µ j h i ( µ j , z = (cid:44)
0. Note, that in Theorem 16 we required, that the branching-o ff curveis convex, but switching of either sign does not change the geometry of overall picture.If any of the above steps fails, the interval µ j is being subdivided and we repeat computation foreach element of subdivision. Finally, if all tasks are completed, using the methods from [20, 7]we check, that the pieces of h i glue into a smooth function. (cid:3) Graphs of µ → h i ( µ,
0) are shown in Figure 6. We would like to emphasize, that they matchquite well bifurcation points coming from the normal forms found in [14]. Our validation algo-rithm used to prove Theorem 18, by its construction, cannot continue with µ →
0. The thresholdvalue µ ∗ = . · − as well as size of out of plane amplitude ∆ z = − are by our choice acompromise between CPU time needed to obtain the result and the range of µ and z we can cover.In particular, the range of mass parameter [ µ ∗ , µ ∗ ] contains two relevant values µ S J = . · − , µ EM = . · − (43)corresponding to Sun-Jupiter and Earth-Moon systems, respectively. The values listed in (43)are taken as the nearest IEEE-754 double precision numbers to the recent mass measurementsreported in [40, 3]. 24 y C μ - L x y C μ - L x y C μ - L - - - - log ( μ ) e C L Figure 6: Curves of symmetry breaking bifurcation points ( x i ( µ ) , ˙ y i ( µ )) for L , , -Lyapunov families and the Jacobiconstant C µ ( x i ( µ ) , , , , ˙ y ( µ ) , µ has a local minimum. In order to make this minimum more evident, the plot is shown in log-exp scale. emark 19. The L is a special case, because for small µ the maximal order of normal formconstructed in [14] is , and it is divergent when µ → [42]. In this case we observe aninteresting phenomenon. If µ → , the Jacobi integral seems to be not monotone along thecurve of bifurcation points — see Figure 6 right-bottom panel. The computation, which indicatethe presence of a local minimum is non-rigorous but performed in high accuracy (400 bits ofmantissa) floating point arithmetic [36] and using th order ODE solver with the tolerance pertime step set to − .6.6. Continuation and bifurcations of halo orbits Theorem 18 guarantees that for µ ∈ [9 . · − , .
5] a branch of halo orbits near L i can beparametrized by out of plane amplitude z to at least | z | ≤ ∆ z . Numerical simulations [19, 21]strongly indicate, that these branches continue to exist for much larger amplitudes and that theyundergo period-doubling and period-quadrupling and third-order touch-and-go bifurcations. Thenext theorem addresses this issue. Theorem 20.
Consider the CR3BP with µ ∈ { µ S J , µ EM } as defined in (43). There exists a smoothfunction h µ : S → R such that for τ ∈ [0 , π ] the following holds true. h µ ( τ ) = ( x ( τ ) , , z ( τ ) , , ˙ y ( τ ) , is an initial condition for an R–symmetric periodic (halo)orbit. h µ ( τ ) = S ( h (2 π − τ )) – the family is S symmetric.This closed loop of halo orbits intersect the invariant subspace { z = , ˙ z = } at exactly twopoints h (0) and h ( π ) , at which an symmetry breaking bifurcation occurs. Moreover, the branchh µ undergoes period doubling, period quadrupling and third order touch-and-go bifurcations aslisted in Table 1. Proof:
The branch of halo orbits is split into four pieces.1. Proceeding as in the proof of Theorem 18 we validate the existence of two symmetrybreaking bifurcations: one in the region x > x < HC3 we also have, that there is a branch of out-of-plane R -symmetric periodic orbitsparametrized by z ∈ [ − ∆ z , ∆ z ], for an explicit ∆ z > z variable until hand-chosen threshold value z = . z = .
625 is parametrized by x variable.By the well known techniques [7] we can check, that these pieces glue into a smooth curve.By the symmetry we obtain the lower branch of halo orbits. Summarizing, we obtained, thatthe branch of halo orbits is a compact, smooth, one-dimensional manifold without boundary. Itis well known [33], that such a manifold is di ff eomorphic to a circle. By the construction themanifold is S -symmetric, thus the parametrization h can be chosen to preserve this symmetry, aswell.Approximate bifurcation points listed in Table1 (except for k = − ) using high-order ODE solvers from the CAPD library [10] based onhigh precision floating-point arithmetic[36]. Then we checked assumptions of Theorem 14 andTheorem 15 on very small sets (of the size about 10 − ) centred at these approximate bifurca-tion points. Notice, that in one case j = µ = µ S J we could not obtain bounds on thirdorder derivatives sharp enough to check convexity of Jacobi integral along bifurcation curve.We checked, however, the conditions
HC2–HC4 , which in particular means, that there are twocurves of halo orbits of period 1 and 4 intersecting at a single point. (cid:3) emark 21. Although we did not prove it, there is a strong numerical evidence that the pointsh µ (0) belongs to the family of L -Lyapunov orbits. A possible method to close this gap is to applythe method for computation of invariant manifolds of Lyapunov orbits, as proposed in [12]. Projections of the two curves h µ ( S ) for µ ∈ { µ S J , µ EM } resulting from Theorem 20 onto ( x , z )plane are shown in Figure 7. These families form a 2 D –tori in the full phase space. Such torusfor the mass µ = µ S J is shown in Figure 8. s z = - - - L Earth - Moon halo family on the Poincarè section z = - - - - L Sun - Jupiter halo family on the Poincarè section
Figure 7: Projection onto ( x , z )-plane of h µ ( S ), µ ∈ { µ S J , µ EM } as defined in Theorem 20. Each point h µ ( τ ) is an initialcondition for a halo orbit, thus the entire family forms a 2 D -tori in the full phase space – see Figure 8. Numerical simulation [22, 19] shows that L , -halo families continue to exists until a collisionwith one of the primaries — see Figure 9. Hence, on the Poincar´e section Π , ( x , z ) coordinatesof these orbits approach (1 − µ,
0) and ( µ,
0) respectively, while | ˙ y | tends to infinity.We have the following partial result for the L -halo family. Theorem 22.
Consider the CR3BP with µ ∈ { µ EM , µ S J } as defined in (43). There is a smoothfunction h µ : [ − , → R such that for τ ∈ [ − , the following statements hold true. h µ ( τ ) = ( x µ ( τ ) , , z µ ( τ ) , , ˙ y µ ( τ ) , is an initial condition for an R–symmetric periodic(halo) orbit. h µ ( τ ) = S ( h ( − τ )) – the family is S symmetric. h µ (0) is a point of an symmetry breaking bifurcation. The function ˙ y µ ( τ ) has a unique local minimum at τ = satisfying ˙ y µ EM (0) ∈ . , , ˙ y µ S J (0) ∈ . , . (44)5. The branches continue to at least ˙ y µ EM ( ± = . , ˙ y µ S J ( ± = . (45)27 igure 8: 2 D torus of primary L -halo orbits in the Sun-Jupiter system. Moreover, these branches undergo period doubling, period quadrupling and third order touch-and-go bifurcations as listed in Table 2.
Proof:
The validation is split into three steps.1. From Theorem 18 we know, that there is an symmetry breaking bifurcation of halo orbitsfrom L -Lyapunov family. The estimates (44) are taken from this proof. The branching offamily of halo orbits h µ ( z ) = ( x µ ( z ) , , z , , ˙ y µ ( z ) ,
0) is parametrized by z ∈ [ − , · ∆ z , with ∆ z = − . We also check that ˙ y (cid:48)(cid:48) µ ( z ) > | z | ≤ ∆ z .2. The branch is then rigorously continued using Appendix A and parametrized by h µ ( z ) = ( x µ ( z ) , , z , , ˙ y µ ( z ) , z > z (dependent on µ ). We also checked that for z ∈ [ ∆ z , ˆ z ] there holds ˙ y (cid:48) µ ( z ) > h µ (ˆ z ) is parametrized by ˙ y until hand-chosen threshold values (45).Summarizing, in each segment the variable ˙ y is increasing along the branch of periodic orbitswhich makes it possible to re-parametrize the curve as a function h µ ( τ ) = ( x µ ( τ ) , , z µ ( τ ) , , ˙ y µ ( τ ) , τ ∈ [0 , S we obtain the second branch for τ ∈ [ − , j =
0) were validated using Theorem 15 andTheorem 14 in high-precision interval arithmetics. Notice, that in two cases j = µ = { µ S J , µ EM } we could not obtain bounds on third order derivatives sharp enough to check convexityof Jacobi integral along bifurcation curve. We checked, however, the conditions HC2–HC4 ,28 igure 9: Intersection of families of L , , -halo orbits in the Earth-Moon system with the Poincar´e section Π along withstrong resonances, that lead to period-doubling and period-quadrupling bifurcations (marked by squares and disks in z > z < L , , -halo branches: L , -branches are shown on { y = , ˙ y > } section, L is shown on { y = , ˙ y < } section; (b) location of strong resonances on L -halo branch; (c) and (d) locationof strong resonances on L -halo branch. able 1: Bound on Jacobi constant for period k -tupling and touch-and-go bifurcations of L halo families for the massparameters µ EM and µ S J . Multiplicity k = k = j = µ = µ S J , denoted by a star, the convexity condition in Theorem 15has not been checked. We could not obtain bounds on third order derivatives so sharp, which would guarantee that thesecond derivative of Jacobi integral along bifurcation curve is non-zero. j k bound on Jacobi constant ( µ EM ) bound on Jacobi constant ( µ S J )0 1 3.17435[03, 36] 3.03588[11,51]1 4 3.08384097317512242430038839[79,91] 3.01979992774088569676042962[45,59]2 3 3.058886412529835176423178819[4,6] 3.015412945713342018918305256[4,5]3 2 3.0216192479264201986830801047[6,8] 3.00909434962110748263791018519[3,9]4 3 2.999986911642456326104353768[3,8] 3.00589961847845578773000478[64,75]5 3 2.997919501600216512922520[69,71] 3.00584514988954760426886596[07,61]6 3 2.94132864491556775199[28,32] 2.9941342902214648929134[15,26]7 4 2.940683922766931384[68,79] 2.9940756819941148370203478568[3,4] * (cid:3) Remark 23.
The threshold values in (45) have been chosen so that the validated arc of L –haloorbits contains all strong resonances at which period-tupling and touch-and-go bifurcations oc-curs for both values of µ ∈ { µ S J , µ EM } . Numerical simulation shows, that for larger values of ˙ y strong resonances are not present. However, we did not validate this conjecture. A possibleapproach to close this gap and obtain a full picture of what happens to L , –halo families is toperform Levi-Civita regularisation [44] and continue branches of halo orbits in these coordi-nates.6.7. Implementation notes The source code of all programs is available to download from the web page of the corre-sponding author [47]. The programs are written in C ++ -11 and use rigorous ODE solvers and al-gorithms for computation of Poincar´e maps and their derivatives [46, 53] from the CAPD library[10]. All programs the from the were compiled using g ++ -4.9.2 and executed on a computerequipped with Intel Xeon E7-8867 v4 2.40GHz processors (64 cores). Appendix A. Interval Newton method for implicit equations
In this section we provide a method for validated computation of implicit functions. Suchan algorithm is needed to check assumptions
C2–C3 or HC2–HC3 , as well as to compute widebranches of periodic orbits far from bifurcation points. The method is an adaptation of the wellknown
Interval Newton Method (INO) [34, 39] to the case of implicit equations. The mainmodification which significantly improves the method is the use of higher order derivatives. The30 able 2: Bound on Jacobi constant for period k -tupling and touch-and-go bifurcations of L halo families for the massparameters µ EM and µ S J . Multiplicity k = k = j = µ = { µ S J , µ EM } , denoted by a star, the convexity condition in Theorem15 has not been checked. We could not obtain bounds on third order derivatives sharp enough, which would guaranteethat the second derivative of Jacobi integral along bifurcation curve is non-zero. j k bound on Jacobi constant ( µ EM ) bound on Jacobi constant ( µ S J )0 1 3.15211[87,91] 3.03413[67,72]1 4 3.071869946146936057096946[46,52] 3.01887850935398942392012584[04,38]2 3 3.05083503865220863946099978[87,93] 3.014802161770970024871372364[53,86]3 2 3.0229911591596336379776897124[1,5] 3.0091557491590008844773187520[3,6]4 4 3.0158978159595140970471[78,93] n / a5 2 3.017143662542155479781194411[85,94] n / a6 4 3.0190035761795315910790[19,29] n / a7 3 3.077836508502735302[69,73] 3.019542146861187965925873562[20,98]8 4 3.104289978034786552471088092[8,9] * * method itself is quite straightforward but to the best of our knowledge, it has not appeared in theliterature. First we recall the Interval Newton Method [34, Thm. 8.4], [39, Thm. 5.1.7]. Theorem 24 ([34, 39]).
Let f : R n → R n be a C map and let X ⊂ R n be a convex, compact set.For x ∈ int X we define the interval Newton operator byN ( f , x , X ) = x − [ D f ( X )] − I f ( x ) , where by [ D f ( X )] I we mean a convex hull of the set of matrices { D f ( x ) : x ∈ X } . If N ( f , x , X ) ⊂ int X then f has a unique zero in X that belongs to N ( f , x , X ) . If N ( f , x , X ) ∩ X = ∅ then f has no zeros in X. The above theorem can be used to validate the existence of solutions to implicit equationsover an explicit domain.
Lemma 25.
Let f : R m × R n → R n be a C map, Z ⊂ R m be the closure of an open set and letX ⊂ R n be a convex, compact set with non-empty interior. For x ∈ int X we define the intervalNewton operator by N ( f , x , X , Z ) = x − [ D x f ( Z , X )] − I f ( Z , x ) . (A.1) If N ( f , x , X , Z ) ⊂ int X then there exists a C smooth function g : Z → X such that the set ofzeroes { ( z , x ) ∈ Z × X : f ( z , x ) = } coincides with the graph of the function g, i.e. { ( z , x ) ∈ Z × X : f ( z , x ) = } = { ( z , g ( z )) : z ∈ Z } . Proof:
Let us fix z ∈ Z and put f z = f ( z , · ). Then we have N ( f z , x , X ) ⊂ N ( f , x , X , Z ) ⊂ int X . From Theorem 24 for all z ∈ Z there exists a unique x = g ( z ) ∈ X such that f ( z , x ) =
0. Let usobserve, that the condition (A.1) implies, that D x f ( z , x ) is invertible at each point ( z , x = g ( z )).31y the implicit function theorem the function g is smooth, as it solves an implicit equation f ( z , g ( z )) ≡ z ∈ Z . (cid:3) The e ffi ciency of the INO (Theorem 24) is hidden in the fact, that we usually have a verygood approximation (from numerical experiments) for zero, i.e. f ( x ) ≈
0. Then the quantity[
D f ( X )] − I f ( x ) can be very tight, even if the computed bound on derivative D f ( X ) is overesti-mated. This is not the case for parametrized maps, as equation (A.1) contains the term f ( Z , x ).If Z is large then having f ( z , x ) ≈ z ∈ Z is rather unlikely.A straightforward way to overcome this problem is to make a substitution ( z , x ) = s ( z , w ),such that in the new coordinates the function z → w ( z ), which solves the implicit equation f ( s ( z , w )) =
0, is flat. Let us fix ( z , x ) ∈ Z × X , such that f ( z , x ) ≈ D x f ( z , x )is non-singular. We define an a ffi ne substitution by( z , x ) = s ( z , w ) : = ( z , x + w − A ( z − z )) , where A = D x f ( z , x ) − D z f ( z , x ). The map s is invertible because its linear part has de-terminant equal to one and thus zeroes of f are in one-to-one correspondence with zeroes of g : = f ◦ s . In the new coordinates the point ( z , w =
0) is an approximate zero of g . Moreover, D z g ( z , w ) = A is computed exactly). The above idea of changing linearly coordi-nate system has been proposed in [7], but the method for validation of a branch of zeroes wasdi ff erent and based on so-called radii polynomial approach [20, 7].In the remaining part of the section we will show, how to e ffi ciently evaluate all the terms,which appear in the INO for the mapping g . Our test show, that using this approach we couldsignificantly reduce overestimations in evaluation of the INO, which lead to significant advantageof this approach in comparison to direct evaluation of INO for g .The INO for the mapping g on the set Z × W , which contains ( z , w ) reads N ( g , w , Z , W ) = − [ D w g ( Z , W )] − I g ( Z , w ) . In what follows, we will show how we can bound all the terms that appear in this expression. Let X be such that s ( Z , W ) ⊂ ( Z , X ) and denote ∆ Z = Z − z . The term g ( Z , w ) can be bounded bymeans of the mean value theorem g ( Z , w ) ⊂ g ( z , w ) + [ D z g ( Z , w )] I · ∆ Z = f ( z , x ) + [ D z g ( Z , w )] I · ∆ Z (A.2)and the set of matrices D z g ( Z , w ) can be bounded by D z g ( Z , w ) ⊂ [ D z f ( Z , x )] − (cid:2) D x f ( Z , x ) A (cid:3) I ∩ (cid:104)(cid:16) D x f ( Z , x ) D x f ( z , x ) − (cid:17) D z f ( z , x ) (cid:105) . (A.3)By the choice of A , we have D z g ( z , w ) =
0. Therefore we expect that for not very large param-eter radius ∆ Z , the term D z g ( Z , w ) · ∆ Z in (A.2) is a small box around zero, as desired. Note, thatthe above considerations hold true for any matrix A . Therefore the quantities D x f ( z , x ) − and D z f ( z , x ) can be computed using just floating point arithmetic, however their product A mustbe bounded rigorously, if we want to take the intersection in (A.3).The quantity g ( Z , w ) can be also bounded using second order Taylor expansion g ( Z , w ) ⊂ f ( z , x ) + D z g ( z , w ) · ∆ Z + ∆ Z T D zz g ( Z , w ) ∆ Z . (A.4)32ecall, the point ( z , x ) is chosen so that f ( z , x ) ≈ s is chosen so that D z g ( z , w ) ≈ g ( Z , w ) is practi-cally quadratic in the radius of the parameter range ∆ Z . One can take the intersection of directevaluation of g ( Z , w ) in interval arithmetic with the bounds obtained from (A.2) and (A.4).In order to compute the Interval Newton Operator for g we have to bound D w g ( Z , W ). Directevaluation gives D w g ( Z , W ) ⊂ D x f ( Z , X ) . (A.5)Using second order derivatives of f we can obtain another enclosure D w g ( Z , W ) ⊂ D x f ( z , x ) + D ww g ( Z , W ) W + D zw g ( Z , W ) ∆ Z , which can be intersected with (A.5).Numerical experiments we performed show that using the above approach we can usuallyvalidate the existence of solution to the implicit equation on much wider domain Z (withoutsubdivision) than in the original coordinates.In principle, both g ( Z , w ) and D w g ( Z , W ) can be bounded using higher order Taylor expan-sions. This should come along with nonlinear (usually polynomial) substitution s , such that allderivatives of g = f ◦ s with respect to z vanish at ( z , w ) up to desired order r . Then the boundon g ( Z , w ) can be made of order O ( (cid:107) ∆ Z (cid:107) r + ).In the context of the CR3BP we have found, that the second order expansion is very e ffi -cient. Note, that computation of higher order derivatives of Poincar´e maps in a high dimensionalsystem is costly — the complexity of the C r –Lohner algorithm [53] used to integrate varia-tional equations is O ( n r n s ), where n is the dimension, s is the order of Taylor method and r ≥ r in ahigh-dimensional system is very expensive. Secondly, the bounds on higher order derivatives ofPoincar´e map are usually much overestimated than those of lower order. Appendix B. Newton like scheme for finding bifurcation points
A straightforward way to localize period-tupling and touch-and-go bifurcation points ofa family of reversible maps f ν ( x ) is to follow the branch of period-2 points x ( ν ) and look forresonant eigenvalues of D f ν ( x ( ν )). In low-dimensional systems one can look at the stability pa-rameter [21, 22]. This method is very e ffi cient when we want to find a rough approximation tothe bifurcation point. In reversible or hamiltonian case multiple eigenvalues may occur makingcomputation of eigenvalues with high-accuracy quite non-trivial task.In this short section we propose eigenvalue-independent yet e ffi cient scheme for finding veryaccurate approximation to bifurcation points. In what follows we focus on reversible Hamiltoniansystems and use the notation from Section 4, but the idea applies to any family of reversible maps.Following Section 4 we assume that ( p , p , p , q = ∈ Fix( R ) is an approximate period-2point for the Poincar´e map P , which is close to 1 : k resonance. We would like to refine it by aNewton-like scheme. Since we have n + p , p , p ) we need n + (cid:101) P ( p , p , p ) : = π q P ( p , p , p , = g Hk ( p , p ) = , (B.1)33here g Hk is defined by (34). The first equation selects the points from the curve of fixed points,while the second equation guarantees that the solution is also on the bifurcation curve. In orderto apply the Newton method to the system of equations (B.1), we need to compute derivatives of g Hk . From (34) we have ∂ g Hk ∂ p ( p , p ) = (cid:90) ∂ G Hk ∂ p ( p , p ( p ) + t ( p − p ( p ))) tdt . If the seed point for the Newton method is quite close to the bifurcation point, we may assumethat p ( p ) almost constant and thus ∂ g Hk ∂ p ( p , p ) ≈ ∂ G Hk ∂ p ( p , p ) (cid:90) tdt = ∂ G Hk ∂ p ( p , p ) . The second partial derivative reads ∂ g Hk ∂ p ( p , p ) = (cid:90) ∂ G Hk ∂ p ∂ p ( p , p ( p ) + t ( p − p ( p ))) dt + (cid:90) ∂ G Hk ∂ p ( p , p ( p ) + t ( p − p ( p ))) p (cid:48) ( p )(1 − t ) dt . Again, assuming p ≈ p ( p ) we can approximate ∂ g Hk ∂ p ( p , p ) ≈ ∂ G Hk ∂ p ∂ p ( p , p ) + ∂ G Hk ∂ p ( p , p ) p (cid:48) ( p ) . In order to compute D G Hk we need second order derivatives of P k and of the function p H ( p , p )obtained from the Lyapunov-Schmidt reduction — see (32). The latest can be computed bydi ff erentiation of the identity π q (cid:16) P k ( p , p , p H ( p , p ) , (cid:17) ≡ . Summarizing, the Newton-like iteration for equation (B.1) is given by( p m + , p m + , p m + ) = ( p m , p m , p m ) − r where r is the solution to the linear equation M · r = b with M = ∂ G Hk ( p m , p m ) ∂ p ∂ p + ∂ G Hk ( p m , p m ) ∂ p p (cid:48) ( p m ) ∂ G Hk ( p m , p m ) ∂ p ∂ (cid:101) P ( p m , p m , p m ) ∂ p ∂ (cid:101) P ( p m , p m , p m ) ∂ p ∂ (cid:101) P ( p m , p m , p m ) ∂ p , b = ∂ G Hk ( p m , p m ) ∂ p (cid:101) P ( p m , p m , p m ) . Using the above scheme, finding approximate bifurcation points of halo orbits with accuracy10 − was not a di ffi cult task. Remark 26.
In the computer-assisted proof of Theorem 18 we used similar strategy to localizeapproximate points of symmetry breaking bifurcations. We solved for zeroes of the function ( x , ˙ y ) → (cid:32) π ˙ x P ( x , , , , ˙ y , , ∂π ˙ z P ∂ z ( x , , , , ˙ y , (cid:33) . eferencesReferences [1] G. Alefeld, 1994. Inclusion methods for systems of nonlinear equations—the interval Newton method and mod-ifications. In: Topics in validated computations (Oldenburg, 1993). Vol. 5 of Stud. Comput. Math. North-Holland,Amsterdam, pp. 7–26.[2] G. Arioli, Periodic Orbits, Symbolic Dynamics and Topological Entropy for the Restricted 3-Body Problem , Comm.Math. Phys. 231 1 (2002) 1–24.[3] The Astronomical Almanac, http://asa.usno.navy.mil .[4] A. Abad, R. Barrio, F. Blesa and M. Rodr´ıguez,
Algorithm 924: TIDES, a Taylor series integrator for di ff erentialequations , ACM Transactions on Mathematical Software (TOMS) 39 (1), 5.[5] R. Barrio, F. Blesa and S. Serrano, Bifurcations and chaos in Hamiltonian systems , Int. J. Bif. Chaos, 20 (2010),1293–1319.[6] J. Burgos-Garcia, J.-P. Lessard and J.D. Mireles James,
Spatial periodic orbits in the equilateral circular restrictedfour body problem: computer-assisted proofs of existence , preprint.[7] J.B. van den Berg, J.-P. Lessard and K. Mischaikow,
Global smooth solution curves using rigorous branch following ,Mathematics of Computation, 79 (2010), 1565–1584.[8] M. Breden, J.-P. Lessard and M. Vanicat,
Global bifurcation diagram of steady states of systems of PDEs via rigorousnumerics: a 3-component reaction-di ff usion system , Acta Applicandae Mathematicae, 128 (2013), 113–152.[9] R. Barrio and D. Wilczak, Systematic Computer-Assisted Proof of branches of stable elliptic periodic orbits andsurrounding invariant tori , SIAM J. App. Dyn. Sys., 16 (2017), 1618–1649.[10] Computer assisted proofs in dynamics, a C ++ package for rigorous numerics, http://capd.ii.uj.edu.pl .[11] M.J. Capi´nski, Computer assisted existence proofs of Lyapunov orbits at L2 and transversal intersections of invari-ant manifolds in the Jupiter-Sun PCR3BP , SIAM J. Appl. Dyn. Syst. (2012), no. 4, 1723–1753.[12] M.J. Capi´nski and P. Roldan, Existence of a Center Manifold in a Practical Domain Around L1 in the RestrictedThree Body Problem , SIAM J. Appl. Dyn. Syst. , (2012) 285–318.[13] S.-N. Chow and J. Hale, Methods of Bifurcation Theory , Springer–Verlag, New York, 1982.[14] M. Ceccaroni, A. Celletti and G. Pucacco,
Halo orbits around the collinear points of the restricted three-bodyproblem , Physica D , (2016) 28–42.[15] A. Celletti, G. Pucacco and D. Stella,
Lissajous and Halo orbits in the restricted three-body problem , J. Non. Sci. , (2015) 343–370.[16] M.C. Ciocci and A. Vanderbauwhede, Bifurcation of periodic points in reversible di ff eomorphisms , Proceedings ofthe Sixth International Conference on Di ff erence Equations: new progress in di ff erence equations, (2004) 75–93.[17] A. Dhooge, W.J.F. Govaerts and Y.A. Kuznetsov, MATCONT: A MATLAB package for numerical bifurcation anal-ysis of ODEs , ACM Trans. Math. Soft., 29 (2003), pp. 141–164.[18] E.J. Doedel,
AUTO: Software for Continuation and Bifurcation Problems in Ordinary Di ff erential Equations , Tech.report, Applied Mathematics, California Institute of Technology, Pasadena, CA, 1986.[19] E.J. Doedel, V.A. Romanov, R.C. Pa ff enroth, H.B. Keller, D.J. Dichmann, J. Gal´an-Vioque and A. Vander-bauwhede, Elemental periodic orbits associated with the libration points in the Circular Restricted 3-Body Problem ,Int. J. Bif. Chaos, 17 (2007), 2625–2677.[20] M. Gameiro, J.-P. Lessard and K. Mischaikow,
Validated continuation over large parameter ranges for equilibriaof PDEs , Math. Comput. Simulation, 79 (2008), 1368–1382.[21] G. G´omez and J.M. Mondelo,
The dynamics around the collinear equilibrium points of the RTBP , Physica D ,(2001) 283–321.[22] K.C. Howell,
Three-dimensional, periodic, halo orbits , Celestial mechanics 32, (1984) 53–71.[23] `A. Jorba and M. Zou,
A software package for the numerical integration of ODE by means of high-order Taylormethods , Experimental Mathematics 14 (2005), 99–117.[24] T. Kapela,
N-BODY choreographies with a reflectional symmetry - computer assisted existence proofs , EQUADIFF2003, Proceedings of the International Conference on Di ff erential Equations, Hasselt, Belgium 22–26 July 2003, p.999-1005.[25] T. Kapela and P. Zgliczy´nski, The existence of simple choreographies for the N-body problem - a computer assistedproof , Nonlinearity 16 (2003) 1899–1918.[26] W.S. Koon, M.W. Lo, J. Marsden and S.D. Ross,
Dynamical Systems, the Three-body Problem and Space MissionDesign , Springer-Verlag New York, (2007).[27] Y.A. Kuznetsov and V.V. Levitin,
CONTENT: A Multiplatform Continuation Environment , Tech. report, CWI,Amsterdam, The Netherlands, 1996.[28] H. Kokubu, D. Wilczak and P. Zgliczy´nski,
Rigorous verification of cocoon bifurcations in the Michelson system ,Nonlinearity, (2007), 2147–2174.
29] J.S.W. Lamb, Reversing symmetries in dynamical systems, PhD Thesis, Universiteit van Amsterdam 1994.[30] J.-P. Lessard,
Rigorous verification of saddle-node bifurcations in ODEs , Indagationes Mathematicae 27 (2016),1013–1026.[31] J.-P. Lessard, E. Sander and T. Wanner,
Rigorous continuation of bifurcation points in the diblock copolymerequation , Journal of Computational Dynamics 4 (2017), 71–118.[32] D. Michelson, Steady solutions of the Kuramoto–Sivashinsky equation,
Physica D, , (1986), 89–111.[33] J. Milnor, Topology from the Di ff erentiable Viewpoint. , Princeton University Press, 1997.[34] R.E. Moore, Interval Analysis. , Prentice Hall, Englewood Cli ff s, N.J., 1966.[35] V.M. Falkner and S.W. Skan, Some approximate solutions of the boundary layer equations , Phil. Mag., (1931),865–896.[36] L. Fousse, G. Hanrot, V. Lef`evre, P. P´elissier and P. Zimmermann, MPFR: A Multiple-precision Binary Floating-point Library with Correct Rounding , ACM Trans. Math. Softw. (2007), Article No. 13.[37] F.J. Muoz-Almaraz, E. Freire, J. Galn, E. Doedel and A. Vanderbauwhede, Continuation of periodic orbits inconservative and Hamiltonian systems , Phys. D 181, (2003) 1–38.[38] M. Net and J. S´anchez,
Continuation of Bifurcations of Periodic Orbits for Large-Scale Systems , SIAM. J. App.Dyn. Sys., 14 (2015), 674–698.[39] A. Neumeier,
Interval methods for systems of equations.
Cambrigde University Press, 1990.[40] E.V. Pitjeva and E.M. Standish,
Proposals for the masses of the three largest asteroids, the Moon-Earth mass ratioand the Astronomical Unit , Celest. Mech. Dyn. Astr., 103 (2009), 365–372.[41] P. Pilarczyk,
Computer assisted method for proving existence of periodic orbits , Topol. Methods Nonlinear Anal.,Vol. 13 (1999), 365–377.[42] G. Pucacco, private communication.[43] O.E. R¨ossler,
An equation for continuous chaos , Phys. Lett. A 57 (5), (1976) 397–398.[44] V. Szebehely,
Theory of Orbits: The Restricted Problem of Three Bodies , Academic Press, 2013.[45] D. Sto ff er and U. Kirchgraber, Possible chaotic motion of comets in the Sun Jupiter system - an e ffi cient computer-assisted approach , Nonlinearity , (2004) 281–300.[46] I. Walawska and D. Wilczak, An implicit algorithm for validated enclosures of the solutions to variational equationsfor ODEs , App. Math. Comp., , (2016) 303–322.[47] D. Wilczak, personal homepage – a reference for auxiliary materials.[48] D. Wilczak,
Chaos in the Kuramoto–Sivashinsky equations – a computer-assisted proof , J. Di ff . Eq. , (2003)433–459.[49] D. Wilczak and P. Zgliczy´nski, Heteroclinic Connections between Periodic Orbits in Planar Restricted CircularThree Body Problem - A Computer Assisted Proof , Commun. Math. Phys. , (2003) 37–75.[50] D. Wilczak and P. Zgliczy´nski,
Heteroclinic Connections between Periodic Orbits in Planar Restricted CircularThree Body Problem. Part II. , Commun. Math. Phys. , (2005) 561–576 .[51] D. Wilczak and P. Zgliczy´nski,
Period Doubling in the R¨ossler SystemA Computer Assisted Proof , Found. Comp.Math. , (2009) 611–649.[52] D. Wilczak and P. Zgliczy´nski, Computer assisted proof of the existence of homoclinic tangency for the H´enon mapand for the forced-damped pendulum , SIAM J. Appl. Dyn. Syst. 8, Issue 4, pp. 1632–1663 (2009).[53] D. Wilczak and P. Zgliczy´nski, C r –Lohner algorithm , Schedae Informaticae, , (2011) 9–46.[54] C. Wul ff and A. Schebesch, Numerical Continuation of Symmetric Periodic Orbits , SIAM J. App. Dyn. Sys., 5(2006), 435–475.[55] P. Zgliczy´nski,
Steady states bifurcations for the Kuramoto-Sivashinsky equation - a computer assisted proof , J.Computational Dynamics, 2 (2015), 95–142., J.Computational Dynamics, 2 (2015), 95–142.