PARAMETERIZATION OF SLOW-STABLE MANIFOLDS AND THEIR INVARIANT VECTOR BUNDLES: THEORY AND NUMERICAL IMPLEMENTATION

The present work deals with numerical methods for computing slow stable invariant manifolds as well as their invariant stable and unstable normal bundles. The slow manifolds studied here are sub-manifolds of the stable manifold of a hyperbolic equilibrium point. Our approach is based on studying certain partial differential equations equations whose solutions parameterize the invariant manifolds/bundles. Formal solutions of the partial differential equations are obtained via power series arguments, and truncating the formal series provides an explicit polynomial representation for the desired chart maps. The coefficients of the formal series are given by recursion relations which are amenable to computer calculations. The parameterizations conjugate the dynamics on the invariant manifolds and bundles to a prescribed linear dynamical systems. Hence in addition to providing accurate representation of the invariant manifolds and bundles our methods describe the dynamics on these objects as well. Example computations are given for vector fields which arise as Galerkin projections of a partial differential equation. As an application we illustrate the use of the parameterized slow manifolds and their linear bundles in the computation of heteroclinic orbits.


Introduction.
We examine the question of numerically computing the stable and unstable normal bundles of some invariant manifolds of a differential equation. Critical to this work is the study of certain invariance equations which lead to numerical schemes for computing the invariant manifolds and bundles. Motivation for this work comes from considering stable/unstable manifolds associated with equilibria in finite dimensional Galerkin approximations of a parabolic partial differential equation. This is a situation where the dynamics of a high dimensional stable/unstable manifold is captured by a smaller number of "slow" variables. The main significance of the slow stable manifold is this: generic orbits in the stable manifold approach the equilibrium along the slow stable manifold to leading order. In this context, the invariant stable normal bundle of the slow invariant manifold provides control over the full manifold in directions normal to the slow manifold. Hence, the bundles we compute give a relatively simple, yet accurate approximation of the dynamically most important part of the full stable manifold.
Invariant manifolds and their normal bundles are fundamental objects of study in dynamical systems theory, going at least back to the seminal work of [12,13]. Given a particular dynamical systems one may be interested in numerical computation of a particular invariant stable/unstable normal bundle associated with an invariant manifold, and indeed a rich literature is devoted to this topic. A comprehensive survey is beyond the scope of the present work and we refer the interested reader to the works of [15,16,14,17,18,7] and the references therein.
Our approach is based on the parameterization method of [4,5,6]. In particular we follow [6] in defining the "slow manifold" associated with a hyperbolic equilibrium. The possibility of numerically computing these slow manifolds using the parameterization method was already mentioned in [6]. In addition to implementing the scheme described in [6] we extend the parameterization method in order to compute invariant stable and unstable vector bundles associated with the slow manifold. We observe p 0 The slow stable manifold of p 0 together with a depiction of the normal bundle at a point p 1 along the slow stable manifold of the equilibrium p 0 . The normal bundle decomposes into invariant stable and unstable normal bundles. The decomposition of the normal space is depicted at the point p 1 and the stable/unstable normal spaces are denoted by N s p 1 , N u p 1 and colored red and blue respectively. The orange surface depicts the full stable manifold of p 0 , of which the green curve (the slow stable manifold) is a submanifold. The essential feature of the stable normal bundle is that the full stable manifold is tangent to the stable normal bundle along the slow manifold. that chart maps for the invariant vector bundles solve a linearized version of the invariance equations discussed int [4,5,6]. We develop formal series expansions for solution of this linear equation and implement numerical methods for computing the manifolds and bundles in specific example problems. We also discuss an application of this framework to the numerical computation of connecting orbits in high dimensional Galerkin approximations of a partial differential equations.
The geometric objects studied in the present work are illustrated in Figure 1.1. In this figure p 0 represents an equilibrium of a differential equation x = f (x) with f : R n → R n . The green curve represents the local slow stable manifold W slow (p 0 ), which is an invariant manifold tangent to linear space spanned by the eigenvectors corresponding to the the slowest stable eigenvalues (see Section 2.2 for more details). For each point p 1 on W slow (p 0 ) we attach a frame for the ambient space R n . The union of all of these frames will be a frame bundle (i.e. the frames are connected to one another "smoothly"). The frame attached to each point decomposes into tangent (to the slow stable manifold W slow (p 0 )) and normal subspaces.
Moreover, the normal subspace is further decomposed into "stable" and "unstable" normal subspaces N s p1 and N u p1 , in the sense that the linearized flow about W slow (p 0 ) takes tangent vectors to tangent vectors, stable normal vectors to stable normal vectors, and unstable normal vectors to unstable normal vectors. The end result is that attached to each point of W slow (p 0 ) we obtain a system of coordinates for R n that is compatible with the dynamics.
Then for all p 1 ∈ W slow (p 0 ) the normal subspace N s p1 is tangent to the (full) stable manifold W s (p 0 ), as illustrated in Figure 1.1. The stable normal bundle {N s p1 : p 1 ∈ W slow (p 0 )} thus provides an approximation of the full stable manifold of p 0 in a neighborhood of the slow manifold. Since the accuracy of the stable normal bundle remains high at considerable distance from the equilibrium, it can act as a computationally convenient boundary condition in computing connecting orbits, which generically enter along the slow direction(s), see Section 4.4 and Remark 4.2.
Additionally, the unstable invariant bundle provides information on the most expanding directions about W slow (p 0 ), and, perhaps more importantly, provides "good coordinates" which can be used in more theoretical studies of the full stable manifold of p 0 . For example, we expect these invariant bundles to provide appropriate coordinates for a graph transform description of infinite dimensional manifolds in a parabolic PDE setting. Furthermore, an accurate approximation of the stable bundle may be of interest when computing the Evans function, when analyzing the stability of a stationary front of a PDE (described by a connecting orbit of the associated ODE). Finally, these coordinates appear suitable for describing generic orbits that pass close to the equilibrium and could for example be useful when studying periodic orbits in "fast-slow" systems. In the present paper we focus on using the stable bundles for approximating the full stable manifold (see Section 4), and we leave the other, somewhat more speculative aspects, for future research.
Modus Operandi: The main theme of the sequel is that the invariant manifolds and invariant bundles discussed above are characterized by certain invariance equations. Solving these invariance equations via high order Taylor/power matching schemes leads to efficient numerical methods for computing the invariant objects.
We point out two strong connections to related work. First, our approach is the analogue for slow stable manifolds of the Floquet normal form associated with a periodic orbit. Many authors have used the Floquet normal form as a computational tool in dynamical systems theory, and we refer the interested reader to [8,9] for more through discussion of the literature. Second, there is a strong link to the work in [1], where stable manifolds with resonant eigenvalues are studied using conjugation to normal forms via power series. The main difference is that we enrich the purely algebraic description in terms of power series with a dynamic interpretation, namely in terms of dynamically invariant (expanding and contracting) bundles.
The organization of this paper is as follows. In Section 2 we review the background material necessary for the remainder of the paper and study the conjugacy equations for the slow stable manifold and its invariant bundles. In Section 3 we develop a power series formalism for the equations of Section 2. Section 4 is devoted to numerical implementation and applications. Indeed, the utility of the methods of the present work is best illustrated through example computations. The methods facilitate the study of hyperbolic equilibria having a small small number of "slow" stable eigenvalues and a large number of "fast" stable eigenvalues (i.e. eigenvalues with real part close to and far from zero respectively). This situation occurs naturally in two everyday applications: one is in the setting of the fast-slow dynamical systems of geometric perturbation theory. A second application is the numerical study of equilibria solutions of parabolic partial differential equations. In the present work we focus on this second setting, and only remark that applying the methods of the present work in the context of geometric singular perturbation theory provides an interesting subject for further study.
For parabolic partial differential equations the typical situation is that a hyperbolic equilibrium will have a finite-dimensional unstable manifold, say of dimension k, and an infinite dimensional stable manifold of co-dimension k (due to the temporal smoothing of a parabolic flow). When some spectral/Galerkin projection of the PDE is truncated to a high enough order the resulting system of ordinary differential equations has a stable manifold of finite but large dimension. However (again due to the smoothing) we expect that the stable eigenvalues "decay rapidly," so that a slow manifold governs the local stable dynamics. In Section 4 we provide a number of example computations in the context of Fisher's partial differential equation model of the spatio-temporal spread of ecological information. We compute some slow stable manifolds and invariant frames for finite dimensional projections of the Fisher's Equation.
We also illustrate how the slow manifolds and frame bundles can be used in order to compute heteroclinic connecting orbits. The computer programs used in this paper are available for download at [24].

