The numerical computation of unstable manifolds for infinite dimensional dynamical systems by embedding techniques
Abstract.
In this work we extend the novel framework developed in [9] to the computation of finite dimensional unstable manifolds of infinite dimensional dynamical systems. To this end, we adapt a setoriented continuation technique developed in [10] for the computation of such objects of finite dimensional systems with the results obtained in [9]. We show how to implement this approach for the analysis of partial differential equations and illustrate its feasibility by computing unstable manifolds of the onedimensional KuramotoSivashinsky equation as well as for the MackeyGlass delay differential equation.
Adrian Ziessler, Michael Dellnitz and Raphael Gerlach
Department of Mathematics
Paderborn University
33095 Paderborn, Germany
1. Introduction
Setoriented numerical methods have been developed in the context of the numerical treatment of dynamical systems (e. g. [11, 12, 16, 18]). The basic idea of these methods is to cover the objects of interest – for instance invariant sets or global attractors – by outer approximations which are created via multilevel subdivision or continuation techniques. They have been used in several different application areas such as molecular dynamics ([42]), astrodynamics ([13]) or ocean dynamics ([17]).
The computation of (global) invariant manifolds has attracted a lot of interest in recent years. Approaches based on geometric concepts can be used to approximate (up to twodimensional) unstable manifolds of vector fields (see [33] for an overview). The approximation by geodesic level sets, for instance, produces a regular mesh that consists of geodesic circles by solving appropriate boundary value problems (e.g., [31, 32]). Topological methods are used to create outer approximations of invariant manifolds. In [33] a survey of different approaches for computing global stable or unstable manifolds of vector fields is given.
In this paper, we will focus on the setoriented continuation method that has been developed in [10] and we will show how to use it for the computation of unstable manifolds of infinite dimensional dynamical systems. We will, in particular, approximate unstable manifolds for semiflows of Banach spaces (cf. [20, 22, 3, 5]). Our approach relies on the results obtained in [9], where the subdivision technique developed in [11] has been extended to the infinite dimensional context. The underlying idea is to compute low dimensional invariant sets of infinite dimensional dynamical systems by utilizing embedding techniques for infinite dimensional systems [24, 40]. To this end, a delay embedding technique has been used in order to obtain a onetoone representation of the (infinite dimensional) dynamics by the socalled core dynamical system (CDS). The CDS is a continuous dynamical system on a state space of moderate dimension which we will refer to as the observation space. Then the subdivision scheme from [11] has been applied to the CDS in order to analyze onetoone images of global attractors for the infinite dimensional dynamical system. The feasibility of this approach has been illustrated by computing invariant sets for delay differential equations with constant time delay. Recently, this approach has also been generalized to delay differential equations with state dependent time delay [47]. In [2] lowdimensional transition manifolds of stochastic dynamical systems showing metastable behavior have been parameterized by using a combination of embedding and manifold learning techniques.
In addition to delay differential equations, other application scenarios in which finite dimensional invariant sets arise are certain types of dissipative dynamical systems described by partial differential equations, for instance the KuramotoSivashinsky equation [34, 44], the GinzburgLandau equation, or scalar reactiondiffusion equations with cubic nonlinearity [27]. For all these systems a finite dimensional socalled inertial manifold exists to which trajectories are attracted exponentially fast, e.g., [7, 15, 45].
In this paper we extend the classical setoriented continuation technique developed in [10] to the infinite dimensional context. This method allows us to approximate unstable manifolds of infinite dimensional dynamical systems in observation space. The general numerical approach we are proposing is in principle applicable to infinite dimensional dynamical systems described by a Lipschitz continuous operator on a Banach space. However, in this article we will focus on partial differential equations (PDEs). To this end, we will develop an appropriate numerical realization of the CDS for PDEs and illustrate the efficiency of our method for a onedimensional KuramotoSivashinsky equation. The KuramotoSivashinsky equation has attracted a lot of interest as a model for complex spatiotemporal dynamics and has been derived in the context of several extended physical problems, e.g., phase dynamics in reactiondiffusion systems [34] or small thermal diffusive instabilities in laminar flame fronts [44]. It is also a paradigm for the existence of complex finite dimensional dynamics in a PDE. In addition, we will also illustrate the efficiency of our method for the MackeyGlass delay differential equation by using the numerical realization of the CDS introduced in [9]. In the context of delay differential equations a related – though not setoriented – approach has recently been utilized in [19] for the numerical approximation of unstable manifolds. In this work the authors compute high order Taylor and FourierTaylor approximations of unstable manifolds for equilibria or periodic solutions.
A detailed outline of the paper is as follows. In Section 2 we briefly summarize the results of [9]. We state the main results of [24] and [40] and we describe the construction of the CDS on the observation space. In Section 3 we first briefly review the continuation method developed in [10] for the computation of unstable manifolds of finite dimensional dynamical systems. Then we explain how to use this technique in the context of infinite dimensional dynamical systems. In [10] convergence has been shown for compact subsets of unstable manifolds. Here we will extend this convergence result to the closure of the entire unstable manifold. In particular, we will prove that in this setting the embedding of the local unstable manifold (in infinite dimensional space) is identical to the local embedded unstable manifold (cf. Proposition 4 (b)). A numerical realization for the construction of the CDS for PDEs is given in Section 4. Finally, in Section 5, we illustrate the efficiency of our method for the MackeyGlass delay differential equation and for a onedimensional KuramotoSivashinsky equation. We will, in particular, generate numerical approximations of onetoone copies of unstable manifolds for the KuramotoSivashinsky equation for different parameter values.
2. Review of subdivision and embedding results
Since our setoriented continuation method is based on the framework developed in [9] we start with a short review of the related material.
2.1. Basic definitions and results from embedding theory
We consider dynamical systems of the form
(1) 
where is Lipschitz continuous on a Banach space . Moreover we assume that has an invariant compact set , that is
This assumption is justified by several classical results, where it has been shown that for many dissipative systems on Banach spaces there exist (nontrivial) compact attractors of finite capacity or Hausdorff dimension (see [38, 36, 8, 6, 21]). In order to approximate the set we combine a classical subdivision technique for the computation of such objects in a finite dimensional space with infinite dimensional embedding results (cf. [24, 40]). For the statement of the main result of [40] we need three particular notions: prevalence [41], upper box counting dimension and thickness exponent [24].
Definition 2.1.

