Computing (Un)stable Manifolds with Validated Error Bounds: Non-resonant and Resonant Spectra

We develop techniques for computing the (un)stable manifold at a hyperbolic equilibrium of an analytic vector field. Our approach is based on the so-called parametrization method for invariant manifolds. A feature of this approach is that it leads to a posteriori analysis of truncation errors which, when combined with careful management of round off errors, yields a mathematically rigorous enclosure of the manifold. The main novelty of the present work is that, by conjugating the dynamics on the manifold to a polynomial rather than a linear vector field, the computer-assisted analysis is successful even in the case when the eigenvalues fail to satisfy non-resonance conditions. This generically occurs in parametrized families of vector fields. As an example, we use the method as a crucial ingredient in a computational existence proof of a connecting orbit in an amplitude equation related to a pattern formation model that features eigenvalue resonances.


Mathematics Subject Classification
Stable and unstable manifolds are fundamental building blocks for understanding the global dynamics of nonlinear differential equations.Since closed-form analytic expressions for stable/unstable manifolds are rarely available, considerable effort goes into developing numerical techniques for their approximation, see e.g., Krauskopf et al. (2005), Haro and da la Llave (2006), Beyn and Kless (1998) and Castelli et al. (2015) and the references therein.One powerful tool for studying invariant manifolds (stable, unstable, strongly (un)stable and even more general ones) is the parameterization method of Cabré et al. (2003aCabré et al. ( , b, 2005)).The parameterization method is based on formulating certain operator equations (invariance equations) which simultaneously describe both the dynamics on the manifold and its embedding.The method has been implemented numerically to study a variety of problems involving stable and unstable manifolds of equilibria and fixed points (Mireles-James 2013;Mireles-James and Mischaikow 2013;Mireles James and Lomelí 2010;Haro 2011;van den Berg et al. 2011), stable/unstable manifolds of periodic orbits for differential equations (Guillamon and Huguet 2009;Huguet and de la Llave 2013;Castelli et al. 2015), quasiperiodic invariant sets in dynamical systems (Haro and da la Llave 2006) and more recently in order to simultaneously compute invariant manifolds with their unknown dynamics (Canadell and Haro 2014;Haro et al. 2014), to mention just a few examples.
The last reference is a book which provides many other examples and much fuller discussion of the literature.
In addition to facilitating efficient numerical computations, the functional analytic framework of the parameterization method also provides a natural setting for a posteriori analysis of errors.The works of van den Berg et al. (2011), Mireles-James and Mischaikow (2013) and Mireles-James (2015) exploit this a posteriori analysis and implement mathematically rigorous numerical validation methods for the stable and unstable manifolds.The term "validation" here expresses the fact that the computations provide explicit bounds on all approximation errors involved.
Approximate parametrizations are often computed by substituting a power series ansatz into the invariance equation and deriving a sequence of homological equations for the power series coefficients.These homological equations are then solved recursively to any desired order.In this paper we employ an alternative methodology for solving the invariance equation.We recast the infinite system of homological equations as a nonlinear zero finding problem on a Banach space of geometrically decaying sequences, and we implement a parametrized Newton-Kantorovich argument in the style of Yamamoto (1998).
There are three advantages to using a Newton method.First, we note the works of van den Berg et al. (2011), Mireles-James and Mischaikow (2013) and Mireles-James (2015) assume that certain non-resonance conditions between the eigenvalues are satisfied.The main goal of the present work is to weaken this assumption.We build on the theory of Cabré et al. (2003aCabré et al. ( , b, 2005) ) and develop computer-assisted methods for rigorous error bounding even in the face of a resonance.More precisely we develop validation schemes which apply at any co-dimension one resonance between the eigenvalues.The formulation of the resonant as well as non-resonant zero finding problems on an infinite sequence space unifies the presentation and implementation as well as the necessary a posteriori analysis.
Second, the zero finding methodology based on a numerical Newton method for the truncated problem leads to improved numerical performance even in the case where validated numerics are not desired.The Newton iteration can always be started from the linear approximation of the manifold by its eigenvectors; however, one has the option of improving the convergence by starting the iteration from a polynomial obtained by solving a few of the lower-order homological equations recursively.Once the iteration begins the order of the polynomial approximation is roughly doubled at each step.The freedom one has in choosing the lengths of the eigenvectors is exploited in order to guarantee that the Taylor coefficients decay at the desired rate.A similar numerical Newton scheme for computing invariant manifolds in discrete time dynamical systems (without rigorous validation) was used in Mireles James and Lomelí (2010).Indeed this kind of Newton scheme is always needed when the parameterization method is applied in Fourier space, see again Haro et al. (2014) and the references therein.
Third, existing rigorous continuation methods (van den Berg et al. 2010;Breden et al. 2013) are directly applicable to our novel approach, thus putting the rigorous computation of branches of connecting orbits within direct reach.Of course, from the viewpoint of continuation, all three advantages are crucial.Indeed, while continuing along one-parameter families of (un)stable manifolds resonances are encountered generically and are thus unavoidable.
Our use of parameterized Newton-Kantorovich arguments is guided by the work of a number of authors on the so-called method of radii polynomials.This method provides a tool kit for solution and continuation of zero finding problems in infinite dimensions.In this method one derives certain polynomials whose coefficients encode information about the approximate solution, the choice of approximate inverse of the derivative, the local regularity properties of the problem and the choice of Banach space on which to work.Once these givens are fixed the roots of the radii polynomials yield not only existence and uniqueness results, but also tight error estimates and isolation bounds for the problem.Radii polynomial methods have been applied successfully to a number of problems in dynamical systems, partial differential equations and delay equations, and we refer the interested reader to Day et al. (2007), Breden et al. (2013), van den Berg et al. (2010), van den Berg et al. (2011), Hungria et al. (2016) and van den Berg et al. (2015) for more discussion and references.
Let us now be more explicit about the above-mentioned resonance conditions (see Sects.2.2, 2.3 and2.4 andCabré et al. 2003a, b, 2005 for full details).We consider a stationary point of an analytic vector field u = g(u), with u ∈ R n .Focusing on the stable manifold, we denote the eigenvalues with negative real parts by λ 1 , . . ., λ d , where d is the dimension of the manifold.A resonance is a non-trivial relation for some ĩ ∈ {1, . . ., d} and ki ∈ N for i = 1, . . ., d.The trivial, excluded case is k = e ĩ , where we use the notation e i for the i-th unit vector.If there are no resonances between the stable eigenvalues (the "non-resonant" case), then there is an analytic conjugacy between the flow of u = g(u) on the local stable manifold and the linear flow θ i = λ i θ i , i = 1, . . ., d.Since the conjugacy map u = P(θ ) is analytic, P can be written as a convergent power series in θ .
When there is a resonance between the stable eigenvalues, the conjugacy map P as constructed above is no longer analytic and cannot be expressed as a convergent power series.One way to resolve this obstacle is to change the flow θ = h(θ ) in parameter space to a nonlinear "normal form," in such a way that the conjugacy is again analytic.In this paper we consider two types of resonances in particular, namely the co-dimension one resonances.The first type of resonance is a single "regular" resonance k and no other resonances.The second type is a double real eigenvalue ( k = e i for some i = ĩ) with geometric multiplicity one, and no other resonances.We note that a double eigenvalue is not a resonance in the strict sense, but we nevertheless use this "uniform" terminology in the current paper.These are the only co-dimension one resonances; hence, they are the types that are encountered generically in one-parameter continuation.For this reason we restrict our attention to these resonance types as our examples.We note that a completely analogous approach works for resonances of higher co-dimension, but here omit the details.
In order to illustrate the application of our methods we discuss three example problems in detail.First we analyze the stable manifold of the origin in the wellknown Lorenz equations.We use this model system to scrutinize our method in the non-resonant case and show how the structure of the vector field is directly reflected in the bounds used for validation.Second we tune the parameter in the Lorenz system to obtain double stable eigenvalues at the origin to showcase our method in this context.
In the final example the validated computation of stable and unstable manifolds is used as ingredient for the rigorous computation of connecting orbits.In particular we consider also the case of regular resonant eigenvalues.Specifically we compute connecting orbits in the system which arises as amplitude equations for the pattern formation model as shown in Doelman et al. (2003).Here U = U (t, x), t ≥ 0 and x ∈ R 2 .Equation ( 2) arises in the study of the interplay between trivial, hexagonal and roll patterns near the onset of instability of the zero solution.The parameters γ, μ and β are related via γ = μ β 2 .See Doelman et al. (2003) and van den Berg et al. (2015) for further details for the relation between ( 1) and (2).
Following the approach of van den Berg et al. (2015) we prove the existence of a heteroclinic connection between the hexagon and ground states.The proof is based on rigorous numerics for a boundary value problem, where the validated manifolds are In Fig. 1 we depict the verified connecting orbit of (1) as well as the corresponding stationary transition layer between hexagonal spots and the uniform state of (2).
We remark that an interesting future extension of this work would be to apply the ideas developed here to the validated computation of local stable/unstable manifolds of periodic orbits.The parameterization method has already been adapted to this context.See Guillamon and Huguet (2009), Huguet and de la Llave (2013) and Castelli et al. (2015) for more details and numerical examples.Using the techniques of the present work it should be possible to validate these computations even in the presence of a resonance between the Floquet multipliers.This will be the subject of a future study.
Finally we remark that the references mentioned in this introductory discussion are far from exhaustive, and a comprehensive overview of the literature is beyond the scope of the present work.In recent years a number of authors have developed numerical validation procedures which provide mathematically rigorous a posteriori error bounds on approximations of invariant manifolds associated to various kinds of invariant sets.We refer the interested reader to CAPD (2015), Capinksi and Simo (2012), Johnson and Tucker (2011), Wittig et al. (2010) and Wittig (2011) for fuller discussion of methods other than those presented here.
The outline of the paper is as follows.Section 2 is dedicated to the description of the general setup of our approach.In Sect. 3 we give more details on how we derive the zero finding problem.In Sect. 4 we transform it to an equivalent local fixed point problem to be solved by a parametrized Newton-Kantorovich-type argument.In Sect. 5 we illustrate the performance of our method with the three examples described above.The code implementing these examples can be found at the webpage Code page (2015).