Conjugacy Equations for Stable Invariant Manifolds and Their
Invariant Bundles.
2.1. Background and Notation. We endow R n with the supremum norm We consider a real analytic vector field f : R n → R n and assume that f generates a globally defined flow Φ : R n × R → R n . Let x 0 ∈ R n , t > 0, and let γ : [0, t] → R n denote the solution of the initial value problem i.e. γ(t) = Φ(x 0 , t). The Fréchet derivative of Φ with respect to x 0 , which we denote by D x Φ(x 0 , t) ≡ M (t) satisfies the non-autonomous linear matrix valued initial value problem Note that for each fixed t the map Φ(x 0 , t) is bijective (the inverse is obtained by backward time flow). Since Φ is a differentiable mapping it follows that Φ is a diffeomorphism, hence M (t) is an isomorphism. We also utilize the fact that if v ∈ R n , then the derivative of Φ(x 0 , t) in the direction v, which we denote by It is worth remarking that this directional derivative m(t) itself satisfies the initial value problem Let p ∈ R n be a hyperbolic equilibrium point for f . In order to simplify the exposition we assume that Df (p) is real diagonalizable. The analysis can be extended to the case where Df (p) is diagonalizable with complex eigenvalues, but the notation becomes more involved. Hence, in order to elucidate the main issues, we opt to restrict to the case of real eigenvalues. For nondiagonalizable linearizations the analysis becomes more involved, see Remark 3.1. Let spec(Df (p 0 )) = {λ 1 , . . . , λ n } denote the eigenvalues of Df (p 0 ). We assume that the stable part of spec(Df (p 0 )) can be partitioned as In particular, the equilibrium is assumed to be hyperbolic. We refer to the subcollections as the "slow stable", "fast stable", "stable", and "unstable" eigenvalues, respectively. The cut-off between L ss and L f s depends on applications and is here made arbitrarily. For each 1 ≤ i ≤ n let ξ i denote a choice of eigenvector associated with the eigenvalue λ i . We take Λ to be the k × k diagonal matrix with the slow stable eigenvalues as diagonal entries and zeros elsewhere, and Ω to be the n × n diagonal matrix having all the eigenvalues as diagonal entries and zeros elsewhere. Let and denote respectively the n × k matrix whose columns are the "slow stable eigenvectors" and the n × n matrix whose columns are all n of the eigenvectors. The columns of Q 0 span R n and the columns of A 0 are linearly independent.
Since Df (p) = Q 0 ΩQ −1 0 it follows that . We refer to these as the slow stable, fast stable, stable, and unstable eigenspaces respectively.

