EFFICIENT COMPUTATION OF LYAPUNOV FUNCTIONS FOR MORSE DECOMPOSITIONS

We present an efficient algorithm for constructing piecewise constant Lyapunov functions for dynamics generated by a continuous nonlinear map defined on a compact metric space. We provide a memory efficient data structure for storing nonuniform grids on which the Lyapunov function is defined and give bounds on the complexity of the algorithm for both time and memory. We prove that if the diameters of the grid elements go to zero, then the sequence of piecewise constant Lyapunov functions generated by our algorithm converge to a continuous Lyapunov function for the dynamics generated the nonlinear map. We conclude by applying these techniques to two problems from population biology.

Observe that R(f ) ⊂ M A .From the perspective of dynamics the most significant difference between a finest Lyapunov function V of the form defined in (1) and a function V A of the form defined in (2) is that the latter may be constant on a larger set of orbits than V .Furthermore, if one chooses a nested sequence of sublattices that converges to Att(X, f ), then V A will converge monotonically to V .
This suggests that we approximate a finest Lyapunov function by approximating V A .However, as is demonstrated in Section 5 it is not necessary to use all attractors in A to define a Lyapunov function for A that is constant on orbits in M A .More precisely, using the lattice structure of A we prove in Theorem 5.2 that a Lyapunov function can be defined using only the join irreducible attractors.
The next computational issue that needs to be addressed is that of identifying a finite bounded sublattice A ⊂ Att(X, f ).This involves approximating the dynamical system f for which we make use of the structure of M A .To explain this structure requires some standard concepts.
Given x ∈ X, consider a function γ x : Z → X, satisfying γ x (0) = x, and γ x (n + m) = ϕ(n, γ x (m)) for all n ∈ Z and all m ∈ Z + .Observe that given our assumptions on f , γ x need not exist.However, if it does its image, also denoted by γ x , is called a complete orbit through x.A set S ⊂ X is invariant under ϕ if and only if for each x ∈ S there exists a complete orbit γ x ⊂ S. The restriction to t ≤ 0 gives the backward orbit γ − x .For a point x ∈ X and for a backward orbit γ − x the orbital alpha-limit set is defined as where cl denotes closure.
Definition 1.1.Let S be a compact invariant set for f .A finite collection M = {M (p) ⊂ S | p ∈ P } of compact, pairwise disjoint, invariant subsets of S, labeled by the poset (P, ≤), is a Morse decomposition in S if for every x ∈ S \ p∈P M (p) and every orbit γ x ⊂ S there exist p, q ∈ P with q < p such that α o (γ − x ) ⊂ M (p) and ω(x, f ) ⊂ M (q).
Each set M (p) is called a Morse set.
As defined by (3), M A is a Morse decomposition [6], and V A is constant precisely on each of the Morse sets.As is made clear below, we identify A by identifying a Morse decomposition.
Observe that, up to this point, the discussion applies to any dynamical system defined on a compact metric space X and that we have reduced to problem to one in which we need to approximate a continuous function f : X → X.An efficient universal means of approximating continuous functions on arbitrary compact metric spaces is probably not possible, thus at this point we restrict the breadth of problems being considered.In particular, for the remainder of this paper we assume that X is a rectangular subset of R d .The justification of this choice is twofold.First, a wide variety of models of dynamics, whether in the form of maps or differential equations, are framed in R d .In addition, many of the specific techniques employed in this paper have been extended to the setting of compact triangulable manifolds via chart maps, which again employ computations on rectangular subsets of R d (e.g. the "Atlas" feature in the Conley-Morse-Database software [12]).Thus, from the perspective of applications working in R d seems to be a reasonable restriction.Second, our approximation method utilizes the framework for computational global dynamics developed in a series of papers [16,1,4,3,18] and is based on an adaptive nonuniform discretization of X.For rectangular sets in R d , the construction of the discretization can be done straightforwardly via bisections resulting in a nonuniform rectangular grid denoted by X .
Having established the goals and the generality of the results of this paper we now turn to an brief description of the contents.Section 2 describes the combinatorial approximation of the dynamics.This takes the form a directed graph F whose vertices are defined in terms of the rectangular grid X and whose edges provide bounds on the images of f .Our focus is on describing the data structures since this obviously impacts the size and dimension of problems that can be considered.In particular, we describe the "TreeGrid" data structure that is used to represent X and compute F and indicate how this is implemented using succinct binary trees.We also provide a brief discussion on the expected memory savings provided by this approach.
Section 3 describes the use of Tarjan's algorithm [24] to identify a linearly ordered collection of sets of grid elements M(p) ⊂ X , p ∈ P called combinatorial Morse sets with the fundamental property that the maximal invariant sets in M(p) determine a Morse decomposition for f .We also discuss the fact that for large applications because of the number of edges involved, storing F in memory is impractical.For this reason we employ a modified algorithm [11] for which the memory requirements are of the same order as the size of X while preserving the run time of Tarjan's algorithm.
Section 4 provides algorithms that take the combinatorial Morse sets as input and determines pairs of sets of grid elements A, A * ⊂ X such that A approximates an attractor and A * approximates the dual repeller.
Section 5 provides the theory and algorithms for the construction of a piecewise constant approximation of a Lyapunov function.We begin by recalling the relationship between the partial order structure of Morse decompositions and the lattice structure of attractors that is used to provide a proof of Theorem 5.2 cited above.We then turn to the construction of an approximate Lyapunov function based on a combinatorial attractor repeller pair (A, A * ).To be more precise about the form of this approximate Lyapunov function, we introduce the following definition.Definition 1.2.Given a combinatorial attractor repeller pair (A, A * ) a function V A : X → [0, 1] is a combinatorial Lyapunov function for f if it satisfies the following conditions: A key step in making the construction of combinatorial Lyapunov functions practical is the introduction of efficient algorithms for evaluating the distance between points in X and the sets A and A * .These algorithms culminate in Proposition 5.7 which provides explicit upper bounds on the memory and time costs for the computation of the combinatorial Lyapunov function based on the dimension of the problem and the approximation associated with the combinatorial Morse decomposition {M(p) ⊂ X | p ∈ P}.
It should be noted that combinatorial Lyapunov functions follow directly from the algorithms described in Section 3. Unfortunately, as one refines the grid X and the approximation of the dynamics F, in general, these combinatorial Lyapunov functions will not converge to Lyapunov functions for f .Thus, from our perspective they should not be viewed as approximate Lyapunov functions for f .However, as we make precise at the end of Section 5 the construction we perform leads to combinatorial Lyapunov functions that converge to a Lyapunov function of f as the approximation of X and F are refined.
Section 6 provides results on the implementation of these techniques to two population models.These models are chosen because they exhibit complex recurrent dynamics with multiple basins of attraction.