A Borel subset of a normed linear space is prevalent if there is a finite dimensional subspace of (the ‘probe space’) such that for each belongs to for (Lebesgue) almost every .

Let be a Banach space, and let be compact. For , denote by the minimal number of balls of radius (in the norm of ) necessary to cover the set . Then
denotes the upper boxcounting dimension of .

Let be a Banach space, and let be compact. For , denote by the minimal dimension of all finite dimensional subspaces such that every point of lies within distance of ; if no such exists, . Then
is called the thickness exponent of in .
With these definitions the main results relevant for our framework are as follows:
Theorem 2.2 ([24]).
Let be a Banach space and compact, with upper box counting dimension and thickness exponent . Let be an integer, and let with
Then for a prevalent set of bounded linear maps there is such that
Note that this result implies that – if is large enough – a prevalent set of bounded linear maps will be onetoone on . This theorem lays the foundation for Robinson’s main result concerning delay embedding techniques:
Theorem 2.3 ([40]).
Let be a Banach space and a compact, invariant set, with upper box counting dimension , and thickness exponent . Choose an integer and suppose further that the set of periodic points of satisfies for . Then for a prevalent set of Lipschitz maps the observation map defined by
(2) 
is onetoone on .
Remark 1.

This result can be generalized to the case where several different observables are evaluated. More precisely, for a prevalent set of Lipschitz maps , , the observation map ,
is also onetoone on , provided that