Parameterization of Stable and Sub-Stable
Manifolds. Throughout this section we enforce the notation established in Section 2.1. The goal of the parameterization method is to study chart maps for invariant manifolds which also satisfy certain conjugacy equations. In the present work we focus on the invariant submanifold of the (local) stable manifold of an equilibrium that is tangent to the slow stable eigenspace E ss at the equilibrium. We call this a sub-stable manifold.
Our goal is to find a smooth function a : B k ν (0) → R n , having a(0) = p, Da(0) = A 0 , and satisfying the conjugacy relationship for all φ ∈ B k ν (0). The geometric meaning of this conjugacy is illustrated in Figure  2  The dynamics in the parameter space generated by exponentiating the matrix of stable eigenvalues Λ. The dynamics in phase space are generated by the flow Φ associated with the vector field f . The diagram "commutes" in the sense that applying first the chart map a and then nonlinear flow Φ is required to be the same as applying the linear dynamics e Λ and then the chart map a. The result is that the dynamics on the local stable manifold are described by the stable linear dynamics.
Note that any function a(φ) satisfying Equation (2.4) is one to one. To see this observe that a is tangent to the slow stable eigenspace at the origin (i.e. Da(0) = [ξ 1 | . . . |ξ k ] = A 0 ) and recall that A 0 is of full rank as its columns are linearly independent. By the implicit function theorem a is of rank k, and hence one-to-one, in some neighborhood B k r (0) ⊂ B k ν (0). Now suppose that φ 1 , φ 2 ∈ B k ν (0) and that a(φ 1 ) = a(φ 2 ). Then for any t ∈ R, Φ[a(φ 1 ), t] = Φ[a(φ 2 ), t] by the uniqueness of the initial value problem. Choose T > 0 so large that e ΛT φ 1 , e ΛT φ 2 ∈ B k r (0). By the conjugacy relation we have that a e Λt φ 1 = a e Λt φ 2 , and since the arguments are in B k r (0), the local immersion gives that e ΛT φ 1 = e ΛT φ 2 . But e ΛT is an isomorphism and we have The utility of Equation (2.4) is limited by the appearance of the flow Φ in the equation. In practice the flow is only known implicitly, i.e., it is determined by solving the differential equation. The parameterization method of [4,5,6] is based on the observation that there is a convenient infinitesimal version of Equation (2.4). This observation is encapsulated in the following Lemma. The proof is elementary, yet helps to motivate the definitions and discussion in the next section on invariant vector bundles. and Then a(φ) satisfies the conjugacy given in Equation (2.4) if and only if a is a solution of the partial differential equation By using (2.6) one can derive that γ is the solution of the initial value problem and differentiate both sides with respect to t in order to obtain Taking the limit as t → 0 gives that a(φ) is a solution of Equation (2.6). This completes the proof.
(i) It follows that if a : B k ν (0) → R n is a smooth function satisfying the linear constraints given by (2.5) as well as the partial differential equation (2.6), then the image of a is an immersed forward invariant disk in W s (p). Moreover a conjugates the dynamics on the range of a to the linear dynamics generated by Λ, in the sense of Equation (2.4). In particular, under the flow Φ the range of a accumulates on p. Hence computing a solution of the PDE subject to the linear constraints gives a means of computing invariant sub-manifolds of W s (p). (ii) Conditions providing for the existence of an analytic map a(φ) satisfying Equation (2.6) under the constraints given by Equation (2.5) are discussed in [4]. In fact it is necessary and sufficient that no "resonance" of the form occurs, in order that there exists an a(φ) conjugating the nonlinear flow to the linear flow generated by Λ. Here λ 1 , . . . , λ k ∈ E ss , λ j ∈ E s , and α = (α 1 , . . . , α k ) ∈ N k is any positive multi-index. Since the non-resonance condition relates only stable eigenvalues and positive multi-indices it reduces to a finite number of conditions. See [4,3] for a more complete discussion.
In the present work we assume a-priori that we have an analytic solution of Equation (2.6) in hand. (iii) Note that the choice of scalings for the eigenvectors in A 0 is free. There are two important remarks in this regard. First one can show that the solution equation (2.6) is unique up to the choice of these scalings [4]. The second point has to do with the power series coefficients of the parameterization map a, a topic discussed in more detail in Section 4.1. The point in this regard is that the choice of the of eigenvector scalings determines the decay rate of the power series coefficient [4]. This fact can be exploited in order to stabilize numerical computations. In particular one often chooses the scalings so that the last coefficient computed has magnitude below some prescribed tolerance.

