Analytic continuation of local (un)stable manifolds with rigorous computer assisted error bounds

# Analytic continuation of local (un)stable manifolds with rigorous computer assisted error bounds

William D. Kalies Email: wkalies@fau.edu Shane Kepley S.K. partially supported by NSF grants DMS-1700154 and DMS - 1318172, and by the Alfred P. Sloan Foundation grant G-2016-7320 Email: skepley@my.fau.edu J.D. Mireles James J.M.J partially supported by NSF grants DMS-1700154 and DMS - 1318172, and by the Alfred P. Sloan Foundation grant G-2016-7320 Email: jmirelesjames@fau.edu
July 12, 2019
###### Abstract

We develop a validated numerical procedure for continuation of local stable/unstable manifold patches attached to equilibrium solutions of ordinary differential equations. The procedure has two steps. First we compute an accurate high order Taylor expansion of the local invariant manifold. This expansion is valid in some neighborhood of the equilibrium. An important component of our method is that we obtain mathematically rigorous lower bounds on the size of this neighborhood, as well as validated a-posteriori error bounds for the polynomial approximation. In the second step we use a rigorous numerical integrating scheme to propagate the boundary of the local stable/unstable manifold as long as possible, i.e. as long as the integrator yields validated error bounds below some desired tolerance. The procedure exploits adaptive remeshing strategies which track the growth/decay of the Taylor coefficients of the advected curve. In order to highlight the utility of the procedure we study the embedding of some two dimensional manifolds in the Lorenz system.

\SetLabelAlign