As observed in [24], the thickness exponent is always bounded by the (upper) boxcounting dimension, i.e.,
Thus, providing that we know the upper boxcounting dimension of we can in principle always choose a worstcase embedding dimension such that
is satisfied (cf. Theorem 2.3).
2.2. The core dynamical system (CDS)
In [9] the results from Section 2.1 have been used to create a finite dimensional dynamical system that allows to approximate invariant sets for infinite dimensional dynamical systems on Banach spaces. In this section we will briefly review the construction of the CDS.
Let be a compact invariant set of the infinite dimensional dynamical system (1). We denote by the image of under the observation map , that is
(3) 
where according to Theorem 2.2 or according to Theorem 2.3. We then construct the CDS
(4) 
with as follows: At first we define on the set by
where is the continuous map satisfying
(5) 
In fact, this is possible since is invertible as a mapping from to . Observe that by this construction is an invariant set for . By using a generalization of Tietze’s extension theorem by Dugundji [14] we extend to a continuous map with (see Fig. 1).
Now we are in the position to extend the CDS to :
Proposition 1.
There is a continuous map satisfying
For the proof the reader is referred to [9].
2.3. Computation of embedded attractors via subdivision
By our construction, the dynamics of the CDS on is topologically conjugate to that of on . In what follows we will briefly review the main statements of [9]. We will use (4) in order to approximate the embedded invariant set (cf. (3)) via subdivision. To this end, we give a brief review of the related subdivision scheme.
Subdivision scheme: Let be a compact set. We define the global attractor relative to by
(6) 
The aim is to approximate this set with the subdivision procedure. By using the subdivision algorithm we obtain a sequence of finite collections of compact subsets of such that the diameter
converges to zero for . Given an initial collection , we inductively obtain from for in two steps.

Subdivision: Construct a new collection such that
(7) and
(8) where .

Selection: Define the new collection by
(9)
The first step guarantees that the collections contains successively finer sets for increasing . In fact, by construction
(10) 
In the second step we remove each subset whose preimage does neither intersect itself nor any other subset in .
Denote by the collection of compact subsets obtained after subdivision steps, that is
Observe that the ’s define a nested sequence of compact sets, that is, . Therefore, for each ,
(11) 
and we may view
(12) 
as the limit of the ’s. Then the selection step is responsible for the fact that the unions approach the relative global attractor:
Proposition 2.
Suppose that satisfies . Then
Observe that in contrast to [11] we have to assume that since the CDS is only continuous and not homeomorphic. Moreover, we note that we can, in general, not expect that . In fact, by construction may contain several invariant sets and related heteroclinic connections. However, if is an attracting set we can prove equality:
Proposition 3.

Let be the global attractor relative to the compact set , and suppose that the embedded attractor satisfies . Then
(13) 
Suppose that is an attracting set with basin of attraction and choose such that and . Define for the continuous maps
and denote the corresponding relative global attractors by , where
Then
where .
Remark 2.
Roughly speaking Proposition 3 (b) states that it is possible to approximate an attracting set for if we perform the computations with appropriately high iterates of .
3. A subdivision and continuation technique for embedded unstable manifolds
In this section we extend the results of [10] for the computation of finite dimensional unstable manifolds of infinite dimensional dynamical systems of the form (1). We will in particular focus on invariant manifolds for semiflows on Banach spaces (cf. [20, 22, 3, 5]).
3.1. Approximation of the local unstable manifold
Let us denote by
(14) 
the unstable manifold of , where is a steady state solution of the infinite dimensional dynamical system (cf. (1)). Furthermore, let us define the embedded unstable manifold by
(15) 
where and is the observation map introduced in Section 2. Observe that, by construction, is an invariant set for (cf. (5)). Our goal is to approximate compact subsets of or even the entire closure via an adaptation of a setoriented continuation method introduced in [10].
Let us denote by the local unstable manifold of the steady state and choose a compact neighborhood such that
Since is a homeomorphism on we can conclude that is compact. Then the following proposition holds:
Proposition 4.

Let be the global attractor relative to . Then

Let us suppose that is a compact attracting set with basin of attraction . Choose such that and . Then
Proof.

Define for the continuous maps
and denote the corresponding relative global attractors by , where
By Proposition 3 (b) we obtain
It remains to show that . Since , it is easy to see that
(cf. (5)). Thus,
∎
Remark 3.

Observe that Proposition 4 (b) states that the embedding of the local unstable manifold (in infinite dimensional space) is identical to the local embedded unstable manifold.