Combinatorial representation of dynamics.
As stated in the introduction, we are interested in computing approximate Lyapunov functions for the discrete time dynamical system obtained by iterating a map f : X → X with X a compact rectangular subset of R d .Note that f need not be injective nor surjective.The computational method utilizes the framework for computational global dynamics developed in a series of papers [16,1,4,3,18] as described in the introduction.First, the phase space is discretized using a finite, rectangular subdivision of X into a union of closed subboxes which have potential overlap only on their boundaries, the collection of which is called a grid and denoted by X .We use the notation | • | : 2 X → R d to denote the realization map for sets of grid elements given by |S| = ∪ ξ∈S ξ.The dynamics is then approximated by a combinatorial multivalued map F : X − → → X which maps grid elements to sets of grid elements.A combinatorial map F can be represented as a directed graph whose vertices are the grid elements in X and which has an edge from ξ to η if and only if η ∈ F(ξ).
The concept of an outer approximation as a method for discretizing the phase space and approximating a discrete dynamical system was originally introduced in [23].
The fundamental property of an outer approximation F is that an orbit of the underlying dynamical system f is contained in grid elements corresponding to a walk through the directed graph of F, in other words for every orbit {x n } with In order to store a grid X and an outer approximation F in an efficient manner for large computations, appropriate data structures are necessary, and we must also take into account the efficiency of the operations on X and F that are required to perform the algorithms.Since one of the main goals of this paper is to explore the overall capability of computing Lyapunov functions, we now describe some of these implementation issues in some detail.
Let Rect denote the collection of grid-aligned rectangular prisms (i.e. products of intervals) in R d .In order to construct a suitable F and X given a formula for f on X, we require a data structure to represent X which provides the following functionality: Note that both of the methods Cover : Rect ⇒ X and Geometry : X → Rect allow for the possibility of overestimating.This allows for a trivial implementation where Cover returns all of X and Geometry returns a bounding box for X.Such a trivial implementation will result in correct, yet trivial results.Tighter enclosures are better.Topological structure is captured by Observe that for each ξ ∈ X , Neighbors(ξ) ⊂ (Cover • Geometry)(ξ) The grids we will consider will satisfy the tightness condition which precludes trivial implementations.
Assuming we can implement a function F : Rect → Rect such that we can produce a combinatorial map F which is an outer approximation (Definition 2.1) via the composition In practice one can usually implement F via rewriting the formula for f using interval arithmetic.We do note, however, that sometimes this procedure results in very large bounding boxes.In fact we can be somewhat more flexible and introduce an enlarged collection of acceptable geometric shapes Geo in R d such that Rect ⊂ Geo.Then we can let Cover : Geo ⇒ X , Geometry : X → Rect, and F : Rect → Geo.This allows the images of F (R) to be tighter around the true image f (R) since they can be chosen from a more refined family of possibilities.For example, we could choose Geo to be sets which can be constructed from a finite union of Rect objects.Or we could let Geo consist of all parallelepipeds.However, in the interests of simplicity, for this paper we will simply take Geo = Rect.
2.1.Binary Tree Based Grids.We implement grids using a binary tree structure.For simplicity we take the phase space X to be a Rect object itself.We take the root of the binary tree to correspond to the Rect object X.We choose a splitting dimension, and subdivide the Rect object corresponding to the root node into two congruent Rect objects via a hyperplane bisection.We take these two Rect objects to correspond to the two children nodes.The left child will correspond to the Rect object with coordinates below the hyperplane, and vice-versa for the right child.We continue this procedure as we please to make a binary tree of any shape.See Figure 1.We impose the constraint that the splitting dimension is a function of the depth alone; in our implementation the subdivision is done round-robin through the dimensions so that the splitting dimension at depth k is k mod d.
The binary tree one constructs via this procedure corresponds to a grid, where the grid elements are taken to be the leaves of this tree.Notice that this binary tree is full, meaning that a node is either a leaf or has precisely two children.In a computer implementation one must refer to these leaves/grid elements in some fashion, so we use contiguous integers (0, 1, 2, • • • , N − 1), where N is the total number of grid elements.We label them according to the order the leaves are visited in a preorder traversal (the ordered depth-first-search of the tree starting at the root).We will call this leaf indexing.In the course of devising algorithms, we also need names for the nodes of the tree as well.For this purpose we will again use the order in which the nodes are visited in a preorder traversal, so the nodes are labelled (0, 1, • • • , M − 1).We will call this node indexing.Note that a leaf has both a leaf index and a node index, which need not be equal, and for a full binary tree we must have M = 2N − 1.
FIGURE 1.This figure illustrates a correspondence between Grid subdivisions and binary trees, as well as the node and leaf indexing scheme we use.
We will assume the binary tree data structure provides the following functionality: TreeGrid Data Structure Interface 1. parent, left, right.Given the node index for a node, return the node index for its parent, left child, or right child (respectively).If there is no such node (i.e. the root has no parent and leaves have no left or right child) then return M .2. GridToTree.Given a leaf index, return the node index of the corresponding leaf.3. TreeToGrid.Given a node index, give the leaf index of the corresponding node if the node is a leaf.If the node is not a leaf, return N .4. bounds.Return the Rect object corresponding to the root node.
Proposition 2.2.The TreeGrid interface can be used to implement a Grid satisfying the tightness condition (4).
Proof.We give implementations of Geometry and Cover.For Geometry, we can determine the Rect object for a grid element, which will be referred to by its leaf index, by calling GridToTree, calling parent repeatedly to find the path to root, and along the way checking if we came from the left or right child by calling left and right and comparing the node indices.This results in path-to-root information, which tells us a sequence of subdivisions which when applied to bounds will result in the correct Rect object to return.
Cover can also be implemented efficiently using the binary-tree structure.When presented with a Rect object to cover, we begin by seeing if it intersects with the Rect corresponding to the root of the tree.If it does intersect, we recurse and ask this question for both of its children.If it does not intersect, we do nothing and continue on with other branches.When we find a leaf node that intersects, we call TreeToGrid and place the appropriate leaf index in the list of results for Cover.
Provided the numerics used in handling computer representations of Rect objects, which involve floating point representations, are performed with care, it is straightforward to see the tightness condition (4) can be achieved.