Setup
We consider the validated computation of a parametrization of the local stable manifold of a hyperbolic fixed point p ∈ R n of a dynamical system induced by a nonlinear ODE using the parametrization method developed in Cabré et al. (2003aCabré et al. ( , b, 2005)).Local unstable manifolds can be obtained by replacing g with −g.We assume that g is analytic, allowing us to look for parametrizations in the analytic category.In particular we assume that g is locally (near p) analytically extendable to the complex plane.As a consequence the coefficients in the power series expansion of the parametrization decay geometrically, at an a priori unknown rate.We come back to the role of this decay rate later, see Remark 2.1 and in particular Sects.5.1.1 and 5.1.2.

The Invariance Equation
We denote by λ 1 , . . ., λ d the eigenvalues with negative real part of the Jacobian Dg( p) at the fixed point p.To fix notation, let there be s pairs of complex conjugate eigenvalues λ 1 , λ 2 , . . ., λ 2s−1 , λ 2s with negative real part and d − 2s real negative eigenvalues λ 2s+1 , . . ., λ d .The corresponding (generalized) eigenvectors are denoted by ξ 1 , . . ., ξ d .We do not assume that the algebraic multiplicity of the eigenvalues is one.For simplicity, in this paper we assume p, λ i and ξ i to be a priori determined analytically.However, it is straightforward to append equations for the equilibrium, as well as for the linearization around it, to the computational part of the analysis.We call the set of stable eigenvalues non-resonant if for every i = 1, . . ., for all k = (k 1 , . . ., k d ) ∈ N d with k = e i (the i-th unit vector).We shall use the notation As a new feature of the present work in comparison with the analysis in van den Berg et al. (2011), Mireles-James and Mischaikow (2013) and Mireles-James (2015) we are able to incorporate resonant cases directly in our novel framework for solving and validating the parametrization of the (un)stable manifold.In particular, we focus on the two types of co-dimension one resonances, namely a single regular resonance and an algebraically double, geometrically simple eigenvalue.We have a regular resonant case when for some ĩ ∈ {1, . . ., d} and a k ∈ N d with | k| ≥ 2. We work out in detail the cases where a regular resonance occurs as the only resonance, and where a double real eigenvalue occurs as the only resonance.In both cases, in the co-dimension one situation, the resonant eigenvalue is real valued.Higher co-dimension resonances, in particular multiple or simultaneous resonances (combinations of regular resonances and/or eigenvalues with higher multiplicity), can be dealt with analogously in our general framework.
In the parametrization method one looks for a map P that conjugates the flow of (3) on the stable manifold to a d-dimensional flow θ = h(θ ) for a suitably simple choice of h.Specifically, let P : C d ⊃ B ν → C n be analytic on the complex polydisc where ν = (ν 1 , . . ., ν d ) with ν i > 0. We commonly refer to its domain as the parameter space.In particular, the map P possesses a d-variate series expansion with a k = (a 1 k , . . ., a n k ) ∈ C n for |k| ≥ 0. We use the usual multi-index notation, where for k The invariance equation for P is given by Additionally and without loss of generality, we prescribe the linear constraints ) Remark 2.1 Note that the choice of the eigenvectors ξ i (i = 1, . . ., d) is not unique.In Lemmas 2.2, 2.3 and 2.5 we analyze the relation between the domain radius ν in (8) and the lengths ξ i .It will turn out that from a numerical perspective we can either vary ν or ξ i .See Remarks 5. 1, 5.2 and Breden et al. (2015), Falcolini and de la Llave (1992) and Cabré et al. (2003a) for a more thorough discussion of this topic.
Concerning the choice for the general polynomial normal form of h(θ ) in (8), note that in the non-resonant and double eigenvalue case it is sufficient to use information on the eigenvalues λ 1 , . . .λ d to do so.In the resonant case we choose a polynomial ansatz informed by spectral information but solve for the corresponding coeffient [τ in (27)].In all cases we shall choose it such that the origin is a (globally) attracting sink in parameter space.We will also see that we can find a subset B ν ⊂ B ν such that the orbits under the flow θ = h(θ ) of initial data in B ν do not leave B ν .We establish the explicit relation between ν and ν in the three cases under consideration in Lemma 2.7.By the conjugation property of P and (9a), a suitable real-valued restriction (see below) of the image of B ν under P thus gives us a parametrization of the local stable manifold W s ( p).We explain the conjugacy of the flows in Sect.2.5.
Let us describe how we can define a real-valued parametrization of the stable manifold starting with a complex valued solution of (8) fulfilling (9).Recall that λ = (λ 1 , . . ., λ 2s , λ 2s+1 . . .λ d ) consists of s pairs of complex conjugate eigenvalues and d − 2s real eigenvalues.We introduce the involutory permutation matrix σ s of the form and I d−2s denoting the d − 2s-dimensional identity matrix.In essence, this involution plays the role of a symmetry.We introduce an involution on d-tuples (including but not limited to C d , R d and N d ) by where v denotes complex conjugation.In particular, we have the involution k → k * on N d .We have ordered the eigenvalues so that λ * = λ, and normalized the (generalized) eigenvectors so that ξ * = ξ .Furthermore, for any variables q = (q k ) k∈N d that allow complex conjugation, we denote (again an involution) Next, consider the set for any compatible choice of ν ∈ R d + , i.e., ν * = ν.The set B sym ν is d-dimensional when interpreted as a real linear space (intersected with the ball B ν ), and in the particular case of all eigenvalues being real this boils down to Moreover, we will find, see Sect. 4, that the coefficients a of P have the symmetry property Under condition ( 14) we obtain the following lemma.
where the last equality uses the invariance of the summation domain under the involution * .
Together with Lemma 2.6 which explains the conjugation property of P in more detail, this establishes that P restricted to B sym ν parametrizes the real local stable manifold of p (see Lemma 2.7 for the relation between ν and ν).
Our goal is to compute a numerical approximation of P together with rigorous bounds on the approximation error and its range of validity B ν by using the method presented in Day et al. (2007).This amounts to first formulating an equivalent zero finding problem on an appropriate Banach space.Second, using an approximate zero, we define a Newton-like fixed point operator T .We establish contractivity of T on a ball around the approximate zero by deriving bounds on the residual, as well as bounds on the derivative that depend polynomially on the radius of the ball.These bounds are used to define so-called radii polynomials as ingredients for a finite set of inequalities encoding the prerequisites for the Banach fixed point theorem.We stress that in this way the radius of the ball on which we obtain contractivity is a variable for which we solve.This is an essential difference of the method in Day et al. (2007) compared to classical Newton-Kantorovich-type arguments.Let us assemble the ingredients to define the zero finding problem.