Parameterization of Invariant Linear
Bundles and Frames. The previous section described the parameterization of local sub-stable manifolds. We now develop a similar theory for parameterization of linear invariant bundles and frames associated with sub-stable manifolds. Our discussion of these invariant bundles is aided by the assumption that we have in hand the parameterization of the manifold discussed above. The fact that the dynamics on the manifold are conjugate to some linear dynamics in the parameter space simplifies substantially the formulation of the conjugacy equations for the invariant bundles.
Throughout this section a : B k ν (0) → R n is a smooth solution of Equation (2.6) satisfying the linear constraints given by Equation (2.5). The next definition specializes the general notion of an invariant vector bundle to suit the needs of the present work. Recall that M (t) is the solution of the variational equation and

Definition 2.3 (Parameterization of a Rank One Invariant Vector Bundle).
Let v : B k ν (0) → R n be a smooth function. We say v parameterizes a rank one, forward invariant vector bundle over a B k ν (0) with exponential rate µ if v satisfies the equation for all t > 0 and φ ∈ B k ν (0). If µ < 0 we say that the vector bundle is exponentially contracted, and exponentially expanded if µ > 0.
The conjugacy described by Equation (2.9) is illustrated graphically in Figure 2.2. For fixed φ ∈ B k ν (0) we think of the vector v(φ) as being attached to the point a(φ) on the manifold image(a). For any t ≥ 0 we evoke Equation (2.4) and have that the flow takes a(φ) to a e Λt φ . Equation (2.9) now asks that the variational flow take the vector v(φ) at a(φ) to the vector e µt v e Λt φ at a e Λt φ . We remark that the notion of invariant vector bundle described by Equation (2.9) is a very special case of the general definition of vector bundle invariance, and one which relies heavily on the assumption that we understand the dynamics on the manifold up to a conjugacy to a linear flow. However we have no need of the general notion in the present work.
We also remark that Equation (2.9) describes bundle invariance in terms of the variational flow M (t), which itself is defined only explicitly by the initial value problem Equation (2.1). In analogy with the discussion of the parameterization of the invariant manifold, we formulate an infinitesimal version of the invariant bundle equation: We refer to Equation (2.10) as the invariant bundle equation.
Taking the limit as t → 0 we infer that v solves Equation (2.10).