Succinct Binary Tree Grid Implementation.
One way of implementing a data structure respecting the TreeGrid interface is to create a binary tree structure by creating "Node" objects in memory that contain pointers to parent, left child, right child, and an integer containing the node index.We could use a size M array to store pointers to these Nodes.GridToTree and TreeToGrid could be implemented by storing arrays of integers of length N and M , respectively.The space usage of such an implementation would be 6M + N ≈ 13N words of storage.A word of storage is the amount of computer memory used to store a pointer.In this era, it is usually 64 bits.Hence this results in a data structure that requires roughly 104 bytes per grid element.
The pointer-based method is quite reasonable in many situations.Its drawback is its space usage.At the expense of an up-front constant in time usage, it is possible to substantially bring down the space requirements by using methods from the succinct data structures literature, in particular succinct balanced parentheses lists and rank-select structures [20,5,13].Specifically, it is possible to construct a succinct TreeGrid data structure providing the required methods in constant time and using only 2N + o(N ) bits of space.Ignoring the o(N ) overhead, this means that for 64-bit words the succinct methods are roughly (13N • 64)/2N = 416 times as space-efficient as a pointer-based method.
We describe how we accomplish this.First, we review the succinct data structures at our disposal.The key features of a succinct balanced parentheses list are that they require only 2n + o(n) bits of space to store n pairs of matching parentheses, and they allow the operations of findopen and findclose, which find opening and closing matching parenthesis, respectively, in constant O(1) time.Meanwhile, a rank-select structure is a succinct data structure that requires n + o(n) bits of space and is capable of answering rank and select queries about an n bit sequence in O(1) time.A rank(i) query asks how many 1 bits occur strictly before some position 0 ≤ i ≤ n in the sequence.A select(i) query asks for the maximal j such that rank(j) = i.We also assume the ability to access the underlying bit sequence in constant time.
We utilize these structures in the following way.Given a full binary tree with N leaves, and hence M := 2N − 1 nodes, we construct an M bit sequence we call the leaf sequence such that [k] = 1 if and only if the node with node index k is a leaf.A leaf sequence completely and succinctly characterizes a full binary tree.Note that the last bit in the leaf sequence is a 1 so that the preorder traversal ends at a leaf.Hence the first M −1 bits of the leaf sequence has an equal number of 0's and 1's.In fact, interpreting a 0 (non-leaf) as an opening parenthesis and a 1 (leaf) as a closing parenthesis, the first M − 1 bits are a balanced parentheses list.Accordingly, we can build a rank-select query structure and a balanced parentheses query structure on top of this for an additional o(M ) space.Proposition 2.3.Let T be a full binary tree with M nodes.The queries findopen, findclose on the M -bit leaf sequence can be used to implement parent, left, and right via the following formulas, where we refer to nodes via their node index: x is not a leaf m (undefined) otherwise ( 7) Moreover, the queries rank and select on the leaf sequence implement TreeToGrid and GridToTree, respectively.
Proof.The key idea is that each pair of matching parentheses in the leaf sequence corresponds to a pairing between a tree node x (the open parenthesis) and the last preordered leaf y on x's left subtree (the closing parenthesis).This is very useful since the next node preordered after the left subtree of x, namely y + 1, is the right child of x.This yields the formulas.The final statement follows directly from the definitions.
Proposition 2.3 implies we can construct the succinct "TreeGrid" data structure requiring only the 2N + o(N ) bits of space advertised above.We should emphasize that this is by no means the only possible construction to achieve these ends.For example, the pioneering paper of Jacobson [13] describes a method using rankselect on level-order bitmaps that is suitable to binary trees such as ours that could achieve 4N + o(N ) bits.Our construction (which we presume is not original) is interesting because it is ultra-succinct [14], that is, it beats the 4N +o(N ) informationtheoretic bounds for representing general binary trees by utilizing special structure (namely, fullness).
In the computation of combinatorial Morse sets, we also have occasion to deal with subgrids, and the strategy our implementation uses is to wrap the above data structure with a bit-sequence which "knocks-out" extra leaves from a full binary tree to achieve an arbitrary binary graph.For an arbitrary binary tree with M nodes, this results in at worst 2M +o(M ) bits for a full tree containing the arbitrary binary tree, and another M + o(M ) bits to represent a rank/select structure on the knock-out sequence.This does not qualify as succinct; we would need to achieve 2M + o(M ) bits worst-case, not 3M + o(M ) bits in the worst-case.However, in practice the subgrids that arise in computational dynamics do not approach this worst case, because nodes with only one child tend to be rare.Accordingly, we consider the best case, where the tree is "almost full".In this situation one sees our strategy achieves roughly 1.5M + o(M ) bits.Hence despite not technically being succinct in the worst case, in the best case, which we claim will be the average case in practice, we obtain ultra-succinct performance.We discuss the performance of our implementation, which is based on the Succinct Data Structures Library [10], in Section 6.