If the steady state is hyperbolic, then is attractive since by assumption its dimension is finite (cf. (14)). In particular, the compact set can be chosen, such that the assumed properties are satisfied.
From now on we assume that the assumptions of Proposition 4 (b) are satisfied. The idea of the continuation algorithm is to globalize the local covering of in order to obtain an approximation of the entire embedded unstable manifold .
3.2. The continuation method
The continuation starts at of the embedded unstable manifold . We choose a compact set containing and we assume that is large enough so that it contains the entire embedded unstable manifold of , i.e.,
(16) 
We remark that this assumption can be relaxed, and we will discuss this point later in the context of the realization of the approximation scheme.
In order to combine the subdivision process with a continuation method, we realize the subdivision using a family of partitions of . We define a partition of to be a finite family of compact subsets of such that
Moreover, we denote by the element of containing . We consider a nested sequence , of successively finer partitions of , requiring that for all there exist such that and for some . A set is said to be of level .
In what follows, we assume that for sufficiently large such that . The purpose is to approximate subsets where and
Here is the CDS (see (4)).
The numerical realization of the continuation algorithm for the approximation of embedded unstable manifolds can be described as follows:
Initialization: Given we choose an initial box , defined by a dimensional generalized rectangle of the form
where for , are the center and the radii, respectively. Choose a partition of and a box such that .

Apply the subdivision algorithm with subdivision steps to to obtain a covering of the local embedded unstable manifold .

Set

For define
(17)
Remark 4.

Observe that the unions
form a nested sequence in , i.e.,
In fact, it is also a nested sequence in , i.e.,

Due to the compactness of the continuation in Step 3 of Algorithm 1 will terminate after finitely many, say , steps. We denote the corresponding box covering obtained by the continuation method by
(18)
In Fig. 2 we illustrate the continuation method described in Algorithm 1.
Intuitively it is clear that the algorithm, as constructed, generates an approximation of the embedded unstable manifold . In particular, we expect that the bigger and are the better the approximation will be.
Proposition 5.

The sets are coverings of for all . Moreover, for fixed , we have

Suppose that is linearly attractive, i.e., there is a and a neighborhood such that
(19) Then the box coverings obtained by Algorithm 1 converge to the closure of the embedded unstable manifold . That is,
Proof.

Applying the subdivision algorithm with subdivision steps to the initial covering , we obtain of , that is,
Here we assume that the subdivision scheme is creating coverings using elements from the partitions . By Proposition 4 (b) and by Proposition 2 the box coverings converge to for . Therefore, for . Since is fixed a continuity argument shows that the sets converge to for , i.e.,

For each Algorithm 1 yields a covering of , and therefore
Suppose there is . Since is compact, it follows that . By definition of , Algorithm 1 generates a pseudo orbit for each , where . That is
Here denotes the element of containing and , i.e., for each continuation steps is covered. Observe that the sequence is monotonically increasing in (cf. Step 3 of Algorithm 1 and Remark 4) and
(20) We first show by contradiction that is unbounded. Suppose that is bounded by some , i.e., . Hence, by monotony of there is such that for all . Using Remark 4 (a) we have
However, by Proposition 5 (a) it follows that which is a contradiction to . Thus, is unbounded.
By assumption is linearly attractive in a neighborhood . Hence, we can use (19) and (20) on the pseudo orbit , where , in combination with the triangle inequality to obtain
Here the last expression converges to zero because and converges to zero for (see (10)). Again we have an contradiction to . It follows that
which yields the desired statement.
∎
Remark 5.

If (16) would not be satisfied, then it can in general not be guaranteed that the continuation method leads to an approximation of the entire set or even . Rather it has to be expected that this is not the case. The reason is that the embedded unstable manifold may ’leave’ but may as well ’wind back’ into it. In this scenario the continuation method, as described above, will not cover all of . This situation is illustrated in Fig. 3 (a).

The assumption in Proposition 5 (b) is, for instance, not satisfied if forms a heteroclinic connection between the steady state solution and another unstable hyperbolic steady state . In fact, in this case the algorithm would also generate a covering of the embedded unstable manifold of . This situation is illustrated in Fig. 3 (b).

If (19) is not satisfied, but is attractive, one can apply the subdivision scheme introduced in [9] to in order to approximate more accurately (cf. Proposition 3 (b)).