Non-resonant Eigenvalues
We distinguish two approaches for solving the invariance equation: the recursive approach and the zero finding approach.The zero finding approach will pave the way to an application of a fixed point argument in the space of power series coefficients.To be concrete, in order to explain the difference between the two approaches, we consider first the case that λ 1 , . . ., λ d are non-resonant.In this case we choose where s is the diagonal matrix containing the stable eigenvalues on the diagonal.
Clearly, the origin is the global attractor for θ = h(θ ).Moreover, By substituting the series expansion (7) into the invariance Eq. ( 8), we derive equations for the series coefficients This leads to the homological equations for all |k| ≥ 2: where b k only depends on a k with | k| < |k| and I n denotes the n-dimensional identity matrix.Note that b k vanishes for |k| = 1; hence, (17) reduces to the eigenvalueeigenvector equation, which is solved by (9b).
Using the initial constraints (9) for a k with |k| = 0, 1 and the fact that λ 1 , . . ., λ d are non-resonant, (17) can be used to compute a k recursively to any desired order (|k| ≤ N ).This is what we refer to as the recursive approach.The recursive approach shows that there is (a priori) a unique solution of (8) satisfying the constraints (9), although the decay of the sequence is not guaranteed a priori.The validation in van den Berg et al. (2011), Mireles-James and Mischaikow (2013) and Mireles-James (2015) relies on analysis in function spaces of so-called N -tails.
Equation ( 17) is derived by writing g(P(θ )) as a power series expansion in θ and matching like powers in the left-and right-hand sides of (8).In particular, where the notation k k means ki ≤ k i for i = 1, . . ., d.For later use, we also introduce the notation k ≺ k for those k k with k = k.Substituting (18) into the invariance Eq. ( 8) leads to the equations By observing that b k is defined by the splitting where b k depends only on a k with k ≺ k, one derives (17).In contrast, in the zero finding approach we omit the splitting (20) and instead interpret ( 19) as a zero finding problem on a space of geometrically decaying sequences {a k = (a 1 k , . . ., a n k )} k∈N d , see Sect. 3 and also Beyn and Kless (1998) for a related approach.
The recursive approach shows that the initial constraints (9) for a k with |k| = 0, 1 determine a k for |k| ≥ 2 uniquely.While a 0 = p is uniquely fixed by the problem at hand, we enjoy the freedom to scale a k for |k| = 1 as those are given by the eigenvectors of Dg( p).The following lemma illuminates how such a scaling influences a k for |k| ≥ 2. Let μ ∈ C d be given such that μ * = μ.We define the scaling μa def = (μ k a k ) k∈N d .In this way scaling preserves the symmetry ( 14).Moreover, for any analytic nonlinearity g, it follows from the power series representation (18) that We have the following invariance of the conjugation map under rescaling by μ.
We note that the above lemma is equivalent to the observation that the scaling θ i → μ i θ i leaves the flow in parameter space invariant, and hence by the conjugacy property, μa solves the homological equations whenever a does.In Remark 5.1 we come back to the practical implications of this scaling invariance.