3.
Extracting recurrent and gradient-like dynamics.Since we do not assume that the continuous map f : X → X generating the dynamics is injective or surjective, X need not be invariant.We compute Morse decompositions with respect to the maximal invariant set in X, As we described in the previous section, orbits of a dynamical system are rigorously contained in walks through the directed graph generated by an outer approximation of that system.In particular this implies that if an orbit {x n } is chain recurrent, then it is contained in a periodic cycle in the directed graph, i.e. there exists a periodic sequence Therefore, the problem of computing a neighborhood of the chain recurrent set can be reduced to finding the set of cyclic vertices in a directed graph, which we will call the recurrent set of F.
Analogous to the chain recurrent set of a dynamical system, the recurrent set of a graph is naturally partitioned by the reachability relation.That is, the equivalence relation ξ ∼ η, if there exists a walk from ξ to η and a walk from η to ξ in F, partitions the recurrent set of F into its recurrent components, which are also known as the strongly connected path components of the digraph.These recurrent components have a partial order induced by the reachability relation in the digraph of F. Corollary 4.2 in [16] implies that once we have computed the recurrent set of F, the geometric realization of the recurrent components are isolating neighborhoods of a Morse decomposition with the same partial order, i.e. we can rigorously extract a Morse decomposition from an outer approximation F : X − → → X .Due to this connection with a Morse decomposition, we denote the recurrent components by M. We emphasize that the recurrent components do not partition the digraph of F, since elements in the recurrent set must be cyclic through a nontrivial walk in F. In contrast, the strongly connected components do partition the digraph.Each recurrent component is a strongly connected component.However, each noncyclic vertex is also itself a strongly connected component but is not a recurrent vertex.
The classical algorithm for extracting strongly connected components of a directed graph given by its adjacency lists (lists of target vertices given a source vertex), is called Tarjan's algorithm [24].For a graph with V vertices and E edges, Tarjan's algorithm requires O(E+V ) time and O(E+V ) space to identify the collection of strongly connected components S. One identifies the recurrent components M as the strongly connected components containing at least one edge; i.e. components with more than one vertex as well as the single vertex components with self loops.
Tarjan's algorithm outputs the strongly connected components in a useful order.Recall that any directed acyclic graph admits a topological sort, which is a total ordering of the vertices such that a < b only if b cannot reach a.This concept applies to strongly connected components of a directed graph since collapsing these components to single vertices renders the directed graph into an acyclic one, the socalled condensation graph.The order in which Tarjan's algorithm produces strongly connected components will correspond to a reverse topological sort of the condensation graph.We will make use of this property later, when we describe algorithms that require iterating through strongly connected components S or recurrent components M in the order of a topological sort or reverse topological sort.
We point out that the O(E) space bounds arise from storage of the adjacency lists and could be quite large.Indeed, E could be as large as V 2 in worst case.Thus the space requirements for Tarjan's algorithm could dominate over the space requirements of the "Grid" data structure.Observe that size of the adjacency lists computed via Equation ( 5) are easily bounded from below by the product of the eigenvalues of Df that are larger than one.In general it will be considerably larger.In other words, storage of the graph can be quite expensive.
There is a solution to this problem.A modified version of Tarjan's algorithm, due to one of the authors and motivated by precisely this application, circumvents the need to store the adjacency lists and provides an O(V ) space bound [11].This algorithm requires that adjacency lists can be retrieved or computed on demand.Moreover, during execution it only asks for each adjacency list precisely once.Consequently it runs in time comparable to Tarjan's algorithm.These properties make it well suited for computing strong components of the directed graph arising from the combinatorial map F of Equation ( 5).This algorithm has been implemented and tested, and is publicly available in Conley-Morse-Database [12].
We have also recently become aware of the strong components algorithm of [15], which assumes adjacency lists are available in external memory and attempts to minimize I/Os.While it was not the author's intention, in fact it can be reinterpreted as an algorithm for computing strongly connected components which requires only O(V ) space provided one is given the ability to compute adjacency lists on demand.Detailed comparisons of up-front constants between these alternatives will be made in [11].
4. Attractor-repeller pairs.In the previous section we observed that graph theoretic techniques provide a naturally efficient way to obtain an approximate Morse decomposition from an outer approximation.To recover an approximate Lyapunov function from this information, we need to describe the relationship between Morse decompositions and attractor-repeller pairs.We assume that the reader is familiar with basic concepts of limit sets, attractors, and repellers in the context of a dynamical system.However, since we are working with semidynamical systems that are not necessarily invertible, we must address some subtleties that do arise; also we must describe these concepts in the context of combinatorial multivalued maps.
Since X need not be invariant, A * also need not be invariant.The phase space X is only forward invariant in general, and A * is the maximal forward invariant set in the complement of an attracting neighborhood U of A. For combinatorial multivalued maps F define the ω-limit set and α-limit set of a set U ⊂ X by The dual repeller A * to an attractor A is defined by A * = α(X \ A).Given an attractor A, the pair (A, A * ) is called an attractor-repeller pair for F. The set of all attractors is denoted by Att(X, F).
An attractor-repeller pair is a Morse decomposition with two Morse sets, the attractor and its dual repeller, with the partial order where the repeller is greater than the attractor.Just as in the case of Morse decompositions, [16,Proposition 5.5] implies that if F is an outer approximation, then the realizations (|A|, |A * |) of an attractor-repeller pair (A, A * ) are attracting and repelling neighborhoods respectively for some attractor-repeller pair (A, A * ) for f .
Recall that in Section 3 we describe an efficient algorithm to obtain the following information from the digraph representing F: 1. the recurrent components M = {M (p) | p ∈ P}, 2. the strongly connected components S, and

Global variables:
A topologically sorted list S of strongly connected components, and an array SCC indicating which strongly connected component each grid element is in.Proof.The correctness of ComputeAttractor requires showing that it computes ω(M ) = Γ + (M ) for a recurrent set M ∈ M.This is the case provided For-wardReachable is correct, since it promises to provide all strongly connected components reachable from M .Correctness of ComputeDualRepeller requires that we show it computes α(X \ A).It is straightforward to verify that α(X \ A) is the subset of grid elements which can reach some recurrent set M in X outside of A. These grid elements comprise the strongly connected components which can reach recurrent sets M which are not in A. Hence the algorithm for ComputeDualRepeller is correct assuming the correctness of BackwardReachable.