Center (#1)

## 1 Introduction

This paper describes a validated numerical method for computing accurate, high order approximations of stable/unstable manifolds of analytic vector fields. Our method generates a system of polynomial maps describing the manifold away from the equilibrium. The polynomials approximate charts for the manifold, and each comes equipped with mathematically rigorous bounds on all truncation and discretization errors. A base step computes a parameterized local stable/unstable manifold valid in a neighborhood of the equilibrium point. This analysis exploits the parameterization method [1, 2, 3, 4, 5, 6]. The iterative phase of the computation begins by meshing the boundary of the initial chart into a collection of submanifolds. The submanifolds are advected using a Taylor integration scheme, again equipped with mathematically rigorous validated error bounds.

Our integration scheme provides a Taylor expansion in both the time and space variables, but uses only the spatial variables in the invariant manifold. This work builds on the substantial existing literature on validated numerics for initial value problems, or rigorous integrators, see for example [7, 8, 9, 10], and exploits optimizations developed in [11, 12, 13].

After one step of integration we obtain a new system of charts which describe the advected boundary of the local stable/unstable manifold. The new boundary is adaptively remeshed to minimize integration errors in the next step. The development of a mathematically rigorous remeshing scheme to produce the new system of boundary arcs is one of the main technical achievements of the present work, amounting to a validated numerical verification procedure for analytic continuation problems in several complex variables. Our algorithm exploits the fact that the operation of recentering a Taylor series can be thought of as a bounded linear operator on a certain Banach space of infinite sequences (i.e. the Taylor coefficients), and this bounded linear operator can be studied by adapting existing validated numerical methods. The process of remeshing is iterated as long as the validated error bounds are held below some user specified tolerance, or a specified number of time units.

To formalize the discussion we introduce notation. We restrict the discussion to unstable manifolds and note that our procedure applies to stable manifolds equally well by reversing the direction of time. Suppose that is a real analytic vector field, and assume that generates a flow on an open subset . Let denote this flow.

Suppose that is a hyperbolic equilibrium point with unstable eigenvalues. By the unstable manifold theorem there exists an so that the set

 Wu\tiny loc(p0,f,r):={x∈Bnr(p0):Φ(x,t)∈Bnr(p0) for all t≤0},

is analytically diffeomorphic to a -dimensional disk which is tangent at to the unstable eigenspace of the matrix . Moreover, as for each . Here is the ball of radius about in . We simply write when and are understood. The unstable manifold is then defined as the collection of all points such that as which is given explicitly by

 Wu(p0)=⋃0≤tΦ(Wu\tiny loc(p0),t).

The first step of our program is to compute an analytic chart map for the local manifold of the form, , such that , is contained the unstable eigenspace, and

 image(P)⊂Wu\tiny loc(p0).

In Section 3 we describe how this is done rigorously with computer assisted a-posteriori error bounds.

Next, we note that is backward invariant under , and thus the unstable manifold is the forward image of the boundary of the local unstable manifold by the flow. To explain how we exploit this, suppose we have computed the chart of the local manifold described above.

We choose a piecewise analytic system of functions , , such that

 ⋃1≤j≤K0γj(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Bd−11(0))=∂P(Bd1(0)),

with

 image(γi)∩image(γj)⊂∂% image(γi)∩∂image(γj),

i.e. the functions , parameterize the boundary of the local unstable manifold, and their pairwise intersections are -dimensional submanifolds. Now, fix a time , and for each , , define by

 Γj(s,t)=Φ(γj(s),t)(s,t)∈¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Bd−11(0)×[0,T].

We note that

 image(P)∪⎛⎝⋃1≤j≤K0image(Γj)⎞⎠⊂Wu(p0),

or in other words, the flow applied to the boundary of the local unstable manifold yields a larger piece of the unstable manifold. Thus, the second step in our program amounts to rigorously computing the charts and is described in Section 4. Figure 1 provides a graphical illustration of the scheme.

Figure 2 illustrates the results of our method in a specific example. Here we advect the boundary of a high order parameterization of the local stable manifold at the origin of the Lorenz system at the classical parameter values, see Section 3. The color of each region of the manifold describes the integration time . The resulting manifold is described by an atlas consisting of 4,674 polynomial charts computed to order 24 in time and 39 in space. The adaptive remeshing described in Section 4.4 is performed to restrict to the manifold bounded by the rectangle .

###### Remark 1 (Parameterization of local stable/unstable manifolds).

Validated numerical algorithms for solving initial value problems are computationally intensive, and it is desirable to postpone as long as possible the moment when they are deployed. In the present applications we would like to begin with a system of boundary arcs which are initially as far from the equilibrium as possible, so that the efforts of our rigorous integrator are not spent recovering the approximately linear dynamics on the manifold. To this end, we employ a high order polynomial approximation scheme based on the parameterization method of [1, 2, 3]. For our purposes it is important to have also mathematically rigorous error bounds on this polynomial approximation, and here we exploit a-posteriori methods of computer assisted proof for the parameterization method developed in the recent work of [14, 4, 15, 16, 12]. These methods yield bounds on the errors and on the size of the domain of analyticity, accurate to nearly machine precision, even a substantial distance from the equilibrium. See also the lecture notes [17].

###### Remark 2 (Technical remarks on validated numerics for initial value problems).

A thorough review of the literature, much less any serious comparison of existing rigorous integrators, are tasks far beyond the scope of the present work. We refer the interested reader to the discussion in the recent review of [18]. That being said, a few brief remarks on some similarities and differences between the present and existing works are in order. The comments below reflect the fact that different studies have differing goals and require different tools: our remarks in no way constitute a criticism of any existing method. The reader should keep in mind that our goal is to advect nonlinear sets of initial conditions which are parameterized by analytic functions.

In one sense our validated integration scheme is closely related to that of [7], where rigorous Taylor integrators for nonlinear sets of initial conditions are developed. A technical difference is that the a-posteriori error analysis implemented in [7] is based on an application of the Schauder Fixed Point Theorem to a Banach space of continuous functions. The resulting error bounds are given in terms of continuous rather than analytic functions.

In this sense our integration scheme is also related to the work of [19, 10] on Taylor integrators in the analytic category. While the integrators in the works just cited are used to advect points or small boxes of initial conditions, the authors expand the flow in a parameter as well as in time, validating expansions of the flow in several complex variables. A technical difference between the method employed in this work and the work just cited is that our a-posteriori analysis is based on a Newton-like method, rather than the contraction mapping theorem.

The Newton-like analysis applies to polynomial approximations which are not required to have interval coefficients. Only the bound on the truncation error is given as an interval. The truncation error in this case is not a tail, as the unknown analytic function may perturb our polynomial coefficients to all orders. We only know that this error function has small norm.

This can be viewed as an analytic version of the “shrink wrapping” discussed in [20]. However, in our case the argument does not lose control of bounds on derivatives. Cauchy bounds can be used to estimate derivatives of the truncation error, after giving up a small portion of the validated domain of analyticity. Such techniques have been used before in the previous work of [15, 16]. The works just cited deal with Taylor methods for invariant manifolds rather than rigorous integrators.

Since our approach requires only floating point rather than interval enclosures of Taylor coefficients, we can compute coefficients using a numerical Newton scheme rather than solving term by term using recursion. Avoiding recursion can be advantageous when computing a large number of coefficients for a multivariable series. The quadratic convergence of Newton’s method facilitates rapid computation to high order. Note also that while our method does require the inversion of a large matrix, this matrix is upper triangular, and hence this inversion can be managed fairly efficiently.

Any discussion of rigorous integrators must mention the work of the CAPD group. The CAPD library is probably the most sophisticated and widely used software package for computer assisted proof in the dynamical systems community. The interested reader will want to consult the works of [9, 21]. The CAPD algorithms are based on the pioneering work of Lohner [22, 23, 24], and instead of using fixed point arguments in function space to manage truncation errors, develop validated numerical bounds based on the Taylor remainder theorem. The CAPD algorithms provide results in the category, and are often used in conjunction with topological arguments in a Poincare section [25, 26, 27, 28, 29, 30] to give computer assisted proofs in dynamical systems theory.

###### Remark 3 (Basis representations for analytic charts).

In this work we describe our method by computing charts for both the local parameterization and its advected image using Taylor series (i.e. analytic charts are expressed in a monomial basis). This choice allows for ease of exposition and implementation. However, the continuation method developed here works in principle for other choice of basis. What is needed is a method for rigorously computing error estimates.

Consider for example the case of an (un)stable manifold attached to a periodic orbit of a differential equation. In this case one could parameterize the local manifold using a Fourier-Taylor basis as in basis as in [31, 32, 33]. Such a local manifold could then be continued using Taylor basis for the rigorous integration as discussed in the present work. Alternatively, if one is concerned with obtaining the largest globalization of the manifold with minimal error bounds it could be appropriate to use a Chebyshev basis for the rigorous integration to reduce the required number of time steps. The point is that we are free to choose any appropriate basis for the charts in space/time provided it is amenable to rigorous validated error estimates. The reader interested in computer assisted proofs compatible with the presentation of the present work – and using bases other than Taylor – are referred to [11, 12, 34, 35, 13, 36]

###### Remark 4 (Why continue the local manifold?).

As just mentioned there are already many studies in the literature which give validated numerical computations of local invariant manifolds, as well as computer assisted proofs of the existence of connections between them. Our methods provide another approach to the computer assisted study of connecting orbits via the “short connection” mechanism developed in [6]. But if one wants to rule out other connections then it is necessary to continue the manifold, perhaps using the methods of the present work. Correct count for connecting orbits is essential for example in applications concerning optimal transport time, or for computing boundary operators in Morse/Floer homology theory.

###### Remark 5 (Choice of the example system).

The validated numerical theorems discussed in the present work are benchmarked for the Lorenz system. This choice has several advantages, which we explain briefly. First, the system is three dimensional with quadratic nonlinearity. Three dimensions facilitates drawing of nice pictures which provide useful insight into the utility of the method. The quadratic nonlinearity minimizes technical considerations, especially the derivation of certain analytic estimates. We remark however that the utility of the Taylor methods discussed here are by no means limited to polynomial systems. See for example the discussion of automatic differentiation in [37]. We note also that many of the computer assisted proofs discussed in the preceding remark are for non-polynomial nonlinearities. The second and third authors of the present work are preparing a manuscript describing computer assisted proofs of chaotic motions for a circular restricted four body problem which uses the methods of the present work.

Another advantage of the Lorenz system is that we exploit the discussion of rigorous numerics for stable/unstable manifolds given in the Lecture notes of [17]. Again this helps to minimize technical complications and allows us to focus instead on what is new here.

Finally, the Lorenz system is an example where other authors have conducted some rigorous computer assisted studies growing invariant manifolds attached to equilibrium solutions of differential equations. The reader wishing to make some rough comparisons between existing methods might consult the Ph.D. thesis [38], see especially Section . For example one could compare the results illustrated in Figure of that Thesis with the results illustrated in Figure of the present work. The manifolds in these figures have comparable final validated error bounds, while the manifold illustrated in Figure explores a larger region of phase space.

We caution the reader that such comparisons must be made only cautiously. For example the validation methods developed in [38] are based on topological covering relations and cone conditions, which apply in a setting. Hence the methods of [38] apply in a host of situations where the methods of the present work – which are based on the theory of analytic functions of several complex variables – breakdown. Moreover the initial local patch used for the computations in [38] is smaller than the validated local manifold developed in [17] from which we start our computations.

The remainder of the paper is organized as follows. In Section 2 we recall some basic facts from the theory of analytic functions of several complex variables, define the Banach spaces of infinite sequences used throughout the paper, and state an a-posteriori theorem used in later sections. In Section 3 we review the parameterization method for stable/unstable manifolds attached to equilibrium solutions of vector fields. In particular we illustrate the formalism which leads to high order polynomial approximations of the local invariant manifolds for the Lorenz system, and state an a-posteriori theorem which provides the mathematically rigorous error bounds. Section 4 describes in detail the subdivision strategy for remeshing analytic submanifolds and the rigorous integrator used to advect these submanifolds. Section 5 illustrates the method in the Lorenz system and illustrates some applications. The implementation used to obtain all results are found at [39].

## 2 Background: analytic functions, Banach algebras of infinite sequences, and an a-posteriori theorem

Section 2 reviews some basic properties of analytic functions, some standard results from nonlinear analysis, and establishes some notation used in the remainder of the present work. This material is standard and is included only for the sake of completeness. The reader may want to skip ahead to Section 3, and refer back to the present section only as needed.

### 2.1 Analytic functions of several variables, and multi-indexed sequence spaces

Let and . We endow with the norm

 ∥z∥=max1≤i≤d|z(i)|,

where is the usual complex modulus. We refer to the set

 Dd:={w=(w(1),…,w(d))∈Cd:|w(i)|<1 for all 1≤i≤d},

as the unit polydisk in . Throughout this paper whenever is understood we write . Note that the -dimensional open unit cube is obtained by restricting to the real part of .

Recall that a function is analytic (in the sense of several complex variables) if for each and , the complex partial derivative, , exists and is finite. Equivalently, is analytic (in the sense of several complex variables) if it is analytic (in the usual sense) in each variable with the other variables fixed, for . Denote by

 ∥f∥C0(D,C):=supw∈D|f(w(1),…,w(d))|,

the supremum norm on which we often abbreviate to , and let denote the set of bounded analytic functions on . Recall that if is a sequence of analytic functions and

 limn→∞∥f−fn∥∞=0,

then is analytic (i.e.  is a Banach space when endowed with the norm). In fact, is a Banach algebra, called the disk algebra, when endowed with pointwise multiplication of functions.

We write for a -dimensional multi-index, where is the order of the multi-index, and to denote raised to the -power. Recall that a function, if and only if for each , has a power series expansion

 f(w)=∑α∈Ndaα(w−z)α,

converging absolutely and uniformly in some open neighborhood with . For the remainder of this work, we are concerned only with Taylor expansions centered at the origin (i.e.  and ). Recall that the power series coefficients (or Taylor coefficients) are determined by certain Cauchy integrals. More precisely, for any and for any the -th Taylor coefficient of centered at is given explicitly by

 aα:=1(2πi)d∫|z(1)|=r…∫|z(d)|=rf(z(1),…,z(d))(z(1))α1+1…(z(d))αd+1dz(1)…dz(d),

where the circles , are parameterized with positive orientation.

The collection of all functions whose power series expansion centered at the origin converges absolutely and uniformly on all of is denoted by . Let denote the set of all -dimensional multi-indexed sequences of complex numbers. For define the norm

 ∥a∥1,d:=∑α∈Nd|aα|,

and let

 ℓ1d:={a∈Sd:∥a∥1,d<∞},

(i.e.  is the Banach space of all absolutely summable -dimensional multi-indexed sequences of complex numbers). When is understood we often abbreviate to and . For any with Taylor series centered at the origin given by

 f(z)=∑α∈Ndaαzα,

let denote the mapping given by

which associates an analytic function, , with the sequence of Taylor coefficients for its power series expansion at . We refer to as the Taylor transform of and note that is both linear, one-to-one, and takes values in . Moreover, we have the trivial bound

 ∥f∥∞≤∥T(f)∥1,

for each . Now, let denote the collection of all functions whose Taylor coefficients are in and note that we have the inclusions

 B1d⊂Bd⊂Cω(D).

In particular, if , then defines a unique analytic function, given by

 f(z)=∑α∈Ndaαzα.

We remark that if then extends uniquely to a continuous function on , as the power series coefficients are absolutely summable at the boundary. So if then is well defined, continuous on , and analytic on .

Finally, recall that inherits a Banach algebra structure from pointwise multiplication, a fact which is critical in our nonlinear analysis in Sections 3 and 4. Begin by defining a total order on by setting if for every and if (i.e. we endow with the lexicographic order). Given , define the binary operator by

 [a∗b]α=∑κ≺αaκ⋅bα−κ.

We refer to as the Cauchy product, and note the following properties:

• For all we have

 ∥a∗b∥1≤∥a∥1∥b∥1.

In particular, is a Banach algebra when endowed with the Cauchy product.

• Let , and suppose that

 f(z)=∑α∈Ndaαzα% andg(z)=∑α∈Ndbαzα.

Then and

 (f⋅g)(z)=∑α∈Nd[a∗b]αzα.

In other words, pointwise multiplication of analytic functions corresponds to the Cauchy product in sequence space.

###### Remark 6 (Real analytic functions in B1d).

If and the Taylor coefficients of are real, then is real analytic on and continuous on .

###### Remark 7 (Distinguishing space and time).

In Section 4 it is advantageous both numerically and conceptually to distinguish time from spatial variables. When we need this distinction we write with the appropriate norm given by

 ∥a∥1,d+1=∞∑m=0∑α∈Nd|am,α|.

In this setting, defines a unique analytic function given by

 f(z,t)=∞∑m=0∑α∈Ndam,αzαtm,

where is distinguished as the (complex) space variable and is the time variable. Analogously, we extend the ordering on multi-indices to this distinguished case by setting if and as well as the Cauchy product by

 [a∗b]m,α=∑j≤m∑κ≺αaj,κ⋅bm−j,α−κ.

### 2.2 Banach spaces and linear algebra

The validation methods utilized in this work are based on a set of principles for obtaining mathematically rigorous solutions to nonlinear operator equations with computer assistance referred to as the radii polynomial approach. A key feature of this philosophy is the characterization of a nonlinear problem in the space of analytic functions as a zero finding problem in sequence space. Specifically, our methods will seek a (Fréchet) differentiable map in and require (approximate) computation of this map and its derivative.

For our purposes, we are interested in bounded linear operators defined on . Let denote the vector space of bounded linear operators from to itself, which we shorten to , equipped with the operator norm induced by . For this discussion we utilize the notation with space/time distinguished. To avoid confusion over indices, we denote indices for linear operators inside square brackets and components of vectors outside square brackets. Now, we fix a basis for composed of where

 [ejκ]m,α=(1(j,κ)=(m,α)0otherwise),

and we specify an element , by its action on these basis vectors which we denote by

 Ajk=A⋅ejk

With this notation in place, our first goal is to compute a formula for the operator norm on defined by

 ||A||1=sup||h||=1||A⋅h||1.
###### Proposition 2.1.

For , the operator norm is given by

 ||A||1=sup(j,κ)∈N×Nd∣∣∣∣Ajκ∣∣∣∣1
###### Proof.

We define which is finite since is a bounded linear operator. Suppose is a unit vector which we express in the above basis as

 h=∞∑j=0∑κ∈Ndhj,kejκ.

Then for each we have

 |[A⋅h]m,α|=∣∣ ∣∣∞∑j=0∑κ∈Nd[Ajκ]m,α⋅hj,κ∣∣ ∣∣.

Applying this directly for each coordinate in leads to the following estimate

 ||A⋅h||1 =∞∑m=0∑α∈Nd∣∣ ∣∣∞∑j=0∑κ∈Nd[Ajκ]m,α⋅hj,κ∣∣ ∣∣ ≤∞∑m=0∑α∈Nd∞∑j=0∑κ∈Nd∣∣[Ajκ]m,α∣∣⋅∣∣hj,κ∣∣ ≤∞∑j=0∑κ∈Nd∣∣hj,κ∣∣∞∑m=0∑α∈Nd∣∣[Ajκ]m,α∣∣ ≤∞∑j=0∑κ∈Nd∣∣hj,κ∣∣∣∣∣∣Ajκ∣∣∣∣1 ≤C∞∑j=0∑κ∈Nd∣∣hj,κ∣∣ =C

and taking the supremum over all unit vectors in we have . Conversely, for any we may choose such that . It follows that

 ||A||1≥∣∣∣∣A⋅ej,κ∣∣∣∣1>C−ϵ

and we conclude that . ∎

Next, we define specific linear operators which play an important role in the developments to follow. The first operator is the multiplication operator induced by an element in . Specifically, for a fixed vector, , there exists a unique linear operator, , whose action is given by

 Ta⋅u=a∗u (1)

for every . With respect to the above basis we can write explicitly as

 [Tjκa]m,α=(aj−m,κ−α(m,α)≺(j,κ)0otherwise)

which can be verified by a direct computation. The second operator is a coefficient shift followed by padding with zeros, which we will denote by . Its action on is given explicitly by

 [η⋅u]m,α={0if m=0um−1,αif m≥1 (2)

Additionally, we introduce the “derivative” operator whose action on vectors will be denoted by . Its action on is given by the formula

 [u′]m,α={um,αif m=0mum,αif m≥1 (3)

The usefulness in these definitions is made clear in Section 4.

Finally, we introduce several properties of these operators which allow us to estimate their norms. The first is a generalization of the usual notion of a lower-triangular matrix to higher order tensors.

###### Proposition 2.2.

We say an operator, , is upper triangular with respect to if for every . Then, each of the operators defined above is upper triangular. The proof for each operator follows immediately from their definitions.

Next, we introduce notation for decomposing a vector into its finite and infinite parts. Specifically, for fixed we denote the finite truncation of to -many terms (embedded in ) by

 umα={uj,κ(j,κ)≺(m,α)0otherwise, (4)

and we define the infinite part of by . From the point of view of Taylor series, are the coefficients of a polynomial approximation obtained by truncating to temporal terms and spatial terms in the direction, and represents the tail of the Taylor series. With this notation we establish several useful estimates for computing norms in .

###### Proposition 2.3.

Fix and suppose is arbitrary. Then the following estimates hold for all .

 ||Ta⋅u||1≤ ||a||1||u||1 (5) ||η(u)||1= ||u||1 (6)

The proof is a straightforward computation.

### 2.3 Product spaces

In the preceding discussion we considered the vector space structure on and described linear operators on this structure. In this section, we recall that is an algebra, and therefore it is meaningful to consider vector spaces over where we consider elements of as “scalars”. Indeed, an -dimensional vector space of this form is the appropriate space to seek solutions to the invariance equation described in Section 3 as well as IVPs which we describe in Section 4. To make this more precise we define

 X=⎧⎨⎩{u(i)m,α}⊂Cd:∞∑m=0∑α∈Nd|u(i)m,α|<∞ for all 1≤i≤n⎫⎬⎭, (7)

and we recognize that an element defines a unique analytic function in -many variables, taking values in . We recall that the restriction of this function to a single coordinate defines a scalar analytic function with coefficients in . Thus, can be equivalently defined as an -dimensional vector space over given by

 X=ℓ1×ℓ1×…ℓ1n-% copies=(ℓ1)n, (8)

and a typical element takes the form with each . When solving nonlinear problems in , we will typically adopt the notation and point of view in Equation (8). Next, we equip with the norm given by

 ||u||X=max1≤i≤n{||u(i)||1}. (9)

Finally, define multiplication in componentwise. Specifically, if , then each is an -length vector of scalars from and the multiplication defined by

 [u∗v]m,α=([u(1)∗v(1)]m,α,…,[u(n)∗v(n)]m,α) (10)

makes into a Banach algebra. The decomposition of into a finite projection and infinite tail is also defined componentwise.

Let denote the vector space of linear operators on , and suppose . Since is a finite dimensional vector space over , it follows that for any fixed basis of over , we can identify with some square matrix, , so that the action of on a vector is left multiplication by which has the form

 Q=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝q(11)q(12)…q(1n)q(21)q(22)…q(2n)⋮⋮⋱⋮q(n1)q(n2)…q(nn)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

It is a standard result that the sup norm defined in Equation (9) induces the operator norm given by

 ||A||X=max1≤i≤n{ri: ri=n∑j=1∣∣∣∣q(ij)∣∣∣∣1}.

where are bounded linear operators and we recall that denotes their operator norm.

### 2.4 A-posteriori analysis for nonlinear operators between Banach spaces

The discussion in Section 2.1 motivates the approach to validated numerics/computer assisted proof adopted below. Let and consider a nonlinear operator (possibly with only densely defined). Suppose that we want to solve the equation

 Ψ(f)=0.

Projecting the components of into sequence space results in an equivalent map on the coefficient level. The transformed problem is truncated by simply restricting our attention to Taylor coefficients with order for some . We denote by the truncated map. The problem is now solved using any convenient numerical method, and we denote by the appropriate numerical solution, and by the infinite sequence which results from extending by zeros.

We would like now, if possible, to prove that there is an near , which satisfies . Should we succeed, then by the discussion in Section 2.1, the function with Taylor coefficients given by is a zero of as desired. The following proposition, which is formulated in general for maps between Banach spaces, provides a framework for implementing such arguments.

###### Proposition 2.4.

Let , be Banach spaces and be a Fréchet differentiable mapping. Fix and suppose there are bounded linear operators , , with one-to-one. Assume that there are non-negative constants, , satisfying the following bounds for all :

 (11) ||Id−AA†||X≤Z0 (12) (13) (14) Y0+(Z0+Z1)r+Z2r2

Then there exists a unique so that .

###### Proof.

Consider the nonlinear operator defined by

 T(x)=x−AF(x).

Since is one-to-one, is a zero of if and only if is a fixed point of . The idea of the proof is to use the Banach fixed point theorem to establish the existence of a unique fixed point in .

Let Id denote the identity map on , suppose , and note that . Then

Taking this together with assumptions (12), (13), and (14) we obtain the bound

 supx∈¯¯¯¯¯¯¯¯¯¯¯Br(¯a)||DT(x)||X≤Z0+Z1+Z2r. (16)

Now, if , then applying the bound (11) and invoking the Mean Value Theorem yields the estimate

 ≤||T(x)−T(¯a)||X+||T(¯a)−¯a||X ≤Y0+(Z0+Z1)r+Z2r2

where the last inequality is due to Equation (15). This proves that maps into itself. In fact, sends into , by the strict inequality.

Finally, assume and apply the bound of Equation (16) with the Mean Value Theorem once more to obtain the contraction estimate

 ||T(x)−T(y)||X ≤supx∈¯¯¯¯¯¯¯¯¯¯¯¯Br(¯a)||DT(x)||X⋅||x−y||X ≤(Z0+Z1+Z2r)||x−y||X <(1−Y0r)||x−y||X (18)

where the second to last line follows from another application of Equation (15) and the last line from noticing that . Since , the Contraction Mapping Theorem is satisfied on . By the strict inequality of Equation (17) we conclude that has a unique fixed point and it follows that is the unique zero of in . ∎

###### Remark 8.

A few remarks on the intuition behind the terms appearing in the proposition are in order. Intuitively speaking, occurs when , , are small, and is not too large. Here measures the defect associated with (i.e.  small means that we have a “close” approximate solution). We think of as an approximation of the differential , and as an approximate inverse of . Then measure the quality of these approximations. These approximations are used as it is typically not possible to invert exactly. Finally is in some sense a measure of the local “stiffness” of the problem. For example is often taken as any uniform bound on the second derivative of near . The choice of the operators is problem dependent and best illustrated through examples. Finally we remark that it is often unnecessary to specify explicitly the space . Rather, what is important is that for all we have that and that .

###### Remark 9.

Following [40, 41, 42, 43], we exploit the radii polynomial method to organize the computer assisted argument giving validated error bounds for our integrator. In short, this amounts to rewriting the contraction mapping condition above by defining the radii polynomial

 p(r)=Z2r2+(Z0+Z1−1)r+Y0

and noting that the hypothesis of Proposition 2.4 in equation (15) is satisfied for any such that . It follows that the minimum root of (if it exists) gives a sharp bound on the error, and if has distinct roots, , then on the entire interval . The isolation bound is theoretically infinite, as the solutions of initial value problems are globally unique. However the width of the interval provides a quantitative measure of the difficulty of a given proof, as when this difference is zero the proof fails.

## 3 The parameterization method for (un)stable manifolds

The parameterization method is a general functional analytic framework for analyzing invariant manifolds, based on the idea of studying dynamical conjugacy relationships. The method is first developed in a series of papers [1, 2, 3, 44, 45, 46]. By now there is a small but thriving community of researchers applying and extending these ideas, and a serious review of the literature would take us far afield. Instead we refer the interested reader to the recent book [37], and turn to the task of reviewing as much of the method as we use in the present work.

Consider a real analytic vector field , with generating a flow , for some open set . Suppose that is an equilibrium solution, and let denoted the stable eigenvalues of the matrix . Let denote a choice of associated eigenvectors. In this section we write , for the unit ball in .

The goal of the Parameterization Method is to solve the invariance equation

 f(P(s))=λ1s1∂∂s1P(s)+…+λdsd∂∂sdP(s), (19)

on , subject to the first order constraints

 P(0)=pand∂∂sjP(0)=ξj (20)

for . From a geometric point of view, Equation (19) says that the push forward by of the linear vector field generated by the stable eigenvalues is equal to the vector field restricted to the image of . In other words Equation (19) provides an infinitesimal conjugacy between the stable linear dynamics and the nonlinear flow, but only on the manifold parameterized by . More precisely we have the following Lemma.

###### Lemma 3.1 (Parameterization Lemma).

Let be the linear flow

 L(s,t)=(eλ1ts1,…,eλdtsd).

Let be a smooth function satisfying Equation (19) on and subject to the constraints given by Equation (20). Then satisfies the flow conjugacy

 Φ(P(s),t)=P(L(s,t)), (21)

for all and .

For a proof of the Lemma and more complete discussion we refer to [47]. The flow conjugacy described by Equation (21) is illustrated pictorially in Figure 3. Note that is the flow generated by the vector field

 ddtsj=λjsj,1≤j≤d,

i.e. the diagonal linear system with rates given by the stable eigenvalues of . Note also that the converse of the lemma holds, so that satisfies the flow conjugacy if and only if satisfies the infinitesimal conjugacy. We remark also that is (real) analytic if is analytic [2, 3].

Now, one checks that if satisfies the flow conjugacy given Equation (21), then

 P(B)⊂Ws(p),

i.e. the image of is a local stable manifold. This is seen by considering that

 limt→∞Φ(P(s),t)=limt→∞P(L(s,t))=p%foralls∈B⊂Rd,

which exploits the flow conjugacy, the fact that is stable linear flow, and that is continuous.

It can be shown that solutions of Equation (19) are unique up to the choice of the scalings of the eigenvectors. Moreover, on the level of the power series representation of the solution, the scaling of the eigenvectors determines the decay rates of the Taylor coefficients of . Proofs are found for example in [3]. These facts are used to show that, once we fix the domain of the praameterization to , the solution parameterizes a larger or smaller local portion of the stable manifold depending only on the choice of the eigenvector scalings. In practice this freedom in the choice in the scalings of the eigenvectors is exploited to stabilize numerical computations. See for example [16].

The existence question for Equation (19) is somewhat more subtle. While the stable manifold theorem guarantees the existence of stable manifolds for a hyperbolic fixed point, Equation (19) provides more. Namely a chart map which recovers the dynamics on the invariant manifold via a flow conjugacy relation. It is not surprising then that some additional assumptions are necessary in order to guarantee solutions of Equation (19).

The necessary and sufficient conditions are given by considering certain non-resonance conditions between the stable eigenvalues. We say that the stable eigenvalues are resonant if there exists an so that

 α1λ1+…+αdλd=λjfor some 1≤j≤d. (22)

The eigenvalues are non-resonant if the condition given in Equation (22) fails for all . Note that since , , all have the same sign, there are only a finite number of opportunities for a resonance. Thus, in spite of first appearances, Equation (22) imposes only a finite number of conditions between the stable eigenvalues. The following provides necessary and sufficient conditions that some solution of Equation (19) exists.

###### Lemma 3.2 (A-priori existence).

Suppose that are non-resonant. Then there is an such that

 ∥ξj∥≤ϵfor each 1≤j≤d,

implies existence of a solution to Equation (19) satisfying the constraints given by Equation (20).

A proof of a substantially more general theorem for densely defined vector fields on Banach spaces (which certainly covers the present case) is found in [48]. Other general theorems (for maps on Banach spaces) are found in [1, 2, 3]. We note that in applications we would like to pick the scalings of the eigenvectors as large as possible, in order to parameterize as large a portion of the manifold as possible, and in this case we have no guarantee of existence. This is motivates the a-posteriori theory developed in [48, 49, 16], which we utilize in the remainder of the paper.

Finally, we note that even when the eigenvalues are resonant it is still possible to obtain an analogous theory by modifying the map . As remarked above, there can only be finitely many resonances between . Then in the resonant case can be chosen a polynomial which “kills” the resonant terms, i.e. we conjugate to a polynomial rather than a linear vector field in . Resonant cases are treated in detail in [1, 15]. Of course all the discussion above goes through for unstable manifolds by time reversal, i.e. considering the vector field .

### 3.1 Formal series solution of equation (19)

In practical applications our first goal is to solve Equation (19) numerically. Again, it is shown in [3] that if is analytic, then is analytic as well. Based on the discussion of the previous section we look for a choice of scalings of the eigenvectors and power series coefficients so that

 P(s)=∑α∈Ndpαsα, (23)

is the desired solution for .

Imposing the linear constraints given in Equation (20) leads to

 p0=pandpαj=ξj for 1≤j≤d.

Here denotes the zero multi-index in , and for are the first-order multi-indices satisfying . The remaining coefficients are determined by power matching. Note that

 λ1s1∂∂s1P(s)+…+λdsd∂∂sdP(s)=∑α∈Nd(α1λ1+…+αdλd)pαsα.

Returning to Equation (19) we let

 f[P(s)]=∑α∈Ndqαsα,

so that matching like powers leads to the homological equations

 (α1λ1+…+αdλd)pα−qα=0,

for all . Of course each depends on in a nonlinear way, and solution of the homological equations is best illustrated through examples.

Example: equilibrium solution of Lorenz with two stable directions. Consider the Lorenz system defined by the vector field where

 f(x,y,z)=⎛⎜⎝σ(y−x)x(ρ−z)−yxy−βz⎞⎟⎠. (24)

For there are three equilibrium points

 p0=⎛⎜⎝000⎞⎟⎠,andp±=⎛⎜ ⎜⎝[]c±√β(ρ−1)±√β(ρ−1)ρ−1⎞⎟ ⎟⎠.

Choose one of the three fixed points above and denote it by . Assume that has two eigenvalues