Note that (19) is satisfied if the observation map is biLipschitz with Lipschitz constant and is linearly attractive with .
4. Numerical realization of the CDS for partial differential equations
As discussed in the introduction, dynamical systems with infinite dimensional state space, but finite dimensional attractors arise in particular in two areas of applied mathematics, namely dissipative partial differential equations and delay differential equations with small constant time delay. In this article we will focus on the PDE case, and we will present one specific realization of the maps and for this situation.
More precisely, we will consider explicit differential equations of the form
(21) 
where is in some Banach space and is a (nonlinear) differential operator. We assume that the dynamical system has a welldefined semiflow on .
In order to numerically realize the construction of the map described in Section 2.2, we have to work on three tasks: the implementation of , the implementation of , and the realization of the timemap of (21), denoted by , respectively. For the latter we will rely on standard methods for forward time integration of PDEs, e.g., a fourthorder time stepping method for the onedimensional KuramotoSivashinsky equation [29]. Observe that the numerical realization of strongly depends on the underlying PDE. The map will be realized on the basis of Theorem 2.2 or Remark 1, respectively. For the numerical construction of the continuous map we will present a new method that uses statistical information from previous computations. This step is in particular crucial for the continuation step, since we want to restart the algorithm with initial conditions that satisfy the identities
(cf. (5)) at least approximately.
From now on we assume that upper bounds for both the box counting dimension and the thickness exponent are available. This allows us to fix according to Remark 1.
4.1. Numerical realization of
In [9] has been defined on the basis of Theorem 2.3, where function evaluations at a fixed time have been used. Thus, in the case of scalar delay differential equations has been defined as the delay coordinate map
(22) 
where is the constant timedelay of the underlying DDE. In principle it would also be possible to observe the evolution of a partial differential equation by a delay coordinate map. However, from a computational point of view this would be very inefficient. The reason is that for the realization of the map one would have to reconstruct functions from time delay coordinates. Thus, for each point in observation space one would essentially have to store the entire corresponding function. To overcome this problem, in this work we will present a different approach. In what follows, we will assume that the function can be represented in terms of an orthonormal basis , i.e.,
(23) 
where the are elements from a Hilbert space (e.g., ). Then our observation map will be defined by projecting a function onto coefficients of its Galerkin expanion. For the approximation of , i.e.,
(24) 
we want to use an optimal basis (i.e., as small as possible) in the sense that it contains the ’most characteristic’ data from an ensemble of functions. The notion ’most characteristic’ implies the use of an averaging operation. Furthermore, this basis has to be capable of representing the solution of the underlying PDE (21) with a small error. Both requirements can be addressed by the proper orthogonal decomposition (POD) (c.f. [43, 1, 23]), also known as the principal component analysis or the KarhunenLoève transformation. In order to compute a basis for and , we first generate timesnapshots of a longtime simulation for some set of initial conditions of the underlying PDE. Then we approximate the PODbasis via the singular value decomposition (e.g., [4, 35, 46]). Using this basis we then approximate by (24) where denotes the th PODcoefficient at time .
Given the basis and using the fact that this basis is orthogonal, we then define the observation map by choosing different observables
(25) 
This yields
(26) 
Observe that is linear and bounded and hence, for sufficiently large, Theorem 2.2 and Remark 1, respectively, guarantee that generically (in the sense of prevalence) will be a onetoone map on .
4.2. Numerical realization of
In the application of the continuation scheme for the computation of embedded unstable manifolds described in Section 3 one has to perform the continuation step
(see (17)). Numerically this is realized as follows: At first is evaluated for a large number of test points for each box . Then a box is added to the collection if there is a least one such that .
Remark 6.
In practice the test points can be chosen according to several different strategies: In low dimensional problems one can choose them from a regular grid within each box . Alternatively one can select the test points from the boundaries of the boxes. In our computations we use a Monte Carlo sampling.
By Section 4.1 the state space for the CDS is given by points where are the PODcoefficients. For the evaluation of at a test point we need to define the image , that is, we need to generate adequate initial conditions for the forward integration of the PDE (21). In the first step of the continuation method we proceed as follows. Given the PODbasis for and , we simply construct initial conditions near the unstable (hyperbolic) steady state by defining the map as
(27) 
Observe that by this choice of and both conditions and are satisfied for each test point (see (5) and (26)). Then we use the timemap of the underlying PDE to obtain a function . If is not sufficiently large, then will in general not be satisfied anymore. Therefore, it is possible that initial conditions for , , generated by (27) are not even close to the unstable manifold. This is not acceptable since the requirement for all (see (5)) is crucial in order to compute a reliable covering of the embedded unstable manifold.
To enforce this equality at least approximately we extend the expansion and construct initial functions by
(28) 
Here only the first PODcoefficients are given by the coordinates of points inside . Thus, it remains to discuss how to choose the PODcoefficients . The idea is to use a new heuristic strategy that utilizes statistical information obtained in the previous continuation step: Suppose we want to evaluate for a large number of test points in a box . By the continuation step (cf. Algorithm 1), there must have been at least one such that for at least one test point . For all these points we can compute the PODcoefficients by
Then we sample the box with all points for which additional information is available. However, the number of these points might be too small, such that is not discretized sufficiently well and we have to generate additional test points. For this, we first choose a certain number of points at random. Then we extend these points to elements in as follows: We first compute componentwise the mean value and the variance of all PODcoefficients , for . This allows us to make a Monte Carlo sampling for the additional coefficients of for , i.e.,
Finally, we compute initial functions of the form
By this construction we expect in each continuation step to generate initial functions that satisfy an approximation of the identity for all . We will illustrate our statistical approach in the next section for the onedimensional KuramotoSivashinsky equation.
5. Numerical results
In this section we present results of computations carried out for the MackeyGlass delay differential equation and the KuramotoSivashinsky equation, respectively. For the numerical realization of the CDS for delay differential equations the reader is referred to [9].
5.1. The KuramotoSivashinsky equation
We start with the wellknown
KuramotoSivashinsky equation in one spatial dimension which is given by
(29)  
This equation has been studied extensively over the past years. It has, for instance, been used to model phase dynamics in reactiondiffusion systems [34] or small thermal diffusive instabilities in laminar flame fronts [44]. Following [25, 30], we normalize the KS equation to an interval length of and set the damping parameter to the original value derived by Sivashinsky, i.e., . Then equation (29) can be written as
(30)  
In equation (30) we introduce a new parameter , where denotes the size of a typical pattern scale (cf. (29)). In [25, 30] numerical and analytical studies were made by varying over a finite interval, showing the complex hierarchy of bifurcations. We are in particular interested in computing the unstable manifold of the trivial unstable steady state for different parameter values . In order to use our algorithm developed in Section 3 it is crucial to have a good estimate of the dimension of the invariant set and , respectively (cf. Section 2). In [39] it has been shown that the dimension of the inertial manifold of (29) for is , i.e., each invariant set has finite dimension. However, these estimates are very pessimistic and we expect that we will obtain onetoone images of the unstable manifold for smaller related embedding dimensions .
In what follows, the observation space is defined through projections onto the first PODcoefficients. For each parameter value we compute the PODbasis by using the snapshotmatrix obtained through a longtime integration with the initial condition
(cf. Section 4.1).
5.1.1. The traveling wave
For the parameter value the KuramotoSivashinsky equation possesses a stable traveling wave solution (cf. Fig. 4 (a)). Due to the symmetry imposed by the periodic boundary conditions there are two waves traveling in opposite directions [30]. Correspondingly the CDS possesses for each traveling wave a limit cycle in observation space (cf. Fig. 4 (b)). We expect that the dimension of the embedded unstable manifold is approximately two since different initial conditions result in trajectories in observation space that are rotations of each other about the origin.
By choosing the embedding dimension we restrict the initial functions in the first continuation step to the subspace that is spanned by the first seven PODmodes and since , we expect to obtain a onetoone image of . In Fig. 5 we show a comparison of the numerical realization of the continuous map with initial functions generated by (27), i.e., where for all test points (red) and the statistical approach discussed in Section 4.2 (blue). By using only PODcoefficients, we construct initial functions that by far do not satisfy .
We choose and initialize a fine partition of for . Next we set . In addition, we define a finite time grid , where , and mark all boxes that are hit in each time step (a similar approach has been used in [28]). This strategy will be used for each example in this section.
In Fig. 6 (a)(d) we illustrate successively finer box coverings of the unstable manifold as well as a transparent box covering depicting the complex internal structure of the unstable manifold. Observe that the boundary of the unstable manifold consists of two limit cycles which are symmetric in the first PODcoefficient . This is due to the fact that the KuramotoSivashinsky equation with periodic boundary conditions (30) possesses symmetry.