ForwardReachable:
Next, we show ForwardReachable is correct and executes in O(E) time, which follows from that disjointness of strongly connected components.Indeed we have that Adjacencies(u), which takes O(card(F(u)) time, is called at most once per vertex u ∈ X .Moreover, "Insert SCC[v] into A", which takes O(1) time, is executed at most card(F −1 (u)) times per vertex u ∈ X .This counts each edge of the directed graph at most twice, so that the sum of this time can thus be bounded by a constant factor times the number of edges of the directed graph corresponding to F, which gives a total complexity of O(E) time.Correctness of ForwardReachable follows from induction given (1) that S is processed in topologically sorted order and (2) the observation that if a strongly connected component S belongs to the forward reachable set of A, then either S ∈ U or there exists a strongly connected component T which directly reaches S, that is, a vertex of S is in the adjacency list of a vertex in T , and hence T necessarily occurs before S in the topological sort of S. A similar argument shows that BackwardReachable is correct and computes in worst case O(E) time.
Since card(P) < E, it follows that the complexities of ComputeAttractor and ComputeDualRepeller are dominated by the complexities of their calls to For-wardReachable and BackwardReachable, and hence are both O(E) as well.Since AttractorRepellerPairs calls ComputeAttractor and ComputeRepeller card(P) times, its complexity is bounded by O(card(P) • E), as claimed.
5. Lyapunov functions.Algorithm (1) takes Morse sets in a Morse decomposition and produces a collection of attractors.However, there exists a deeper structural relationship between Morse sets and attractors that we now describe.We use this relationship to restrict the set of attractors used to define a Lyapunov function for the Morse decomposition.As discussed in the Introduction, this Lyapunov function is defined in terms of Lyapunov functions for attractor repeller pairs.Thus we address the algorithmic issues of constructing a combinatorial Lyapunov function in this limited setting, before extending it to the full Morse decomposition.

A Lyapunov function for Morse decompositions.
We continue to study the dynamics generated by a continuous map f : X → X.Throughout this section we set S := Inv(X, f ).For the remainder of this subsection we restrict our attention to f : S → S. Since we are focussing, for the most part, on the relationship between attractors and Morse decompositions and since Att(X, f ) = Att(S, f ), this is not a serious restriction.
As indicated in the Introduction Att(S, f ) is a bounded distributive lattice.To be more precise given A, A ∈ Att(S, f ) the lattice operations are defined by and the minimal and maximal elements are 0 = ∅ and 1 = Inv(X, f ).
Let M = {M (p)} be a Morse decomposition of S labelled by the poset (P, ≤).Given a finite poset (P, ≤), the set of all down sets of P is defined by A direct calculation shows that O(P) is a finite, distributive lattice under the operations of union and intersection.In [19] it is shown that the map µ : O(P) → Att(S, f ) defined by is a lattice embedding, where W u (M ) is the unstable set of M .The image A = µ(O(P)) is then a finite bounded sublattice of Att(S, f ) which is isomorphic to O(P).
In this way the attractor A is determined uniquely from the Morse decomposition M and the partial order (P, ≤).
We are interested in Lyapunov functions that are compatible with Morse decompositions which leads to the following definition.
is called a Lyapunov function for (X, M, P, ≤).
From the definition of Morse decomposition, one can check that for every attractor A ∈ A each Morse set is contained in either A or A * , and in addition Following the original ideas of Conley, for any set of weights β A > 0 with A∈A β A = 1 and any choice of Lyapunov functions V A for attractor-repeller pairs (A, A * ) the function is a Lyapunov function for the Morse decomposition M with order (P, ≤), see Definition 5.1 and also [16,Equation 7]).However, using all attractors in A to define such a Lyapunov function is more than necessary, and to develop a computationally efficient algorithm it is desirable to find a smaller collection of attractors that will generate a Lyapunov function.The optimal such set is ultimately determined by the specific structure of the poset (P, ≤), but there is a natural subset that works which can be defined for all partial orders and which typically yields a significant reduction in the number of attractors that need to be considered.The elements in O(P) can be written as unions of down sets of the form ↓ p = {q ∈ P | q ≤ p}, which are called the join-irreducible elements of the lattice O(P).The set of join irreducibles is denoted by J ∨ (O(P)).Therefore each attractor in A can be written as a union of attractors of the form A p = ν(↓ p), i.e.J ∨ (A) = {A p | p ∈ P}.Theorem 5.2.Let S be a compact, invariant set, let M be a Morse decomposition for S, labeled by a poset (P, ≤), and let A be the corresponding sublattice of attractors whose join irreducible elements are A p .For each p ∈ P let β p > 0 with p∈P β p = 1 and V p be a Lyapunov function for the attractor-repeller pair (A p , A * p ). Then is a Lyapunov function for (S, M, P, ≤).
Proof.From the definition of A r = ν(↓ r) we have that Hence V (M p ) = p ≤r β r .If q < p, then c q = V (M q ) = q ≤r β r < p ≤r β r = V (M p ) = c p .This establishes conditions (i) and (ii) in Definition 5.1.Now suppose x ∈ S \ ∪ p∈P M (p).Since M is a Morse decomposition, there exist p, q with q ≤ p such that ω(x) ⊂ M (q) and α o (γ − x ) ⊂ M (p).Then M (q) ⊂ A q , M (p) ⊂ A * q , and x ∈ S \ (A q ∪ A * q ), which implies that V q (f (x)) < V q (x).Since each function V p is decreasing, V (f (x)) < V (x) This establishes condition (iii) in Definition 5.1 and completes the proof.
Note that in Theorem 5.2 the partial order (P, ≤) on the Morse sets is given.Any coarser order obtained by adding more relations is also admissible, but such a coarser order relation results in a loss of dynamical information about the system.In particular, the set of join irreducible elements may grow, thus requiring additional attractors to be included in the definition of V , i.e. in (10).

Approximating Lyapunov functions for attractor-repeller pairs.
Let A ∈ Att(X, f ).We begin by describing the the standard analytical construction of a Lyapunov function for an attractor-repeller pair (A, A * ).The distance potential dp A : X → [0, 1] defined by Moreover, dp A is Lipschitz since A and A * are disjoint and compact.We incorporate the dynamics of f through the family of functions which maximize the distance potential over forward orbits.The function is a Lyapunov function for (A, A * ), as is shown in [16].Following [16,3] this function can be approximated by a function on grid elements as follows.Let A ∈ Att(X , F).With respect to the attractor repeller pair where x ξ is the centroid of ξ.Then the function ν * A : X → [0, 1] defined by maximizes the value of dp A over the complete forward image of ξ.Finally we consider the function V A : X → [0, 1] defined recursively by Lemma 3.2 and Theorem 3.3 in [3] imply that if (A, A * ) is an attractor-repeller pair for F that well approximates the attractor-repeller pair (A, A * ), then V A is a combinatorial Lyapunov function for (A, A * ) which well approximates V A .We make this statement more specific in Section 5. 4.
In practice, computing the distance potential dp A using distances between sets of grid elements is not computationally efficient.However, we now describe an O(N log N ) time algorithm, where N is the number of grid elements in X , to compute a different combinatorial distance potential cdp A : X → R given a combinatorial attractor-repeller pair (A, A * ).We show it also approximates the true distance potential (11) arbitrarily closely as we refine the grid.Our algorithm works in the case where the distance metric d is chosen to be Manhattan distance and the grid implements the "TreeGrid" interface discussed above so that there is an underlying binary tree.We remark that this algorithm is a vast improvement over a naive approach where one calculates via brute force the minimum distance from each grid element of X to each grid element in the attractor and in the repeller, which would require O(N 2 ) time.