Remark 2.5 (Rank-Invariant Bundles and Invariant Frames).
If v 1 , . . . , v are linearly independent solutions of the invariant bundle equation associated with rates µ 1 , . . . , µ then we say that parameterizes a rank-forward invariant bundle. If = n we say that V parameterizes an invariant frame bundle.
One can check by differentiating Equation (2.6) that the tangent vectors The equations for q i , i = 1, . . . , k are decoupled due to the linearity and diagonalization of the flow in parameter space. In addition the tangent vectors are seen to have "asymptotic direction" given by by considering the linear constraints in Equation (2.5). Hence, in addition to solving the invariance equation (2.10) with rates λ 1 , . . . , λ k , the tangent vectors q i (φ) accumulate on the slow stable eigenspace under the flow. Motivated by this remark we ask if it is possible to find more parameterized rank one, forward invariant bundles, satisfying similar asymptotic direction conditions? More precisely we are interested finding a collection of n rank one parameterized vector bundles q i ; each having exponential rate given by the eigenvalue λ i and having asymptotic direction given by the associated eigenvector ξ i . The following Theorem enumerates some basic properties of the desired collection. Then: • For all φ ∈ B k ν (0) the collection of vectors q 1 (φ), . . . , q n (φ) span R n . • For all t ≥ 0 and φ ∈ B k ν (0) the derivative of the flow along the slow stable manifold factors as M (t) = Q e Λt φ e Ωt Q −1 (φ). (2.12) Proof. Note that it follows from the assertions (in particular (2.12)) that q i (φ), φ ∈ B k ν (0) is a rank one forward invariant vector bundle with exponential rate λ i , see Definition 2.3. Hence, from (2.9) we see that by applying M (t) to Q(φ) we obtain the useful conjugacy = e λ1t q 1 e Λt φ | . . . |e λnt q n e Λt φ = Q e Λt φ e Ωt . (2.13) , and we recall that the eigenvectors ξ i are linearly independent. Since the q i are smooth functions, it follows that there is an 0 < r ≤ ν so that linear independence hold for all φ in the neighborhood B k r (0) ⊂ B k ν (0) (i.e. the determinant of a matrix is a smooth function of its entries).
Chooseφ ∈ B k ν (0), c ∈ R n , and suppose that Then for any t ≥ We now choose T > 0 so large that e ΛTφ ∈ B k r (0). Then Q(e ΛT φ) is an isomorphism by the arguments above, and we conclude that e Ωt c = 0, hence c = 0 since e ΩT is also an isomorphism. This shows that q 1 (φ), . . . , q n (φ) are linearly independent in B k ν (0). Finally, we note that Equation (2.12) follows from Equation (2.13) once we know that Q(φ) is invertible.
Remark 2.7. There is a close analogy to the structure of the equations that appear in the description of normal bundles of periodic orbits via Floquet theory (with the slow stable manifold replacing the periodic orbit). In view of this analogy, we refer to Equation (2.14) as the slow manifold Floquet normal form Equation. We refer to [8,9] for details on the computational implementation of Floquet theory to find invariant manifolds and bundles associated to periodic orbits. We note that although (2.14) is a matrix equation, it can be solved column by column if the nonresonance conditions discussed in Section 3 are fulfilled.

Remark 2.8 (Stable Manifold Tangent to the Stable Normal Bundle Along the Slow Manifold
). An application of the Fenichel theory can be used in order to show that the slow-stable manifold W slow loc (p 0 ) parameterized by a has its own stable manifold. Since points on image(a) accumulate at p 0 it follows that the stable manifold of W slow loc (p 0 ) is a subset of W s (p 0 ). In fact a dimension count shows that the union of W slow loc (p 0 ) and its local stable manifold are a local stable manifold for p 0 . Moreover the stable normal bundle parameterized by q 1 , . . . , q k gives the linear approximation of the stable manifold of W slow loc (p 0 ). In particular, the stable normal bundle is tangent to the stable manifold of p 0 along W slow loc (p 0 ), and is a quadratically good approximation of W s (p 0 ) in a neighborhood of the slow stable manifold. We do not delve into a proof of these assertions (which involve arguments similar to the ones laid out in detail above). Rather, we provide a numerical illustration in Section 4.3.