Double Eigenvalues
We now consider a single repeated (real) eigenvalue with geometric multiplicity one (and no other resonances).Assume without loss of generality that λ 2s+1 = λ 2s+2 , with the rest of the eigenvalues being distinct and non-resonant.We choose ξ 2s+1 to be a (real valued) eigenvector for the double eigenvalue, and ξ 2s+2 a (real valued) generalized eigenvector such that Furthermore, we choose the flow in parameter space to be (cf.Jordan normal form) so that it is compatible with ( 23) through ( 8).Note that in the context of parameterdependent families of vector fields this approach can be numerically stabilized by using versal transformations, see for example Chow et al. (1994).In the current work no numerical instabilities occur however, as we do not consider continuation problems (although this is the subject of ongoing research).We observe, again, that the origin is the global attractor for θ = h(θ ), and that h(θ * ) = h(θ ) * .The corresponding version of ( 19) is for all |k| ≥ 2 with k 2s+2 ≥ 1.We note that the additional term (k 2s+1 + 1)a k+e 2s+1 −e 2s+2 occurs at the same "level of recursion" as the term a k , i.e., |k + e 2s+1 − e 2s+2 | = |k|.While this necessitates caution in the recursive computation of a k (which can still be done with the correct ordering), it does not introduce difficulties for our interpretation of (25) as a zero finding problem.Concerning the influence of scaling a k with μ ∈ C d (μ * = μ) we obtain the following.If one chooses a rescaling with μ 2s+1 = μ 2s+2 , the normal form (24) needs to be adapted accordingly.

Regular Resonant Eigenvalues
As discussed before, the two co-dimension one resonances are the (real) double eigenvalue (Sect.2.3) and the regular resonant eigenvalue (the other co-dimension one resonance) given by for some ĩ > 2s (i.e., λ ĩ ∈ R) and a k ∈ N d with | k| ≥ 2 and k * = k.Here we assume that no other resonances occur.
It follows from the discussion in Sect.2.2 that for the choice h(θ ) = s θ the mathematically equivalent recursive and the zero finding approaches are bound to fail.More precisely, assuming (26) to be fulfilled, it is apparent from (17) that the equation for a k is in general not solvable.To resolve this, we modify h to the nonlinear "normal form" where τ ∈ R is to be determined later.This alters (19) for k k (note that kĩ = 0) to Note that | k| ≥ 2 and thus |k + e ĩ − k| = |k| + 1 − | k| < |k|.This implies that ( 28) is amenable to recursive solving for k k, as in the non-resonant case.Obviously, for k = k Eq. ( 28) can be solved only if τ satisfies a solvability condition, which will be discussed below.We note that the choice (27) for h represents the simplest effective one.It facilitates an application of bordered matrix techniques (Govaerts 2000) to solve (28) in an efficient way while also providing a means to obtain uniqueness of a k , see below.The defining properties of the regular co-dimension one resonance imply that k * = k and e * ĩ = e ĩ .Hence, for any τ ∈ R we infer that h(θ * ) = h(θ ) * .However, since τ is a priori unknown, one needs to derive that τ is real indeed, see Eq. ( 31).
We now turn our attention to solving (28) for k = k.As we will see below, there is a unique τ making (28) for k = k solvable for a k , and by appending a suitable additional constraint we can make this solution unique.Using that kĩ = 0, we rephrase (28) for k = k as with b k depending only on a k with k ≺ k.Note that by construction a e ĩ ∈ ker(A), as it is an eigenvector of Dg( p) corresponding to the eigenvalue λ ĩ .The following lemma shows that appending an additional constraint on a k brings about unique solvability of (29) for the pair (τ, a k ).

Lemma 2.4
Let A ∈ R n,n have an algebraically simple eigenvalue 0 with eigenvec- is non-singular.Proof 2.3 The proof of this lemma can, for example, be found in Kuznetsov (2004), p. 174.
We choose a real vector ζ ∈ ker(A T ).From Lemma 2.4 we see that there is a unique pair (τ, a k ) satisfying This pair (τ, a k ) solves ( 29), and in addition Moreover, one derives that showing that the value of τ is independent of the choice of ζ .Furthermore, since b k is real for real a by its definition (20), τ is real whenever a is (ζ can be chosen to be real valued since λ ĩ is real).In this case we also obtain a scaling result.

Explicit Dynamics in Parameter Space
Let us explain how the invariance Eq. ( 8) encodes the conjugation of the flows of u = g(u) and θ = h(θ ) that we denote by (t, u) and (t, θ) for concreteness.In particular we explain when a restriction to a smaller ball B ν is in order.The following lemma contains the key observation and makes the role of P as a conjugation of flows precise.
In resonant cases we may need to restrict θ to a subset of B ν to ensure that P(θ ) ∈ W s ( p).For example, one may choose a smaller ball B ν ⊂ B ν to ensure (t, θ) ∈ B ν for all θ ∈ B ν and t ≥ 0. To be able to formulate a criterion for all three cases simultaneously, we introduce (34) The following lemma establishes both a uniform and a pointwise criterion.The former establishes a parametrization of W s loc ( p), whereas the latter is convenient when analyzing a specific orbit.Proof 2.6 For the non-resonant case this is a direct consequence of Lemma 2.6 and the fact that i (t, θ) ≤ θ i for all t ≥ 0. For the regular resonant case, the explicit flow for θ = h(θ ) is given by For i = ĩ we have | i (t, θ)| ≤ θ i for all t ≥ 0. For the resonant coordinate we infer This proves part (c) and part (a) for the regular resonant case follows from the inequality 0 (y 1 , y 2 ) ≤ 0 (|y 1 |, |y 2 |).Finally, the proof for the double eigenvalue case follows by putting ĩ = 2s + 1, k = e 2s+2 and τ = 1 in the above arguments.

The Zero Finding Problem
In this section we derive the zero finding problem on the space of geometrically decaying series coefficients whose solution corresponds to a solution P of ( 8) via ( 7).The functional analytic setup is close to the one utilized in Hungria et al. (2016) with the main difference lying in the convolution structure.

Spaces and Norms
As we work with analytical parametrizations P of the form (7), we consider the complex sequence spaces with ν ∈ R d + .Note that if a j = (a j k ) k∈N d ∈ W ν for j = 1, . . ., n, then P is well defined in B ν .If we define for w, w ∈ W ν the convolution operation the Banach space (W ν , • ν ) becomes a Banach algebra with multiplication * .In particular, Taking scalar constraints necessary in the resonant case into account we define the space where l 0 denotes the number of extra variables (as well as the number of constraints).
In particular l 0 = 0 in the non-resonant case and for double eigenvalues, whereas l 0 = 1 for (single) regular resonances.We have chosen this somewhat convoluted general notation to deal with all cases in one framework.The notation naturally allows incorporation of multiple (simultaneous) resonances by taking l 0 > 1.We denote elements x ∈ X ν l 0 by To avoid using projection operators we use the somewhat awkward looking, but compact, notation {x −l } l 0 l=1 for the scalar part of x.We endow X ν l 0 with the norm