Dijkstra's Algorithm.
A key subroutine of our distance potential algorithm will be Dijkstra's algorithm [8] for computing shortest paths on a graph W with weighted edges.Recall that Dijkstra's algorithm proceeds by seeding an initial vertex set S with a priority value 0, and placing these vertices in a priority queue, with the highest priority being given to those vertices with the smallest priority value.A vertex u with highest priority is processed by marking it as processed, and placing each unprocessed neighbor v of u into the priority queue with a priority value equal to the priority value of u plus the weight of the edge from u to v. When a vertex is inserted into the priority queue multiple times, the effect is the vertex will have the smallest priority with which it has been inserted.This continues until all vertices have been processed.The Dijkstra distance of a vertex v to the set S, which we will denote d W (v, S), is the priority value attached to the vertex v when it is processed.If the algorithm is implemented with a so-called Fibonacci heap, then this algorithm, Dijkstra(S, W), requires O(E + V log V ) time and O(V ) space [9], and returns the Dijkstra distances d W (v, S) for all v ∈ W.Here the log V cost is associated with the sorting necessary in the priority queue.
Our algorithm for computing a combinatorial distance potential requires a weighted graph W to which Dijkstra's algorithm is applied.To this end, we utilize the binary tree structure underlying the grid X .In particular, we realize the vertices of W as the nodes in the binary tree underlying X .We will abuse notation and regard X as a subset of the vertex set.What is left is to specify the edges of W and their weights.To this end, first recall that to each such node v ∈ W there is an associated rectangular box in R d , and we will use the notation |v| to denote this geometric realization.
We let (u, v) be an edge in W if one of two conditions holds: • Condition (A): u and v have a parent-child or child-parent relationship.
• Condition (B): The rectangles |u| and |v| are congruent and share a codimension-1 face.
In either case, we assign a weight to the edge (u, v) which is the Manhattan distance between the midpoints of |u| and |v|.
The motivation for this weighted graph comes from the following observation.Call a path in R d that is comprised of a finite sequence of axis-aligned treks a Manhattan curve.Paths in the graph W correspond to Manhattan curves, and the weights on the edges correspond to the Manhattan distances of the treks of which they are comprised.Hence the Dijkstra distance on this graph gives an upper bound for the Manhattan distance between the grid elements.Shortly, we will show that the amount of overestimate can be bounded.
We omit the algorithm for producing the weighted adjacency list of a vertex v in the weighted graph W and content ourselves with only a few comments.It is straightforward to give the adjacencies due to Condition (A).Condition (B) can be re-expressed it in terms of the path-from-root information, i.e. the sequence of left/right directions which one must follow from the root of the tree to find the vertex.Expressing it this way provides a method to construct an efficient algorithm.Given path-from-root information, we split it into d lists such that the k-th list has the path information at depths {i | i mod d = k}.The idea is that the k-th list has the path information corresponding to choices made for splitting dimension k.If for each such list we construct a sequence of digits by a digit 0 for a left move and a 1 for a right move, we obtain d binary sequences.These binary sequences we may then reinterpret as nonnegative integers written in binary notation.Call these the coordinates of the tree node, and indeed for a uniformly subdivided grid they would be precisely the natural integer coordinates one would put at a fixed depth.Supplemented with depth, the coordinates provide an alternative characterization of the tree node.Using this terminology, Condition (B) states that u and v are the same depth and that the coordinates of u and v are identical in all but one dimension, and in the dimension the coordinates differ, they differ by precisely one.With this insight it becomes straightforward, albeit somewhat tedious, to produce an efficient algorithm that relies on tree-moves for computing the weighted adjacency list of a vertex v.We find the coordinates for v, find all of the at most 2d coordinates which are neighboring in the sense of Condition (B), and locate them via tree moves left and right starting at the root by following the directions of the binary digits of the new coordinates in round-robin fashion.
Proposition 5.3.The Dijkstra distance on W approximates the Manhattan distance in the following sense.For each pair of grid elements ξ 1 , ξ 2 ∈ X , and for each where δ is the maximal diameter of the grid elements, i.e. δ := max{diam (ξ) | ξ ∈ X }.
Proof.Let x ξ1 , x ξ2 be the centroids of ξ 1 , ξ 2 , respectively.Observe ), since the Dijkstra distance in W between ξ 1 and ξ 2 measures the Manhattan distance of some Manhattan curve between x ξ1 and x ξ2 .We need only bound how much d W (ξ 1 , ξ 2 ) may be overestimating d(x ξ1 , x ξ2 ).In order to do this, we give a path in W which we can show does not overestimate by too much.We construct this path as follows.Let D be the largest number so that every leaf in the binary tree for X is at least depth D. Let η 1 and η 2 be the unique ancestors at depth D of ξ 1 and ξ 2 in the binary tree for X .Note that we might have ξ 1 = η 1 or ξ 2 = η 2 in this construction.Consider the path in W that starts at ξ 1 , climbs through parents until depth D is reached at η 1 ∈ W, walks to adjacent neighbors at depth D, which are defined by Condition (B), in the most straightforward way to reach η 2 ∈ W, and then descends through the appropriate path of children in order to reach ξ 2 .Let x η1 , x η2 be the centroids of η 1 , η 2 , respectively.Observe that d W (η 1 , η 2 ) = d(x η1 , x η2 ) since the adjacencies given by Condition (B) of W at depth D provide a uniform grid.From this we may realize the bound The latter inequality holds since the points x ξi , x ηi for i = 1, 2 are both contained in a rectangle of diameter less than δ.For each k = 0, 1, 2, • • • , the sum of the weights of parent-child edges leading from depth D + kd to depth D + (k + 1)d is at most δ/2 1+k since it consists of d treks whose lengths most be less than the diameter of a depth D + kd tree node.Hence a geometric series argument gives us d W (ξ 1 , η 1 ) < ∞ k=1 δ/2 1+k = δ and similarly d W (ξ 2 , η 2 ) < δ.This gives us d W (ξ 1 , ξ 2 ) − d(x ξ1 , x ξ2 ) < 4δ, which can be combined with |d(x 1 , x 2 ) − d(x ξ1 , x ξ2 )| ≤ δ to give the desired result.Proof.Proposition 5.3 gives d W (ξ, η) ≤ 5diam (X ) + d(x ξ , y) for all η ∈ S and all y ∈ |S| which implies min η∈S d W (ξ, η) ≤ 5diam (X ) + min y∈|S| d(x ξ , y).Hence follows by the same reasoning.

Combinatorial Distance Potential.
Given A ∈ Att(X , F) and the corresponding attractor-repeller pair (A, A * ), we can now define our formula for the combinatorial distance potential cdp A : X → R, for all ξ ∈ X , We give the following algorithm to compute cdp A .
Output: A distance potential cdp A and also the Dijkstra distance d W (A, A * ) Realize the weighted graph W corresponding to the grid X .
To this end we provide the following algorithm.

Algorithm 3 Computing Distance Potential Star
Global: The topologically sorted list S of the strongly connected components

DistancePotentialStar
Given: An attractor-repeller pair (A, A * ) Output: The combinatorial distance potential star as defined in Equation ( 18) and also the Dijkstra distance d W (A, A * ).
} end for end if end for Proposition 5.6.The algorithm DistancePotentialStar correctly computes the combinatorial distance potential star defined in Equation (18), with time complexity O(dN + N log N + E) and space complexity O(N ).Here N is the number of grid elements and E is the number of edges in the directed graph corresponding to the combinatorial map F of Equation (5).
Proof.The correctness of the algorithm follows from the observation that cν * P satisfies the following recursive relation: Iterating through S \ (A ∪ A * ), which consists only of singleton components, in a reverse topological sort ensures that cν * A (η) is always computed before cν * A (ξ) whenever η ∈ F(ξ), so this recursive definition is properly implemented.Proposition 5.5 gives us a baseline resource usage of O(dN + N log N ) time and O(N ) space.The time cost apart from this is proportional to the number of adjacencies examined, and no adjacency ξ → η is ever examined more than once; this results in the extra O(E) term in the time complexity.The extra space complexity is due to storage of the result and can be subsumed into the big-oh notation (in fact, the distance potential can be overwritten).5.5.

