Discretization strategies for computing Conley indices and Morse decompositions of flows

Conley indices and Morse decompositions of flows can be found by using algorithms which rigorously analyze discrete dynamical systems. This usually involves integrating a time discretization of the flow using interval arithmetic. We compare the old idea of fixing a time step as a parameters to a time step continuously varying in phase space. We present an example where this second strategy necessarily yields better numerical outputs and prove that our outputs yield a valid Morse decomposition of the given flow.


INTRODUCTION
While the numerical approximation of ordinary differential equations has a long history, the systematic study of how to compute time invariant structures began in the 1980s. Rigorous computations of these structures is an even more recent phenomenon (see [1,8] and references therein). These latter efforts can be roughly divided into two approaches: direct computation of invariant sets, e.g., periodic orbits, heteroclinic and homoclinic orbits, invariant manifolds, and a more indirect approach based on identification of isolating neighborhoods. The advantages of the direct approach are clear: it leads to efficient algorithms and provides precise numerical bounds on solutions. The disadvantage is that, in general, it does not provide information about the global structure of the dynamics and the results are sensitive to changes in parameters. As indicated below, the indirect approach overcomes these disadvantages, however, the construction of appropriate methods of approximation has proven to be a significant challenge.
For the sake of simplicity consider a smooth ordinary differential equation that generates a flow ϕ : R × R d → R d . The focus of this paper is on obtaining a finite representation of ϕ from which information concerning the structure of the dynamics can be obtained. A theoretical resolution to this challenge is as follows. To perform computations, we make use of a discretization in time. Choose h > 0 and define the diffeomorphism ϕ h : R d → R d by In analogy with the setting of flows, the maximal invariant set in N under ϕ h is defined by Inv(N, ϕ h ) := {x ∈ N | ϕ n h (x) ∈ N, ∀n ∈ Z} from which one can define isolating neighborhoods and isolated invariant sets.
A fundamental result [14,Theorem 1] guarantees that S ⊂ R d is an isolated invariant set for ϕ if and only if S is an isolated invariant set for ϕ h , independent of the choice of h > 0. Furthermore, according to [14,Theorem 2] the Conley index of S computed using ϕ h determines the Conley index for S under ϕ. Thus from a mathematical perspective no information is lost by studying the map ϕ h as opposed to the flow ϕ.
To obtain a representation of ϕ h with which we can compute, we discretize the phase space. For the sake of simplicity of presentation, we assume that we are interested in understanding the global dynamics of ϕ or equivalently ϕ h restricted to a rectangular isolating neighborhood X ⊂ R d . This allows us to discretize X using a cubical grid X (see Section 2) with diameter less than some given ε > 0. The dynamics of ϕ h can be encoded using a multivalued map F ρ : Making use of the fact that F ρ can be interpreted as a directed graph, a wide variety of efficient algorithms have been developed for finding isolating neighborhoods and Morse decompositions for ϕ h (see [2], [11]) and for computing the associated Conley indices [10,15,21]. Furthermore, by choosing finer grids, i.e., letting ε → 0 and making better approximations, i.e., letting ρ → 0, one can recover any isolated invariant set [11]. More generally, one can find all attractor-repeller pairs or equivalently all Morse decompositions and in the limit recover Conley's fundamental decomposition theorem [3,11]. Thus, from a theoretical perspective this approach allows us to recover the local and global dynamics associated with isolated invariant sets generated by ordinary differential equations.
However, in concrete applications there are significant technical obstructions. Observe that there are three explicit parameters in the above mentioned approach: the time step h, the grid size ε, and the accuracy of the approximation of the dynamics ρ. Furthermore, the dynamics that can be extracted is quite sensitive to the particular choices of these parameters. For example, for a fixed ε, if h is too small then ξ ∈ F ρ (ξ ) for every ξ ∈ X and every ρ ≥ 0. In this case no isolating neighborhood can be extracted and hence there is no interesting dynamics that can be resolved. Notice that the number of grid elements increases rapidly as an inverse function of ε and thus it becomes computationally prohibitive to choose ε too small. Observe that if v(x) = 0 for all x ∈ ξ , then it is reasonable to assume that for a sufficiently large h, ξ ∈ F ρ (ξ ), at which point one may hope to be able to start extracting interesting dynamical features. Clearly, the size of the parameter h depends on v , the magnitude of the vector field. The naive conclusion is that h should be chosen to be large. However, this leads to other difficulties. We need to evaluate ϕ h which in practice requires us to integrate the differential equation (1) over time h. Simple Gronwall inequalities imply that, in general, the numerical bounds on errors associated with this integration will grow exponentially in h. Suggesting that large h leads to exponentially large ρ, which in turn leads to large images of F ρ . Larger images of F ρ imply that less information related to the dynamics of ϕ can be extracted.
To address these issues, in this paper we consider a more general approach to encoding the dynamics. Let τ : X → (0, ∞) be a continuous function and consider the map ϕ τ : X → R d defined by ϕ τ (x) := ϕ(τ(x), x). In principle this allows us to vary the integration time over the phase space, which allows us to compensate for the varying magnitude of the vector field and sensitivity of propagation of numerical errors. A brief outline of the paper is as follows.
Section 2 presents numerical examples that contrast results using fixed time steps ϕ h and spatially varying time steps ϕ τ . The goal is not to provide new results concerning the dynamics of a particular ordinary differential equation, but rather to demonstrate the potential usefulness of this approach.
We use the results of Section 2 as motivation for the central mathematical results of this paper. As indicated above, [14,Theorem 1] guarantees that on the level of isolated invariant sets the dynamics of ϕ h captures the dynamics of ϕ. Since this is not always true for a ϕ τ with varying time steps we prove Lemma 3.18, which provides a criterion under which we obtain the desired isolation.
In Section 3.2, we show that computing the Conley index of an isolated invariant set S for ϕ τ and ϕ provides the Conley index for S under the flow ϕ. This generalizes [14,Theorem 2]. In Section 3.3, we prove that a Morse decomposition for ϕ τ is a Morse decomposition for the flow ϕ if all the Morse sets are invariant under ϕ, thereby showing that our algorithm does compute what we are interested in.
Varying the time step in space was first proposed in [7], but with less attention for the mathematical background, which we present here.
This article uses the following notations: R ≥0 := {x ∈ R | x ≥ 0}, R + := {x ∈ R | x > 0}. We recall definitions where they are needed for our statements. Their proofs require some background in dynamical systems, which can be found in [22] and other textbooks about the field.