Computation of the Linear Bundles of the Slow Stable Manifold: Power Series Solution of Equation 2.14.
We begin the discussion by assuming that we have a power series solution of Equation (2.6). More precisely let a : R k → R n be the parameterization of the slow manifold given by Assume that the coefficients of a solve Equation (2.6) "order by order" in the sense of power series. This procedure is discussed in detail in Section 4.1 for the discretized Fisher's Equation.
For more discussion of the computation of the parameterization coefficients the interested reader is also referred to the derivations in [3,20,21,23,22,10,11]. For now we simply assume that the coefficients a α are known explicitly. Let Df : R n → GL(R n ) denote the Jacobian matrix of the vector field as a function of position and. Let Df [a(φ)] have power series representation with A 0 = Df (0), and A α ∈ Mat n×n for all α ∈ N k . Again we remark that in applications the coefficients in the expansion of Df [a(φ)] are worked out via automatic differentiation once the expansion of a and the specific form of the vector field f are given. An example is discussed in detail in Section 4.2 for the truncated Fishers Equation.
Again, for the moment we assume that the coefficients A α are known. Proceeding formally, we expand the unknown solution of Equation (2.14) as a power series 12 with Q 0 = Q(0). Then the Floquet Equation (2.14) becomes Matching like powers of φ gives from which we isolate the Q α term in order to obtain the linear homological equation for the matrix coefficient Q α . Here we have used that fact that A 0 = Df (p) and defined the matrices Recall that the differential is diagonalized by Making the change of variables Introducing the notation we solve Equation (3.4) column by column and obtain the matrix equations Recalling that Ω is the diagonal matrix of eigenvalues, we now assume that the nonresonance conditions hold. Then, solving Equation (3.5) row by row, we obtain These recursion relations let us compute Q to any desired order. Remark 3.1.
1. Since λ 1 , . . . , λ k all have the same (negative) sign, the non-resonance conditions (3.6) reduce to finitely many conditions. Explicitly, resonances can occur only for |α| ≤ (λ n + |λ m |)/|λ 1 |. 2. In the current paper we assume that resonances as do not occur, see (4.15) and (3.6). This allows us to analytically conjugate to linear flows e Λt and e Ωt on the slow stable manifold and normal bundle, respectively. When resonances occur, then a parametrization of the invariant normal bundle that analytically conjugate to a linear flow are in general no longer available. One either needs to work in less smooth classes of parametrizations, or conjugate to a more complicated flow. In particular, one may attempt to conjugate the flow on the normal bundle to a nonlinear flow in parameter space, along the lines of the analysis in [1]. In that case however, decoupling of the bundle into rank-1 invariant bundles is not available.

Numerical Computations for Fisher's Equation.
In this section, we illustrate our theoretical construction with a numerical implementation. Consider the partial differential equation subject to the Neumann boundary conditions where u n (t) are smooth functions of time. Proceeding formally, we obtain the system of ordinary differential equations The corresponding mapping F : R N +1 → R N +1 given by defines a real analytic (quadratic) vector field on R N +1 . With this notation the truncated differential equation is simply u = F (u).

Formal Computation of the Stable/Unstable
Manifolds. Throughout the following section we denote by u = (u 0 , . . . , u N ) a vector in R N +1 . Define the cosine convolution product truncated to order N to be the mapping [ * ] N : Note that the mapping is bilinear and commutative. Define also the diagonal linear operator L : R N +1 → R N +1 componentwise by the expression Then our vector field is rewritten as Note that the derivative of F at a acting on a vector h ∈ R N +1 is given by (4.10) Suppose now that p ∈ R N +1 is an equilibrium solution. For the eigenvalues of the linearization of F around p we use the notation as introduced in Section 2.1. We seek a function Note that Matching like powers for |α| ≥ 2 leads to for all |α| ≥ 2. Isolating from Equation (4.12) terms of order α and moving them to the left hand side of the equation leads to otherwise.
Comparing this with Equation (4.10) we see that Equation (4.13) can be rewritten as (4.14) We refer to Equation (4.14) as the homological equation for a. Note that each α the right-hand side of (4.14) depends only on terms of lower order. Hence we can solve recursively to any finite order as long as the non-resonance condition are satisfied for all multi-indices α and every 1 ≤ j ≤ m. Note that for m + 1 ≤ j ≤ n there is no possible solution of Equation (4.15). Note also that for |α| large enough there is no solution of Equation (4.15) and hence we have only a finite number of non-resonance conditions to satisfy. Remark 4.1 (Efficient Numerical Algorithms). We note that instead of solving Equation (4.14) term by term it is also possible to treat Equation (4.13) as a system of nonlinear equations and apply a Newton scheme. The scheme can be run with the linear approximation (4.11) as the initial guess. When a high order approximation of the manifold is desired an iterative scheme is often more efficient than solving recursively term by term due to the large number of multi-indices. The interested reader is referred to the numerical implementation at [24].