5.2.4.
Computing Lyapunov functions for attractor-repeller pairs.Finally we replace cν * A in formula (15) for V A to obtain the function which approximates V A .Lemma 3.2 in [3] is the number of edges in the directed graph corresponding to the combinatorial map F of Equation (5).
Proof.The proof is essentially the same reasoning as the proof of Proposition 5.6; we see that iterating through singleton strongly connected components in reverse topological sort correctly implements the recursive definition of Equation (19).

Computing Lyapunov functions for Morse decompositions.
Our goal is to compute combinatorial Lyapunov functions for F, i.e. piecewise-constant, combinatorial Lyapunov functions that approximate continuous Lyapunov functions for f .Given an outer approximation F, the recurrent set R has finitely many recurrent components M. Corollary 4.2 in [16] implies that f is gradient-like on the complement of |R|, and the realizations of the recurrent components are isolating neighborhoods for a Morse decomposition for f .This means that the maximum amount of recurrence/gradient-like information we can obtain from F is a Morse decomposition M with partial order (P, ≤).
Theorem 5.2 provides a formula for a Lyapunov function for a Morse decomposition as a weighted average of Lyapunov functions for attractor-repeller pairs generated from the join-irreducible attractors.This leads us to define the following piecewise-constant, combinatorial Lyapunov function for where the weights β p are defined according to the distance from the attractor A p to its dual repeller A * p .Specifically, In the next section we discuss the reasoning behind this specific choice of weights and describe how well these functions approximate continuous Lyapunov functions for f .The following algorithm computes the final combinatorial Lyapunov function cV.
n → ∞, the corresponding functions cV n converge uniformly to a continuous Lyapunov function for f .First, it should be clear that without any further assumptions about how well F n approximates the images of f , nothing can be said.For the purposes of clarity, we consider the minimal outer approximation of f on X n given by Clearly, this outer approximation is not attainable in practical computations, but the notion of a convergent sequence of outer approximations is developed in [18].We do not include the details of such convergent sequences here, however the results that we obtain for sequences of minimal outer approximations will also hold for convergent sequences of outer approximations using Proposition 5.4 and Corollary 5.6 in [18].
We begin by considering Lyapunov functions for attractor-repeller pairs.Proposition 5.5 in [16] establishes a one-to-one correspondence between well-separated attractor-repeller pairs for f and attractor-repeller pairs for the minimal outer approximation F o , and Corollary 5.6 in [18] establishes the same result for convergent sequences of outer approximations.Specifically, for c > 0 let is sufficiently small, then for each A ∈ Q c there exists a unique A ∈ Att(X , F) with the property that We denote this correspondence by Π : Q c → Att(X , F).This leads to the following theorem.
Theorem 5.9.If A n is a sequence of attractors arising from a convergent sequence of outer approximations F n : X n → X n with diam (X n ) → 0 as n → ∞, and there exists an attractor A for f such that (|A n |, |A * n |) → (A, A * ) in the Hausdorff metric, then cV An converges uniformly to V A defined in equation (12).Moreover, given A such a sequence A n always exists.
Proof.Corollary 5.4 states that the combinatorial distance d W (ξ, S) and the geometric, Manhattan distance d(x ξ , |S|) satisfy Lemma 3.7 of [3] establishes that the distance potentials dp An and cdp An defined in equations ( 13) and (17) respectively, as well as the functions V An and cV An defined in equations ( 15) and (19) respectively, satisfy similar estimates so that there exists The second issue is due to spurious recurrent components, and could be handled by a strategy of grouping components as well.In this case we would group components together so that the distances between combinatorial attractors and their dual repellers of the join irreducible pairs were not too small.Pursuing these two strategies both computationally and theoretically will be the subject of future work.
6. Implementation and Results.The above algorithms have been implemented and are freely available in the C++ software package Conley-Morse-Database [4] [12].In this section we discuss the performance of the Succinct Grid and Pointer Grid described in Section 2 and computational results for Lyapunov functions for two different dynamical models.We also report a breakdown of the memory usage among subcomponents of our algorithms and discuss scalability prospects.6.1.Performance of Succinct Grid vs Pointer Grid.The Conley-Morse-Database software provides C++ classes PointerGrid and SuccinctGrid implementing both pointer-based and succinct-tree based grids described in Section 2. SuccinctGrid makes heavy use of the software package SDSL (Succinct Data Structures Library) [10] which provides a C++ library of succinct data structures including succinct trees and rank-select structures.PointerGrid, on the other hand, is a simple pointer-based scheme implemented according to the description given earlier.We measure that the succinct-tree based grid implementation achieves an efficiency of 4.58 bits per grid element, which represents a 180-fold increase in space efficiency compared to the pointer grid method we have implemented, which requires 104 bytes per grid element.This is more than the asymptotically achievable 3 bits per grid element (or 1.5 bits per binary tree node -see the discussion at the end of Section 2.2) due to the o(N ) usage of the succinct data structures taking approximately 34% of the space in the example we measured.There is the expected time/space trade-off: SuccinctGrid is the slower alternative by a factor of 3. We refer the reader to Table 2 below.