Zero Finding Problem: Non-resonant Case
We define an operator f nonres on X ν 0 whose zeros correspond to analytic solutions P of (8) subject to the linear constraints (9).Based on (21), we set where x = a and Dg( p)ξ i = λ i ξ i .Given a vector field g and once ξ i is chosen, f nonres is fixed, and via (7) zeros of (39) correspond to a parametrization of the local stable manifold of p.

Zero Finding Problem: Double Eigenvalue
Based on (25), we define the operator f double on X ν 0 by where x = a, and ξ i are the (generalized) eigenvectors, as discussed in Sect.2.3.

Zero Finding Problem: Regular Resonant Case
In the regular resonant case, we choose a setup with τ as an extra unknown and (30) as appended equation.Based on (28) we define the operator f regres on X ν 1 by where x = (τ ; a) and ζ ∈ ker((Dg( p) − λ ĩ I n ) T ) is fixed (and real valued).We recall that we use negative indices to number scalar constraints; hence, we have slightly abused notation here by using k = −1 to denote the scalar part of f regres , whereas everywhere else k ∈ N d .

Fixed Point Operator and Radii Polynomials
We are now equipped with operators f nonres , f double and f regres given by (39) defined on X ν 0 , (40) defined on X ν 0 , and (41) defined on X ν 1 , respectively.We simply use the notation f defined on X if this does not lead to confusion.The zero of f corresponds to a parametrization of the stable manifold, a fact we still need to make more precise in due course.We follow the setup of Day et al. (2007) and Hungria et al. (2016) and derive an equivalent fixed point operator, whose contractivity on a ball around an approximate solution we establish using the so-called radii polynomials.Similar approaches have been used in previous works such as Eckmann et al. (1984), Arioli and Koch (2010) and Arioli and Koch (2015) and earlier in Lanford (1982).The structure of the fixed point operator T : X → X is Newton-like, i.e., where A, which is specified below, plays the role of an approximate inverse of D f ( x), with x an approximate solution to f (x) = 0. We extend the symmetry from Sect.2.1 to x ∈ X ν l 0 by setting We define T , and in particular A, in such a way that T allows a well-defined restriction to the closed symmetric subspace X sym given by where The function f defined in the three cases (39), ( 40) and ( 41) respects the symmetry. 123 Proof 4.1 We start by showing that c(a * ) = c(a) * .The vector field u = g(u) is real; hence g(u) = g(u).Considering the definition (18) of c(a), we will exploit the identity g(P(a, θ)) = g(P(a, θ)), where we write P(θ ) = P(a, θ) for clarity.We first observe, by using that the summation domain is invariant under the involution k → k * , that It follows that On the other hand By combining ( 45) and ( 46) we infer that It follows from uniqueness of the Taylor coefficients that c(a * ) = c(a) * .Using that (k • λ) = k * •λ this proves the assertion for the non-resonant case.Additionally, for the case of double eigenvalue and regular resonant cases we observe that the conditions k 2s+2 ≥ 1 and k k are invariant under the involution, since k * i = k i for i > 2s + 1 and k * = k.Finally, the symmetry of the scalar part f −1 in (41) follows from ζ = ζ and k * = k.
We now can make the correspondence between zeros of f and parametrizations of local stable manifolds more precise.
To define the linear operator A that appears in (42), we start by defining a finite dimensional projection π m on X .For m = (m 1 , . . ., m d ), we introduce the index set (47) Note that we require m * = m in the case of s complex conjugate eigenvalue pairs, so that k ∈ I m if and only if k * ∈ I m .Using this notation we set The range of π m can be identified with C M where M = l 0 + n d i=1 (m i + 1).The symmetry operation descends to an involution s on C M through An element x F in the range of π m can be lifted to X through the trivial embedding We define the "Galerkin" projection We now assume an approximate zero x F of f m has been computed (e.g., using a standard Newton method), and that A m is a numerical approximation of the inverse of the Jacobian, i.e., A m ≈ (D f m ( x F )) −1 .We denote by A s m def = 1 2 A m + s A m s its symmetrized version, which has the property Furthermore, replacing A m by its symmetrization A s m is not strictly necessary, since the symmetry of the fixed point, derived in Lemma 4.5 below, can also be obtained from the a priori uniqueness through Lemma 4.2.
We define the linear operator A by By construction, A conserves the symmetry, as expressed by the following lemma.We aim to show that T is contraction on a small ball around an approximate zero x = ι x F ∈ X .It follows from Lemma 4.5 that if x is symmetric ( x * = x) or almost symmetric, then the unique fixed point of T in the ball is a symmetric zero of f (provided A is injective).
Lemma 4.5 Assume that T is a contraction on the ball B x (r ) with B x (r )∩X sym = ∅.
Then T has a unique fixed point x in B x (r ), and x * = x.If, in addition, A s m is invertible, then x is a zero of f , and hence corresponds a parametrization of the real stable manifold.
Proof 4.5 The first part follows from the Banach fixed point theorem.Since T leaves X sym invariant, T is a contraction mapping on B x (r ) ∩ X sym = ∅, hence its fixed point lies in X sym .If A s m is injective, then the fixed point x corresponds to a zero of f , and the rest of the proof follows directly from Lemma 4.2.
As explained in Remark 4.1, the assertions of Lemma 4.5 also hold if one replaces A s m by its unsymmetrized analogue A m , provided it is invertible.We now describe how we show the contractivity of T on a suitable ball.Let us introduce the operator A † given by which acts as an approximation to the derivative of f at ι x F .Note that in particular A A † ≈ Id.For later use we note that DT (x)y can be split as Here Id denotes the identity on X .We continue by defining bounds that will be used to prove contractivity of T .
The crux of this definition is that the bounds Y on the residue and Z on the derivative of T can be constructed explicitly, see the examples in Sect. 5. We note that the inclusion of the scalar r in ( 55) trivially scales the bounds on DT by a factor r .We include this factor here to keep the notation compatible with earlier papers (going back to Yamamoto 1998).The radius r of the ball is not fixed a priori and the l 0 + n radii polynomials p(r ) are used in the following parametrized version of the Newton-Kantorovich theorem.
Lemma 4.6 Let r > 0 be such that p −l (r ) < 0 for l = 1, . . ., l 0 and p j (r ) < 0 for j = 1, . . ., n.Then T is a contraction on B x (r ) and there is a unique fixed point x ∈ B x (r ) of T .
Remark 4.2 We emphasize two consequences.
1.The l 0 + n conditions p(r ) < 0 reduce the validation of zeros of the operators f defined on the infinite dimensional spaces X to a finite set of inequalities that can be checked rigorously using interval arithmetic.Note that the inequality p(r ) < 0 is to be understood component-wise.2. We can translate p(r ) < 0 to a statement about the error of the image of the parametrization in phase space.Denote by P m (θ ) = k∈I m a k θ k the approximate parametrization corresponding to x, and by P(θ ) = k∈N d a k θ k the exact parametrization corresponding to x.Then for all θ ∈ B sym ν and all j = 1, . . ., n This gives control over the tail coefficients in the exact parametrization.Note that it is this information that is crucial to the method in van den Berg et al. (2011) for deriving a posteriori bounds on the derivative of the truncation error in the parametrization (see van den Berg et al. (2011), Section 5.2, Eqn. ( 87)).As a consequence, the analysis in van den Berg et al. ( 2011) is applicable with the radius r obtained from the current approach.
To derive the bounds ( 54) and ( 55) good analytical control on the operator f and its derivative D f at x is essential.In Sects.5.1 and 5.2 we illustrate the "mechanics" involved in the derivation of these bounds.