A NUMERICAL EXAMPLE
2.1. Combinatorial representations. Given d > 0 and a tuple s ∈ R d + representing a box size, the space R d is the union of the following cubes (more precisely, products of intervals): Let X ⊂ Y be a finite subset. For a set A ⊂ X , its realization is |A| : Our goal is to compute a meaningful Morse decomposition within a compact set X := |X |.
A map F from X to its power set is written as F : X ⇒ X . We consider it as a multivalued map on X or as a directed graph which has an edge from ξ ∈ X to ξ ∈ X iff ξ ∈ F(ξ ). The image of A ⊂ X is F(A) := ξ ∈A F(ξ ) ⊂ X . For a subset Z ⊂ X, we let int X (Z) be the interior of Z with respect to the subspace topology of X ⊂ R d .
We apply the algorithms and software described in [2] and [4] for finding Morse decompositions (Definition 3.17) as follows. The software constructs a combinatorial enclosure F : X ⇒ X of a discrete dynamical system f : R d → R d . Then it constructs the strongly connected path components {N (p) | p ∈ P} of the directed graph F and builds the directed graph MG(F) with vertices P and a directed edge from p to q if there are G ∈ N (p), H ∈ N (q) and a path from G to H in F. This describes the dynamics of the discrete dynamical system f because of the following Theorem 4.1 from [11]. Theorem 2.3. Let X = |X | be an isolating neighborhood for f and S = Inv(X, f ). Then: In the rest of this section, we give an example flow ϕ and show that this algorithm yields a finer output when using f (x) = ϕ(τ(x), x) than when using f (x) = ϕ(h, x). The justification that our outputs are indeed Morse decompositions for ϕ is formulated in Theorem 3.20. Every norm . in this article is Euclidean, i.e., Fixed time step. The following ordinary differential equation is particularly challenging to analyze with a fixed time step h.
The equation has a fixed point (0, 0) and limit cycles with radius 1 and √ µ around the fixed point. This can be seen by its representation in polar coordinates: The norm of the vector field v increases quickly away from the origin because Far away from the origin, the solutions behave like the solutions ofṙ(t) = r(t) 5 . This equation is solved by where r 0 = r(0). Note that the function r is only defined for t < ( Hence, also our integration algorithm fails when the input is a box If h is chosen larger, each box ξ ∈ X near the boundary of X is assigned an arrow in F h to all of the other boxes because the algorithm fails to find an enclosure of ϕ h (ξ ) (which does not even exist). This would lead to a very large invariant part Inv(X , F h ) and should therefore be avoided. But when choosing h small enough, the directed graph F h has a lot of small strongly connected components whose corresponding invariant set of ϕ is empty (so-called spurious Morse sets). They have to occur when the time step h is so small that ϕ(h, ξ ) ∩ ξ = ∅. But this happens easily near the origin where v(x) is small. The outputs for µ = 2 are shown in Figure 1(b). The algorithm as described in [4] subdivided each dimension of |X | into 2 9 intervals of equal length. Each colored region is a combinatorial Morse set N (p), p ∈ P. The index set P had 57 675 elements. All but two of the combinatorial Morse sets found are empty.
It took 27 seconds on a laptop with an Intel i5 CPU to find the combinatorial Morse sets, but the Morse graph could not be computed because P was too large for the memory. Also subdividing each dimension into 2 12 intervals did give similarly bad outputs (after 2396 seconds).

2.3.
Variable time step. We use the following heuristic for the time step function τ. For each subdivision level, the Euclidean norm s is the diagonal of each of the congruent boxes ξ ∈ X . Using parameters D > 1 and δ > 0, we define the continuous function (4) τ :

FIGURE 2.
Output for the same system as in Figure 1, but using time step function τ from Equation (4) as proposed in Section 2.3.
The idea is that the distance between x and ϕ(τ(x), x) should be around D box diagonals -using a first-order approximation ϕ(τ(x), Since τ is usually not constant within a box ξ , the value τ(ξ ) is an interval. We can find an enclosure of this interval by considering τ as a function on intervals and replacing x by the box ξ when calculating in interval arithmetic. The software library CAPD [5] is used to construct a combinatorial enclosure F τ : X ⇒ X of ϕ τ . We use D = 4 and δ = 0.1 in the following numerical example. The varying time step strategy proposed above yields the finest Morse decomposition of the flow in S = Inv([−3, 3] 2 , ϕ) = { x | x ≤ √ µ } (more precisely, the finest one which does not contain empty invariant sets). The output is shown in Figure 2.
Finding the Morse decomposition and the Morse graph took 29 seconds on the same hardware as in the first example, but with each dimension subdivided into only 2 8 intervals. The algorithm additionally checks the correctness of the computed Morse decomposition. This verification and its output in Figure 2(c) are described in Remark 3.21.

THEORETICAL JUSTIFICATION
For the remainder of the paper we let Y be a locally compact separable metric space. Let ϕ : R × Y → Y be a flow. We show when certain invariants for ϕ τ also yield the corresponding invariants for ϕ.
Even though there is an inverse ϕ −h of ϕ h for any flow ϕ, there need not be an inverse for ϕ τ .
Theorem 3.3. Let S ⊂ Y be compact. Then the following three conditions are equivalent: (i) S is an isolated invariant set with respect to ϕ.
(ii) For every continuous map τ : Y → R + , S is an isolated invariant set with respect to ϕ τ . (iii) There is a number h > 0 such that S is an isolated invariant set with respect to ϕ h .
Proof. Assume condition (i) holds and fix a continuous map τ : Y → R + . Choose N, an isolating neighborhood for S with respect to ϕ. Obviously ϕ τ (S) ⊂ S. Using Lemma 3.2, for every x ∈ S there is an x ∈ S such that ϕ τ (x ) = x. This means S ⊂ ϕ τ (S). Hence S is invariant with respect to ϕ τ .
To see that S is an isolated invariant set with respect to ϕ τ , consider the map One easily verifies that each x ∈ S has a compact neighborhood V x such that σ (V To show the opposite inclusion take x ∈ Inv(M, ϕ τ ) and let γ : Z → Y be a solution of ϕ τ through x. Let x n := γ(n) and t n := τ(x n ). Then τ (x) = ϕ(t n , ϕ n τ (x)) = ϕ(t n , x n ). By definition of T , we have t n ≤ T . Since x n ∈ M, we have σ (x n ) ≥ T . It follows that ϕ([0,t n ], x n ) ⊂ N for all n, and consequently x ∈ Inv(N, ϕ) = S. Thus implication (i) =⇒ (ii) is proven. Implication (ii) =⇒ (iii) is obvious because we can always take τ to be a constant positive function. Implication (iii) =⇒ (i) is part of Theorem 1 in [14].
The following follows immediately from the implication (iii) =⇒ (i).  There is no full analogue of Corollary 3.4 since we cannot replace condition (iii) in the theorem above by the statement (iii') There is a continuous function τ : X → R + such that S is isolated invariant with respect to ϕ τ . An example for which (iii ) =⇒ (i) is wrong can be constructed by considering a system with a limit cycle and using a time step function τ which sends a certain point x on the limit cycle to itself under ϕ τ by letting τ(x ) be the period of the orbit. With an appropriate choice of τ for x near x , this yields an isolated fixed point x for ϕ τ , but the set S = {x } would not be invariant for ϕ. We defer a detailed discussion of this example to Remark 3.22, where it is discussed in the context of Morse decompositions.

Conley indices.
Given an isolated invariant set of a flow or of a discrete dynamical system, its Conley index is well-defined. Even though we are interested in the Conley index of a flow (as described in [6]), we propose using algorithms developed for the calculation of Conley indices of maps. Therefore, we need to consider both situations here. Definition 3.6. Let N ⊂ Y be compact and S = Inv(N, ϕ). A pair (P 1 , P 2 ) of compact sets P 2 ⊂ P 1 ⊂ N is an index pair for (S, ϕ) in N if it has the following properties: For an isolated invariant set S, an index pair exists and the following definition does not depend on a particular choice of the isolating neighborhood N and the index pair (P 1 , P 2 ). There are several ways of defining the Conley index for maps. Each definition uses a certain kind of index pair and and equivalence relation on a map on this index pair (e.g., [19,9,16]  A pair (P 1 , P 2 ) of compact sets Let H * denote the homology functor with compact supports as defined in [12,Chapter 9]. We use this version of the homology functor because of its good excision properties. However, note that this homology theory coincides with singular homology on compact neighborhood retracts, in particular compact polyhedra and cubical sets.
Let (P 1 , P 2 ) be an index pair for (S, f ) in N. Define maps on topological pairs: N)).
The induced homomorphism H * (i p ) in relative homology is an isomorphism because of excision (see, [12,Theorem 9.3]). The endomorphism In the next two subsections, we use the following notations for an endomorphism α : V → V . Let dom(α) = V denote the domain of α and gker(α) the generalized kernel defined by 3.2.1. Leray functor. Given an endomorphism α : V → V of graded modules (e.g. homology groups), the Leray functor L assigns to α an automorphism L(α) as follows. Let

Then let
Observe that if α is already an automorphism, then L(α) = α. For a detailed definition of this functor, we refer the reader to [16,Section 4]. We define α ∼ β iff L(α) and L(β ) are conjugate.
Here we let the Conley index be the equivalence class of L(I P ) under this equivalence relation. This definition only depends on S and not on the choice of P or N. [16,Theorem 2.6]. For all the index maps I p considered in this article, the automorphism L(I P ) is conjugate to the identity as is shown later in this section.
The main ingredient for comparing the Conley index of ϕ with the Conley index of ϕ h for some h > 0 is the existence of a common index pair for both dynamical systems. The following lemma was verified in the proof of [14,Theorem 2]. Lemma 3.9. Let S be an isolated invariant set for ϕ h (hence for ϕ). Then there is a pair Q = (Q 1 , Q 2 ) of compact sets such that Q is an index pair for (S, ϕ) and an index pair for (S, ϕ h ).
We also have a slightly more general statement for a time step function τ, but we need an extra assumption about the equality of the invariant sets with respect to ϕ and ϕ τ .  λ (ϕ(t, x)).
For an arbitrary ζ > 0, define subsets of V : The local compactness of Y can be used to observe that for any open neighborhood U of S there is a ζ > 0 such that ϕ([0, T ], H(ζ )) ⊂ U. Applying this observation to U = intV , we conclude the existence of an ε > 0 such that ϕ([0, T ], H(ε)) ⊂ V .
Let M := H(ε). Applying the observation to U = G(ε) shows the existence of a δ > 0 such that ϕ([0, T ], H(δ )) ⊂ G(ε). Define The pair (Q 1 , Q 2 ) is an index pair for (S, ϕ) and for (S, ϕ τ ) in M. It is straightforward to check the properties in Definitions 3.6 and 3.8, hence (i) is proven. To see (ii), we note that is a homotopy from i Q to f Q . Hence the index map I Q is the identity.

Shift equivalences.
Another equivalence relation commonly used on the index map I P is the following one proposed in [19]. If P and P are two index pairs around the same isolated invariant set for a discrete dynamical system, then I P and I P are shift equivalent.
The following theorem is similar to, but slightly stronger than Theorem 3.11.
Theorem 3.14. Let S be an isolated invariant set for ϕ and ϕ τ . Let (P 1 , P 2 ) be an index pair for (S, ϕ τ ) and I P its index map. Then there is an isomorphism CH * (S, ϕ) ∼ = H * (P 1 , P 2 )/ gker(I P ).
Proof. There is a common index pair (Q 1 , Q 2 ) for (S, ϕ) and (S, ϕ τ ) due to Lemma 3.10. Thus, by Theorem 3.13, there are homomorphisms r, s and an n ∈ N such that the following diagram commutes.
Definition 3.16. Let f : Y → Y be a discrete dynamical system. Let y ∈ Y be such that there is a solution Z → Y through y.
Definition 3.17. Given a dynamical system (i.e., a flow ϕ or a map f ) and an isolating neighborhood X with S = Inv(X) (denoting Inv(X, ϕ) or Inv(X, f ), respectively), we call a set of disjoint isolated invariant sets {M p | p ∈ P} together with an acyclic directed graph MG with vertices P a Morse decomposition in X if for every y ∈ S one of the following holds: (i) y ∈ M p for a certain p ∈ P; or (ii) there are p, q ∈ P and a path from p to q in MG, such that α(y) ⊂ M p and ω(y) ⊂ M q .
In order to show that a certain Morse decomposition for ϕ τ is also a Morse decomposition for ϕ, we can apply the following two lemmata. For Lemma 3.18. Let {M p | p ∈ P} be a Morse decomposition in X for ϕ τ with isolating neighborhoods N p , i.e., N p ∩ N q = ∅ for p = q and M p = Inv(N p , ϕ τ ) ⊂ int N p for all p ∈ P. Let p ∈ P. If ( * ) ϕ [0,τ] (N p ) ⊂ X and ϕ [0,τ] (N p ) ∩ N q = ∅ for all q ∈ P \ {p}, then Inv(N p , ϕ) = Inv(N p , ϕ τ ).
Let γ : Z → N p be a solution to ϕ τ in N p . Then for each n ∈ Z: ϕ([0, τ(γ(n))], γ(n)) ⊂ N p by assumption ( * ) in the lemma. Gluing these pieces together yields a trajectory for ϕ in N p . We also show Assume there is a point x ∈ Inv(N p , ϕ τ ) \ Inv(N p , ϕ τ ). Hence, x / ∈ ∪ q∈P M q . Now either α(x, ϕ τ ) or ω(x, ϕ τ ) must lie in some M q with q = p because they cannot lie within the same Morse set. But this means that any solution of ϕ τ through x has to contain points in int N q , which is disjoint from N p . We conclude x / ∈ Inv(N p , ϕ τ ), a contradiction. Overall, we get the following inclusions, where the middle one is trivial.
Each set is the same subset of N p . Therefore also Inv(N p , ϕ) = Inv(N p , ϕ). Lemma 3.19. Let X ⊂ Y be an isolating neighborhood for ϕ τ (and hence for ϕ). Let {M p | p ∈ P} be a Morse decomposition for ϕ τ in X with Morse graph MG and assume that each M p is invariant also for ϕ. Then {M p | p ∈ P} is also a Morse decomposition for ϕ in X with Morse graph MG.
Proof. We only need to show that the Morse graph MG is preserved. From here, consider a point y ∈ Inv(X, ϕ) ⊂ Inv(X, ϕ τ ). There are p, q ∈ P and a directed path from p to q in MG such that α(y, ϕ τ ) ⊂ M p and ω(y, ϕ τ ) ⊂ M q . We show that ω(y, ϕ) ⊂ M q by contradiction. The analogous statement for the α-limit is proven similarly.
Let N q be an isolating neighborhood for ϕ τ and therefore for ϕ around M q , hence Inv(N q , ϕ) = Inv(N q , ϕ τ ) = M q ⊂ int N q . We continue analogously to the proof of Theorem 3.3. We define a function σ : . Now there is a compact neighborhood N q of M q such that σ (x) ≥ T for all x ∈ N q . For all n ∈ N, let γ(n) := ϕ n τ (y). Assume that ω(y, ϕ) contains a point y outside of N q . We construct a subsequence of γ as follows: Since N q is a neighborhood of ω(y, ϕ τ ) and γ(N) ∩ M q = ∅, there is anñ ≥ 0 such that γ(ñ) ∈ N q \ M q . Let s > 0 be such that γ(ñ) = ϕ(s, y). There is an s ≥ s such that σ (ϕ(s , y)) = T , and then ϕ([s , s + T ], x) ⊂ cl(N q \ N q ). Hence, there is an n 0 ≥ñ such that γ(n 0 ) ∈ cl(N q \ N q ).
Going on like this, we construct a subsequence γ(n 0 ), γ(n 1 ), . . . of points in cl(N q \ N q ). It has a converging subsequence whose limit is also in ω(y, ϕ τ ) ⊂ M q . A contradiction.
We are ready to show  In Figure 2(c), each set N (p) is shown with a darker collar around it such that together they form Z(p). In our example, the algorithm successfully checked Z(p) ⊂ X and Z(p) ∩ N (q) = ∅ for any p = q. The additional checks for criterion (B) took only about a second.  The following example shows that when (A) is not fulfilled, we need some kind of check that a Morse set M p for ϕ τ is also invariant under ϕ. For every point θ ∈ S 1 , its limit sets are ω(θ , ϕ) = S 1 and α(θ , ϕ) = S 1 . The sets ∅ and S 1 are the only subsets invariant for ϕ. Define the time step function τ : S 1 θ → sin θ + 2π ∈ R + , which is well-defined and continuous since it has period 2π. Figure 3 shows a plot of τ and a trajectory of ϕ τ .

FINAL REMARKS
Using a time step function varying in space gives a lot of flexibility when applying rigorous algorithms for finding Morse decompositions. Sometimes, the step function τ is the only way to obtain a meaningful Morse decomposition numerically. Additionally, the heuristic presented can be justified more naturally than a specific choice of fixed time step. Up to now, we cannot see from a given equation which strategy is more promising. It would be interesting to find a more systematic approach to find out when to use which strategy.