Computational Examples.
We have tested our algorithms on two examples, the first being two-dimensional and the second being three-dimensional.Both examples computed in a few hours using a single core on a Macbook Pro laptop with 8GB of RAM.For each example we present a Conley-Morse graph, a visualization of the combinatorial Morse sets, and a visualization of the computed combinatorial Lyapunov function.6.2.1.Two-dimensional overcompensatory Leslie Model.The first example is of the map f : R 2 → R 2 given by where we choose parameters θ 1 = 20.0,θ 2 = 20.0,φ = 0.1, and p = 0.7.The phase space region was taken to be X = [0, 74] × [0, 52].This model was first analyzed in Ugarcovici and Weiss [25].Later in [1] the multiparameter system was examined computationally via Conley-Morse theory.Here we compute the Conley-Morse graph using the program SingleCMG and a combinatorial Lyapunov function using Lyapunov.SingleCMG and Lyapunov are both programs in the Conley-Morse-Database project.
Resource requirements are described in Table 2 below.We briefly discuss the results.The Conley-Morse graph in Figure 2 has five nodes, corresponding to five combinatorial Morse sets that were found.The labels on these nodes tell us the Conley index information.We will not discuss how these are arrived at here, cf.[1].We content ourselves with describing intuitively what they indicate.To this end we provide the following table.Note that if an invariant set has a certain Conley index, this does not imply that the invariant set is the set indicated in the table.For example, in the model (23), the invariant set with the Conley index of a stable period-3 orbit is actually a 3-part chaotic attractor.Conley index of an unstable period-T orbit In the phase space picture we can easily see the Morse sets with the Conley indices of an invariant circle and unstable fixed point.The Morse sets with Conley indices of an unstable period-3 orbit and a stable period-3 orbit are small at this resolution so they are highlighted in blue and purple respectively.The Morse set with the saddle-point Conley index corresponds to the origin in the lower left corner, highlighted in red.
We remark that to compute Conley indices of combinatorial Morse sets on the boundary of phase space requires that we extend the phase space slightly into the negative numbers.If we do not, the boundary will cause trivial Conley indices to appear, which we consider an artifact.
In Figure 3 we see what we should expect to see: the combinatorial Lyapunov function is constant on Morse sets, and takes higher values on the Morse sets which are higher in the Conley-Morse graph.In this figure white is the highest value and black is the lowest value.6.2.2.Three-dimensional Leslie/Gower Competition Model.Our second example is a three-dimensional model considered by Cushing et.al. [7].The map is given by with parameters A = 5.0, B = 0.1, C = 0.11, D = 0.8, E = 5.0, F = 0.12, G = 0.1.The phase space region was taken to be X = [0, 20] 3 .We present our results for this model in Figure 5 and Figure 4. Supplementary materials are available for visualizing the 3D combinatorial Lyapunov function [22].6.3.Scalability.Our space requirements have been reduced with the succinct grid data structure in Section 2.2 and a space efficient technique for computing the strongly connected components, see [11].We present these resource usage statistics for the computation on the overcompensatory Leslie model of Equation ( 23) in Table 2.
From this table we observe that had we stored adjacency lists and used pointerbased grids we would have required 3748 + 1398 = 5146 MB.Dispensing with adjacency list requirements, by itself, allows us to bring this down to 3748 MB.In combination with the succinct grid data structure our space requirements are further brought down to 1704 MB.Thus, we have demonstrated a factor 3 of scaling in this paper.
However, we note that we have reduced memory requirements of objects that require random access patterns, and thus by necessity must be in internal memory.Other memory requirements remain, but it appears that many, if not all of them, can be relegated to external memory.See [26] for a survey of external memory algorithms and data structures.For example, the Dijkstra algorithm's priority queue can be implemented in external memory [2].This suggests we may have opened the door to orders-of-magnitude scaling via such methods.We leave such an investigation for future work. References.

14 AAlgorithm 1
. GOULLET, S. HARKER, W.D. KALIES, D. KASTI, AND K. MISCHAIKOW 3. a topological sort of S. Given a recurrent component M ∈ M, define Γ + (M ) to consist of all grid elements in X reachable from M .Since F(Γ + (M )) = Γ + (M ), Γ + (M ) ∈ Att(F).We show in Section 5 that the information in the attractors {A p := Γ + (M (p)) | p ∈ P} are sufficient to produce a combinatorial Lyapunov function.Consequently, we now give an algorithm to compute the combinatorial attractor-repeller pairs of this form using only the information from the connected components and topological sort.Computing Attractor/Dual Repeller Pairs Global variables: A list M of the combinatorial Morse sets, a topologically sorted list S of strongly connected components, and an array SCC indicating which strongly connected component each grid element is in.AttractorRepellerPairs for M ∈ M do A ← ComputeAttractor(M ).A * ← ComputeDualRepeller(A). end for ComputeAttractor Given: M ∈ M Return: Minimal attractor containing M Create an empty set U Insert M into U.A ← ForwardReachable(U) Return A. ComputeDualRepeller Given: Combinatorial attractor A ⊂ S Return: Combinatorial dual repeller of A Create an empty set U for

Proposition 4 . 2 .
Given the recurrent components M = {M (p) | p ∈ P}, Algorithm (1) computes all combinatorial attractor-repeller pairs of the form (A p , A * p ) | A p = Γ + (M (p)) and p ∈ P in time O(card(P)•E), where E is the number of edges in the directed graph corresponding to the combinatorial map F of Equation (5).

Proposition 5 . 5 .
Algorithm DistancePotential computes the combinatorial distance potential defined by Equation 17.The time and space complexity are O(dN + N log N ) and O(N ), respectively, where N = card(X ).Proof.Correctness of the combinatorial distance potential computation follows from the preceding discussion and the correctness of the computation of d W (A, A * ) follows from the observation that d W (A, A * ) = min ξ∈A d W (ξ, A * ).Since W will have approximately (d+1)N edges, the earlier discussion of Dijkstra's algorithm tells us that the time complexity will be O(dN + N log N ).The space complexity of O(N ) reflects the cost of storing the result, and storing the priority queues.
Trivial, Trivial) A trivial index-Conley index of the empty set (x − 1, Trivial, Trivial) Conley index of a stable fixed point (Trivial, x − 1, Trivial) Conley index of a saddle point (Trivial, Trivial, x − 1) Conley index of an unstable fixed point (x − 1, x − 1, Trivial) Conley index of a stable invariant circle (x T − 1, Trivial, Trivial) Conley index of a stable period-T orbit (Trivial, x T − 1, Trivial)

FIGURE 2 .
FIGURE 2. Conley-Morse Graph and a depiction of Morse Sets for the overcompensatory Leslie model of Equation (23) as computed by the SingleCMG program in the Conley-Morse-Database software package.

FIGURE 3 .
FIGURE 3. Combinatorial Lyapunov function for overcompensatory Leslie model of Equation (23) as computed by the Lyapunov program in the Conley-Morse-Database software package.The greyscale intensity represents the value of the computed combinatorial Lyapunov function with white corresponding to the highest vale and black the lowest value.

FIGURE 4 .FIGURE 5 .
FIGURE 4. Conley-Morse Graph and a depiction of Morse Sets for the Leslie/Gower Model of Equation (24) as computed by the SingleCMG program in the Conley-Morse-Database software package.

if end for Proposition 5.7. Algorithm LyapunovSummand correctly
(19)ies that cV A is a combinatorial Lyapunov function for (A, A * ).The following algorithm computes cV A .computes the combinatorial Lyapunov function cV A for the attractor-repeller pair (A, A * ) defined in Equation(19)in O(dN + N log N + E) time and O(N ) space, where N is the number of grid elements and