Applications
In this section we consider three applications to illustrate the performance of our method.We start with the well-known Lorenz equations and compute a 2D local stable manifold at the classical parameter values yielding non-resonant eigenvalues.Subsequently, we consider non-standard parameter values in the Lorenz system to investigate the specific issues in the double eigenvalue case.Finally we consider the case of regular resonant eigenvalues in a system of ODEs originating from a pattern formation model.We compute (un)stable manifolds that serve as ingredient for a connecting orbit computation in this system.

Local Manifolds in the Lorenz System
We consider the Lorenz differential equations (Lorenz 1963) In Sect.5.1.1 we choose the classical parameter values β = 8 3 and σ = 10 and ρ = 28 and describe the validation process including the derivation of the necessary bounds for the 2D local stable manifold of the origin (local Lorenz manifold) in some detail.In particular, we investigate how to choose the various computational parameters involved in the validation process taking different objectives into account.In Sect.5.1.2we use the fact that we have explicit formulas for the stable eigenvalues at the origin to tune the parameters β, ρ, σ such that there are double eigenvalues at the origin.We explain the validation analysis in this context, which serves as preparation for the regular resonant case.

Detailed Analysis for the Local Lorenz Manifold: Non-resonant Eigenvalues
The dimension of the manifold is d = 2, and the phase space dimension is n = 3. Denote the stable eigenvalues by λ 1,2 and the corresponding eigenvectors by ξ 1,2 .As a matter of fact λ 1,2 ∈ R. Hence, we search for a map P : R 2 ⊃ B ν → R 3 together with ν = (ν 1 , ν 2 ) fulfilling ( 8), with g given by the Lorenz vector field (59), and the linear constraints ( 9).As we are in the non-resonant case we can choose ν = ν as discussed in Sect.2.5.We make the power series ansatz (7) with a k = (a with * defined in ( 36).Together with p = 0 ∈ R 3 this completes the ingredients for f nonres defined in (39).
Recalling the definition of the radii polynomials in (56b), we notice that a necessary condition for finding a radius r fulfilling p i (r ) < 0 for i = 1, 2, 3 is that the components of Z 1 derived in (65) be smaller than one.The parameters that are under our direct control are m = (m 1 , m 2 ) and ν = (ν 1 , ν 2 ) or ξ 1,2 , respectively.Note that varying ν 1,2 is in the following precise sense equivalent to varying ξ 1,2 , which implies that we may as well fix either ν or ξ .
Remark 5.1 From the scaling operation analyzed in Lemma 2.2, we notice that This implies that we should either vary the decay rates (domain radius) ν or the (eigenvector) scalings μ to affect â j ν , j = 1, 2, 3 in Z 1 .Below we fix ν = (1, 1), as this is numerically the most stable choice.
Let us describe the dependence of the coefficients of p i (r ) on these computational parameters in more detail by deriving explicit formulas for the radii polynomials.

Derivation of the bounds: details
To define the radii polynomials specified in (56b) we first need to compute the bounds Y j and Z j ( j = 1, 2, 3) defined in ( 54) and ( 55).Assume we have already calculated an approximate zero x = ( â1 , â2 , â3 ) = ι xF where ι is defined in (49), and set A m ≈ (D f m ( xF )) −1 .We start by noting that ( f ( x)) k = 0 for k / ∈ I 2m , since g(u) is quadratic.Using this and the fact that The bounds Y j ( j = 1, 2, 3) are then obtained by computing the finite sums y j ν .
To derive the bounds Z j (r ) defined in (55) we use the splitting (53).Let v, w ∈ B 1 (0).We start by deriving an expansion The explicit expressions for z are listed in Table 1.
Our next goal is to compute a bound The first term to estimate is and by using |w where absolute values are to be understood component-wise.Furthermore, for the tail terms we have To accommodate matrix multiplication we collect all estimates of z k,2 for k ∈ I m in a vector.Namely, we define the vector To be explicit, we use the estimate Then applying the definition of A given in (51), and interpreting A m as an 3n m × 3n m matrix, yields the values for Z k,1 and Z k,2 summarized in Table 2.
All absolute values are to be understood component-wise.The values of k are defined in (63) Finally, by estimating the finite part (k ∈ I m ) and the tail part (k / ∈ I m ) separately, we can compute where we have used the Banach algebra convolution estimate in W ν to bound the tail sums.This completes the ingredients for (56b).
The main influence of ν on these estimates is via the terms âi ν (i = 1, 2, 3).On the other hand these norms are also controlled by the length of the stable eigenvectors ξ 1,2 appearing in definition (39), and also, albeit weakly, by m = (m 1 , m 2 ).Let us analyze this interplay.
Choice of m and ξ 1,2 The considerations that lead to the settings of the computational parameters depend on the goals in the application at hand.One might, for example, be interested in uniformly maximizing the image in phase space of the parametrization, or one might want to capture a particular point in phase space in the image, i.e., maximize the image in a certain direction.In the following we collect several observations for either task.Let us begin with a fundamental restriction; namely, the components Z j 1 in (65) need to be smaller than one.This can be used as a first feasibility check for the validation.
This motivates the following procedure (recall that we have fixed ν = (1, 1), see Remark 5.1):  (15,15)) We see that for a larger number of modes we obtain smaller error bounds.Note in addition that the larger the norm of ξ 1,2 is the bigger the uniform error bound r on B ν gets.Bottom dependence of the norm of the approximate solution on the number of rescalings, hence on the norms ξ 1,2 .These are an indicator for the size of the image in phase space We consider the two cases m = (5, 5) and m = (15, 15), see Fig. 2. First, for m = (5, 5), the validation succeeds with ξ 1,2 = 1 with a radius of r ≈ 10 −7 .Recall that by Remark 4.2 the validation radius r can be seen as an accuracy measure.We thus consider smaller radius r as higher accuracy.After conducting 5 rescalings with factors μ 1 = μ 2 = 7 6 to uniformly maximize the image, we fail to validate.The accuracy r we obtain after 4 rescalings decreased to 10 −5 (for fixed m the (uniform estimate on the) accuracy naturally decreases when increasing the domain).Second, for m = (15, 15) we also succeed to validate for ξ 1,2 = 1 but with smaller uniform error bound r ≈ 10 −15 .We are able to rescale 15 times with the same factors.This increases ξ 1,2 to 10.09 and increases the validation radius to 10 −2 .2. Fast-slow choice of m and μ: maximizing the image in the slow direction Next we consider the cases m 1 = m 2 and/or μ 1 = μ 2 , see Fig. 3.We recall that |λ 1 | ≈ 10|λ 2 |, hence we refer to λ 1/2 as the fast/slow eigenvalue.For most orbits the dynamics close to the origin is dominated by the slow direction.This ≈ 11.8.Bottom dependence of the norm of the approximate solution on the number of rescalings, hence on the norms ξ 1,2 .Note the clear dominance of the â3 ν which reflects the fact that ξ 2 = (0, 0, 1) T is, for example, of interest when computing connecting orbits that approach the equilibrium along the slow direction.Capturing a large portion of the slow direction can thus be desirable.We choose m 2 ≥ m 1 and μ 2 > 1 > μ 1 .First we choose m 1 = m 2 = 15 .We succeed to validate for ξ 1,2 = 1 with a radius r ≈ 10 −15 .Let μ 1 = 6 7 and μ 2 = 7 6 .We obtain 14 successful rescalings with gradually decreasing accuracy (r ≈ 10 −3 after 14 rescalings).If we choose m 1 = 5 and m 2 = 15 we observe qualitatively different behavior of the validation radii.Starting with a success at ξ 1,2 = 1 and radius r ≈ 10 −8 the accuracy increases to 10 −11 in the first 8 rescalings with the factors μ 1 = 6 7 and μ 2 = 7 6 until the norms of ξ 1,2 "align" with the choice of m.Then the accuracy decreases to 10 −3 after 18 rescalings.
The above considerations can serve as a starting point for more elaborate future investigations.One might, for example, devise an optimization scheme in which one takes not only the radius r as an unknown in the radii polynomials but also considers ν or ξ 1,2 , respectively, as variables.
Table 3 Expansion coefficients for (D f ( x + r v)r w − A † r w) k for the double eigenvalue case, using the standard Kronecker δ symbol