Formal Computation of the Invariant Frame Bundles.
In this section we continue to enforce the notation of Section 4.1. In particular, we assume that we have solved the homological equations (4.14) for the Fisher's equation to order K resulting in the approximate parameterization 16 where a α are numerical solutions of Equation (4.14). The essential term in the computation of the invariant bundles is the computation of composition term DF [a(φ)], see (3.1). Here for any u ∈ R N +1 let Y u denote the associated convolution matrix, i.e. since This observation allows us to write i.e., the (N +1)×(N +1) matrix coefficients in the Taylor expansion of DF [a(φ)] have the form A α = L − Y aα . This expansion facilitates the computation of the coefficients of the slow manifold Floquet normal form up to order K. We denote the resulting approximation by where the (N + 1) × (N + 1) matrices Q α are computed as discussed in Section 3. The interested reader should consult the computer codes at [24] for implementation details.

Numerical Computation of Some Invariant
Manifolds and Invariant Bundles. We now consider the results of a number of numerical computations.
Example 1: Projection Dimension N = 2 We begin by illustrating several computations in the N = 2 setting (so that the phase space is three dimensional). In this case the system reduces to The system has two equilibrium solutions which do not depend on the choice of κ. Note that for values of 1 < κ < 4 the origin has two unstable eigenvalues and one stable eigenvalue. (the eigenvalues of the trivial solution are κ − n 2 , n = 0, . . . , N ). We (arbitrarily) fix the value of κ = 2.423 (for readability we give only a few digits of all floating point numbers in this manuscript) and note that the system has another equilibrium solution of We (again somewhat arbitrarily) call λ s1 and λ s2 the slow and eigenvalue, respectively. We compute the full two dimensional stable manifold to order K = 25 and the slow stable manifold associated with the eigenpair (λ s1 , ξ s1 ) to order K slow = 35. This results in polynomials a K : R 2 → R 3 and a K slow : R → R 3 . We also compute the parameterization of the invariant frame bundle associated with a K slow and denote this by Q K . The frame bundle parameterized by Q K is illustrated in Figure 4.1.
In this low dimensional problem we can parameterize the two dimensional local stable manifold W s (p 1 ), the slow stable manifold W slow (p 1 ), and the invariant frame bundle of the slow manifold. It is then possible to check numerically that the stable normal bundle parameterized by the first column of Q K (φ) provides a good linear approximation of the full stable manifold in a neighborhood. To illustrate this we choose a point φ ∈ [−0.8, 0.8] and compute The results are plotted using log-log axes in Figure 4.2 for several values of φ as a function of |h|. We see that the error decreases quadratically as indicated by the fact that the lines all have slope of 2, and the magnitude of the error is rather uniform in φ. In each case the red vector corresponds to the stable normal vector bundle, blue to the unstable vector normal bundle, and green to the tangent bundle. We fix an sphere of radius r = 0.1 and advect the sphere for time T under the flow Φ. The resulting ellipse is plotted at φ 2 . Note that the unstable normal vector at a(φ 2 ) coincides most stretched axis of the ellipsoid, while the stable vector at a(φ 2 ) seems to coincide with the most contracted axis of the ellipsoid. This confirms our intuition that the flow is most contracting/expanding along the stable/unstable normal bundles respectively (even though this is only in fact only true up to quadratic approximation, i.e., the expansion/contracting rates hold exactly only for the variational flow).
In order to illustrate that the frame bundle parameterized by Q describes the most expanding and contracting directions normal to the local slow stable manifold W slow loc we choose φ ∈ B 0.8 (0) ⊂ R (the domain of a and Q) fix T, r > 0, let S r be the sphere of radius r centered at the origin in R 3 and compute We then advect S phase under the flow for time T and obtain S flow = Φ(S phase , T ). Now S flow is a set diffeomorphic to a sphere but stretched out in a nonlinear way by the flow. Note that since a(φ) is an invariant manifold the center of S flow is still on the slow stable manifold. If r and T are not too large then the action of the flow is approximated well by the action of the variational equations and we expect that S flow is ellipsoidal. In fact the expansion and contraction of the sphere is well aligned with the stable/unstable normal invariant bundles. This construction is illustrated numerically in Figure 4.3. The numerical accuracy of the conjugacy relation (and by virtue the accuracy of the expanding/contracting rate conditions) is discussed more quantitatively in the next section. The source code which produces these results is found in the file paperCode_lowDimExampleScript.m which is available at [24].