Analysis for the Local Lorenz Manifold: Double Eigenvalues
In order to analyze the situation for double eigenvalues as discussed in Sect.2.3, we choose for the parameters in the Lorenz system (59) the relation ρ = 1+(σ +1) 2 /(4σ ), leading to double eigenvalues λ = −(σ + 1)/2.Using this data we set up the operator f double specified in (40).Note that the (generalized) eigenvectors ξ 1,2 fulfilling ( 23) can be computed explicitly in this model case.For the general case we refer to methods developed in Alefeld and Spreuer (1986) and Rump (2001).It will be subject of future work to discuss their applicability in the current context.Let us now discuss the influence of the choice of parameters m and ν (or ξ 1,2 ) we first define the radii polynomials as we did in 5.1.1.
To compute the bounds Y j ( j = 1, 2, 3) defined in ( 54) and Z j ( j = 1, 2, 3) defined in (55) we follow the same strategy as in Sect.5.1.1.First the derivation of the Y -bounds is exactly analogous.In deriving z k,i (r ), i = 1, 2 fulfilling in the analogue of (61) we notice that the only difference from the formulas in Table 1 is induced by the additional linear term (k 1 + 1)a k , with k def = (k 1 + 1, k 2 − 1) in f double for |k| ≥ 2 and k 2 > 0. For k 2 = 0 one should read a k ≡ 0. The resulting differences in z k,1 can be read off in Table 3.
To obtain Z k (r ) fulfilling the equivalent of (62) we obtain, in addition to k from (63) and χ 2 from (64), a vector χ 1 ∈ R 3n m such that To be precise, we set χ for all 1 ≤ k 2 ≤ m 2 and j = 1, 2, 3, whereas all other components of χ 1 vanish.To obtain the analogue of Z 1 in (65) we need to bound the sum The following lemma can be used to control the factor k 1 +1 |λ 1 |(k 1 +k 2 ) uniformly for k / ∈ I m .We formulate it in this more general form as it will be reused for more general terms of this type in analyzing resonant eigenvalues in Sect.5.2.
The absolute value |A m | is to be understood component-wise Lemma 5.1 Let m ∈ N d with m ĩ ≥ 1 for an ĩ ∈ {1, . . ., d}.Let Then we have, with Proof 5.1 For any k ∈ I m there is at least one j ∈ {1, . . ., d} such that k j ≥ m j + 1.We distinguish two cases: • k ĩ ≥ m ĩ + 1.Then .
• 0 ≤ k ĩ ≤ m ĩ .We estimate The proof follows from combining these estimates.
Following the same steps as in Sect.5.1.1 we obtain Z k,1 as given in Table 4, while Z k,2 remain unaltered from Table 2.
Hence, one finds (in all the cases discussed below) that Z 1 is given by .
with u = (u 1 , u 2 , u 3 , u 4 ) T .The relation between the parameters γ, μ and β is given by γ = μ β 2 .We refer to Doelman et al. (2003) for a complete description of the seminal asymptotic reduction of the PDE (2) to the ODE system (67), obtained via a combination of spatial dynamics, bifurcation theory and geometric singular perturbation theory.For further details about the particular form (67) one may consult (van den Berg et al. 2015).We take the same viewpoint as in van den Berg et al. (2015) and investigate by rigorous numerical techniques, for fixed parameter value γ , connecting orbits between equilibria of (67).While in van den Berg et al. (2015) the configuration for the coexistence of hexagons and rolls (spots and stripes) is considered, we concern ourselves with the coexistence of the trivial state and hexagonal patterns.In terms of (67) this corresponds to a connecting orbit between the two fixed points , and p 2 = (0, 0, 0, 0) T .As ( 67) is Hamiltonian, a necessary condition for the connection between p 1 and p 2 to exist is for the equilibria to lie on the same energy level, which is the case for γ = − 2 135 .The equilibrium p 2 has resonant eigenvalues (λ 1 = 2λ 2 ).For this reason the case of coexistence of the trivial state and the hexagonal spot pattern was not considered in van den Berg et al. (2015).However, our current approach to validating the (un)stable manifold is well suited for this situation.
To compute the connecting orbit we adapt the approach of van den Berg et al. (2015).Introducing a rescaling factor L > 0 we aim at solving where θ(ψ) def = (ρ cos(ψ), ρ sin(ψ)) with some fixed ρ < 1 (playing the role of the phase condition in this otherwise autonomous problem).The fact that the second component u 2 (1) − P s 2 (θ (ψ)) = 0 can be replaced by the a posteriori check of sign(u 2 (1)) = sign(P s 2 (θ (ψ))) is explained in Lemma 2 in van den Berg et al. (2015) and is related to the Hamiltonian nature of the problem.Using the phase condition (i.e., fixing ρ) and omitting the second component deals on the one hand with the fact that every time shift of a connecting orbit is again a connecting orbit and on the other hand with the fact that the intersection of the two-dimensional unstable and stable manifolds corresponding to the connecting orbit is not transverse in R 4 .In a nutshell, they guarantee the isolation of the zero of F that we set out to find.
To rigorously compute zeros of (69) we choose to discretize the time dependence using a Chebyshev series (other choices, such as splines, are also possible).For details on how this is done in this example we refer the reader to van den Berg et al. (2015) and for the general method to Lessard and Reinhardt (2014).In this paper we focus on the rigorous computation of the maps P s,u and especially P s , as there we encounter eigenvalue resonances.The validated computation of P u is conducted analogously to the Lorenz equation explained above and we do not give any further details below.Furthermore, in order to validate zeros of (69) we need rigorous information on D P s,u (θ ) that we obtain in the same way as explained in Remark 3 in van den Berg et al. (2015).Note that by Remark 4.2 the a posteriori bound δ s,u corresponds to our validation radius r .
We now delve into the validated computation of the stable manifold of the origin, which has resonant eigenvalues.Explicitly, the eigenvalues of Dg( p 2 ) with negative real part are given by λ We note that 2λ 2 = λ 1 , so condition (5) holds with ĩ = 1 and k = (0, 2).Following the approach described in Sect.2.4 we choose the nonlinear normal form to describe the dynamics in the (stable) parameter space.The coefficients c k = c k (a; γ ) in (41) are given by This completes the ingredients for (41) in the particular case of (67).Note that the condition in Lemma 2.7c) reads 0 θ 1 , fulfilling the analogue of (61).The result is summarized in Table 6 in "Appendix 1."Note that z −1 (r ) = 0, as f −1 is linear in a k and k m for our choice of m.
The main deviation from the non-resonant case is given by the off-diagonal linear terms introduced by the resonance.Note that this is analogous to the double eigenvalue case, but with the additional difference given by the presence of the additional unknown τ .To expand the bounds Note that χ 1 i = 0 for i = 1, 2, 3, since f −1 is linear.Then, applying the definition of A given in (51) yields the values for Z k,1 , Z k,2 and Z k,3 summarized in Table 7 in "Appendix 1." Summing up Z k,i component-wise while splitting into a finite part and infinite tail yields Z (r ).The result is shown in Table 8 in "Appendix 1."An important ingredient in computing Z 1 is the control of terms of the form where We use Lemma 5.1 to bound the factor In the specific problem under consideration we are in luck, since τ ≈ 0, hence the contribution from the term (73) is small.This stems from the fact that the exact solution is τ = 0 as we derive analytically in "Appendix 2." However, the strength of our method is that one does not need to determine the coefficients of the normal form beforehand, as they are part of the overall set of unknowns for the nonlinear problem.

Numerical implementation
The implementation of the validation of an approximate zero of (69) can be found at the webpage Code page (2015).There, a complete instruction on how to run the codes can be found.(2015) and the red and green parts are computed using the conjugation maps P u,s harnessing the formula u(t) = P( u,s (t, θ)) with u,s (t, θ) the flow induced by h u,s for θ ∈ Bsym (1,1) (see Lemma 2.6).One time unit corresponds to L = 43.2034 We list the main parameters for the computations, which were chosen after numerical experimentation.As time rescaling parameter we choose L = 43.2034 in (68).For both the stable and unstable manifold we choose as domain radius in (6) ν = (1, 1).
For the validation of the connection we use the approach of van den Berg et al. (2015) in conjunction with the implementation of van den Berg and Sheombarsing (2015).The validated solution profiles are shown in Fig. 4.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 1
Fig. 1 Left co-existing hexagonal and trivial patterns in (2) for β = 1 and μ = − 2 135 .Right corresponding connecting orbit in (1) with u and v components in red and blue, respectively (Color figure online)

Lemma 2 . 1 1
Assume (14)  to be fulfilled.The map P is real-valued on the invariant subspace B Using the properties of the involution (11) in the second equality, together with θ * = θ (on B sym ν ) and a * = a (assumed) in the third, we compute

Fig. 2
Fig. 2 Top change of the validation radii while rescaling the eigenvectors (left m = (5, 5), right m =(15, 15)) We see that for a larger number of modes we obtain smaller error bounds.Note in addition that the larger the norm of ξ 1,2 is the bigger the uniform error bound r on B ν gets.Bottom dependence of the norm of the approximate solution on the number of rescalings, hence on the norms ξ 1,2 .These are an indicator for the size of the image in phase space

Fig. 3
Fig. 3 Top left change of validation radii while rescaling with factor μ = ( 6 7 , 7 6 ) for the choice m 1 = m 2 = 15.We observe qualitatively similar behavior to the uniform scaling.Top right change of validation radii while rescaling with factor μ = ( 6 7 , 7 6 ) for the choice m 1 = 5, m 2 = 15.We observe qualitatively different behavior to the uniform scaling.The maximal accuracy is obtained for ξ 2 ξ 1

Fig. 4
Fig.4Profiles of the u and v component of the rigorously verified connecting orbit.The blue part of either orbit is computed using Chebyshev series with the implementation from van den Berg and Sheombarsing (2015) and the red and green parts are computed using the conjugation maps P u,s harnessing the formula u(t) = P( u,s (t, θ)) with u,s (t, θ) the flow induced by h u,s for θ ∈ Bsym (1,1) (see Lemma 2.6).One time unit corresponds to L = 43.2034

Table 4
Expansion coefficients for Z k