Example 2: Numerical Error in Higher Dimensions:
We now consider the accuracy of our computational schemes in various higher dimensions. Since we work in a higher dimensional phase space we consider parameterizations by two variable polynomials, i.e., we will consider the manifolds and invariant frame bundles defined by the two slowest eigenvalues (2D slow manifolds). These are computed for Fisher's equation truncated at various numbers of modes (i.e. with N varying).
As a measure of the accuracy of the computations we check the validity of the conjugacy equations under numerical integration. More precisely let a K be the K-th order polynomial approximation of the slow manifold associated with the two slowest eigenvalues. Then the quantity with fixed choice of T > 0 gives a measure of how well a K satisfies the flow conjugacy give by Equation on the disk of radius r. The flow Φ may be approximated by any numerical integration scheme we wish.
Similarly let Q K be the K-th order approximation of the invariant frame bundle of a, and for any φ let γ(t) be the solution of the initial value problem  Table 4.1 illustrates numerical accuracy of the parameterization of the two dimensional slow stable manifolds and bundles for a ten dimensional projection of the PDE, at several different values of κ. Note that the observed bundle conjugacy error is usually within three order of magnitude of the observed manifold conjugacy error. Table 4.2 on the other hand illustrates that for a given choice of scaling α and polynomial approximation order K, we can always adjust the size of the parameter domain r so that the Conjugacy errors are near machine precision. Table 4.3 illustrates the effects of increasing the order of polynomial approximation as the dimension of the system increases. Note that with other parameters fixed we are able to achieve almost machine precision in for both the manifolds and bundles by taking K sufficiently high.
Finally we remark that for all conjugacy tests we choose an integration time of T = 0.01. This is in order to focus our attention on the discretization errors associated with the manifold parameterization while minimizing the discretization error due to numerical integration. The interested reader can reproduce the results of this section using the program paperCode_dimensionComps.m which is available at [24].  Here the values of the system parameter κ, the order K of the polynomial approximation, and the scaling α of the slow eigenvectors are varied. The conjugacy relations are checked for one hundred points in a circle or radius r = 0.8 in the parameter domain. The conjugacy relations are integrated over the time interval [0, 0.01] using MatLab's standard ode45 routine.

Numerical Computation of Connecting Orbits.
In this section we discuss numerical computation of a connecting orbit for the truncated Fisher's equation given by Equation (4.5). More precisely suppose that N > 0 is fixed and that p 0 , p 1 ∈ R N +1 are equilibrium solutions of (4.5). Suppose that dim (W s (p 0 )) = m 0 , and that dim (W u (p 1 )) = N + 1 − m 0 + 1 ≡ m 1 , so that there is the possibility of a transverse heteroclinic connecting orbit from p 0 to p 1 . Then there are parameterizations of the local stable/unstable manifolds A 0 : B r0 (0) ⊂ R m0 → R N +1 and A 1 : B r1 (0) ⊂ R m1 → R N +1 having that image(A 0 ) = W u loc (p 0 ) and image(A 1 ) = W s loc (p 1 ).  We seek T > 0, θ ∈ B r0 and φ ∈ B r1 so that Φ A 0 (θ), T = A 1 (φ).
If (θ,φ,T ) is a zero of G then the orbit of A 0 (θ) (and equivalently the orbit of A 1 (φ)) is heteroclinic from p 0 to p 1 .
Such a solution is not isolated as any time shift of a heteroclinic orbit is again a heteroclinic orbit. We introduce spherical coordinates α in B r0 and β in B r1 parameterizing the m 0 − 1 dimensional sphere or radius r 0 and the m 1 − 1 dimensional sphere of radius m 1 − 1 respectively and define the constrained heteroclinic operator F : R N +1 → R N +1 given by F (α, β, T ) = Φ A 0 (α), T − A 1 (β), (4.18) and have that an isolated zero of F corresponds to a transverse heteroclinic connecting orbit from p 0 to p 1 . If we now discretize A 0 , A 1 and the flow Φ then we can evaluate the map F . By integrating the variational equations we compute the derivative of F and implement a Newton procedure to locate roots of F . In practice we compute the flow map Φ and the differential DΦ by numerically integrating the vector field and the variational equations with any reasonable scheme. On the left: the 2d unstable manifold of p 0 (blue), the 2d slow stable manifold of p 1 (red), the 1d unstable manifold of p 1 (blue), and the 2d slow stable manifold of p 2 all projected into three dimensions. The heteroclinic connection from p 0 to p 1 appears as a green arc. On the right: the first component of the connecting orbit in the time domain with the same color coding. Note that the parameterizations "swallow up" the flat portion of the curve.