Eigenpairs for the analysis of complete Lyapunov functions

. A complete Lyapunov function describes the qualitative behaviour of a dynamical system: the areas where the orbital derivative vanishes and where it is strictly negative characterise the chain recurrent set and the gradient-like ﬂow, respectively. Moreover, its local maxima and minima show the stability properties of the connected components of the chain recurrent set. In this study, we use collocation with radial basis functions to numerically compute approximations to complete Lyapunov functions and then localise and analyse the stability properties of the connected components of the chain recurrent set using its gradient and Hessian. In particular, we improve the estimation of the chain recurrent set, and we determine the dimension and the stability properties of its connected components.


Introduction
We consider dynamical systems, whose dynamics are given by an autonomous system of ordinary differential equations (ODEs), x where f: R n ⟶ R n is a C 1 vector field. Complete Lyapunov functions fully classify the qualitative behaviour of the system. A complete Lyapunov function is a scalar-valued function V: R n ⟶ R, which is nonincreasing along solutions of (1). Moreover, it is strictly decreasing on a maximal set (with respect to set inclusion) [1][2][3]. If V is differentiable, then the decrease property can be expressed by V ′ (x) ≤ 0, where V ′ (x): � ∇V(x) · f(x) denotes the orbital derivative, that is, the derivative along solutions of (1). e chain recurrent set C is defined as the subset of R n , where infinitesimal perturbations make the dynamics recurrent, while the dynamics on its complement G are gradient-like. A complete Lyapunov function is constant on each chain transitive component of the chain recurrent set, including the compact connected components [4], and strictly negative where the flow is gradient-like. Recently in [5,6], the existence of a C ∞ complete Lyapunov function V for every ODE with a continuous right-hand side was established, such that C � x ∈ R n |∇V(x) � 0 { }. In this study, we advance the method from [7][8][9][10] to numerically compute approximations to complete Lyapunov functions and locate the chain recurrent set. Furthermore, we additionally introduce and test a novel algorithm to study the stability behaviour of the connected components of the chain recurrent set. e study is organised as follows. In the next section, we describe our numerical method to compute complete Lyapunov functions for system (1) using mesh-free collocation. In Section 3, we introduce a novel algorithm to analyse the stability of the chain transitive components of the chain recurrent set and describe an improved algorithm to localise the chain recurrent set. en, we apply our new methodology to numerous systems from the literature in Section 4 and analyse the results before concluding in Section 5.

Numerical Computation of Complete Lyapunov Functions
To compute approximations to complete Lyapunov functions for system (1), we use the method first presented in [7], but with various tweaks that improve its performance. e method uses collocation with radial basis functions to find a norm-minimal generalised interpolant v to the PDE in a reproducing kernel Hilbert space H. Typically, one sets r(x) � −1 for all x and localises the chain recurrent set by studying where the numerically computed approximative solution v to (2) fails to fulfill v ′ (x) � −1. Note that the PDE (2) is in general ill-posed. Indeed, the PDE does not have a solution since there cannot exist a function decreasing along solution trajectories in a domain containing a connected component of the chain recurrent set, for example, an equilibrium or a periodic orbit. e chosen radial basis function determines the kernel and thus the reproducing kernel Hilbert space. H is a reproducing kernel Hilbert space if there exists a (unique) kernel Φ: D × D ⟶ R such that Φ(·, x) ∈ H for all x ∈ D and g(x) � 〈g, Φ(·, x)〉 H for all x ∈ D and all g ∈ H, where 〈·, ·, 〉 H denotes the inner product of the Hilbert space H. Often, the kernel is defined by an appropriately chosen radial basis function ψ: R + 0 ⟶ R through Φ(x, y) � ψ(‖x − y‖). In [11], it was illustrated that Wendland's compactly supported radial basis functions [12] are particularly well suited for our application, and we use them in all our examples. e norm-minimal interpolant in reproducing kernel Hilbert spaces is a linear combination of the Riesz representers of functionals, and these are easy to compute in a reproducing kernel Hilbert space; we refer the interested reader to [13] for the general theory and [14] for the theory related to the computation of Lyapunov functions. e procedure to compute a complete Lyapunov function approximation v: D ⟶ R, D ⊂ R n compact, is as follows. First, fix a set of collocation points Note that a triangulation of the phase space is not needed. e quality of the approximation will depend on the socalled fill distance of X, that is, h X : � max x∈D min x i ∈X ‖x − x i ‖. Here and elsewhere in this study, ‖ · ‖ � ���� 〈·, ·〉 √ denotes the Euclidian norm. As collocation points, we use the hexagonal grid from [15], which delivers optimal fill distance for a given density of collocation points. e density is determined by a parameter α > 0; the smaller α > 0, the higher the density and the smaller the fill distance. As radial basis function ψ, we set ψ(x): � ψ l,k (c‖x‖), where ψ l,k is a Wendland function with smoothness parameter k ∈ N, l is fixed as l � ⌊n/2⌋ + k + 1, and c > 0 fixes the radius of the support of ψ as c − 1 . is results in a Hilbert space H that consists of the functions of the Sobolev space W k+((n+1)/2) 2 (R n ) , and its norm is equivalent to the usual norm of W k+((n+1)/2) 2 (R n ). By the Sobolev embedding theorems, W k+((n+1)/2) 2 (R n ) ⊂ C k (R n ). Hence, we take k ≥ 2 to assure that the Hessian of v is well defined. We use the C++ tool from [16] to generate ψ 0 � ψ and ψ i (r): � (1/r)(d/dr)ψ i−1 (r), i � 1, 2, 3, in optimised form for evaluation.
For continuous linear operators where the β i are chosen such that holds for all i � 1, 2, . . . , N, and is the unique norm-minimal function in H that satisfies the above conditions. e superscript y in L y i denotes the application of the operator L i to the function y↦ψ 0 (‖y − x‖). e coefficient vector β � (β i ) ∈ R N is determined by solving the system of equations Aβ � α, where the entries of A are given by a ij � L x i L y j ψ 0 (‖x − y‖) and the entries of α by α i � r(x i ). Note that A is symmetric and positive definite; hence, Aβ � α always has a solution, even if the PDE (2) is ill-posed.
Evaluating the L for v, and this function is our approximation to a complete Lyapunov function. More detailed explanations of this construction are given in ( [15], Chapter 3). We use two improvements of the method described above to compute v in some of the numerical examples. First, substituting the dynamical system (1) by with a small parameter δ > 0, often enhances the performance of the algorithm [17]. Note that systems (1) and (6) have the same trajectories [18], but the speed of the solutions to (6) is nearly uniform. Second, starting with r(x) � −1 in (2), one can use the numerical solution v 0 : � v to this problem to generate a new right-hand side r of (2) and solve again for v 1 : � v. Obviously, this procedure can be repeated, and generally, the approximation of the chain recurrent set improves with more iterations; see [9,19]. Note that since A has already been LU or Cholesky-factorised in the first step, subsequently solving the linear system Aβ � α is very fast. Using the numerical solution v j , we first set . . , N and then normalise the values through r( and Y x i is a set of points aligned along the flow of the ODE system at the collocation point x i . In particular, we use where the parameter ρ ∈ (0, 1) determines how far the evaluation points will be placed away from x i and m ∈ N denotes the number of points on each side of the collocation point; recall that α > 0 controls the density of the points in the collocation grid X and thus the distance between the collocation points.

The New Method
In this section, we first discuss how the new method is used for the stability analysis and then for the improved determination of the chain recurrent set.
Let us first discuss our method for the stability analysis of the points in the chain recurrent set C. In the following, V denotes a true complete Lyapunov function and v a numerical approximation to V. It is based on the following observations. For a true complete Lyapunov function V for system (1) as constructed in [6], we have where H V (x) is the Hessian matrix of V at x and (‖O(‖h‖)‖/‖h‖ 3 ) is bounded as h ⟶ 0. Furthermore, if x ∈ C, we have ∇V(x) � 0 and Since the Hessian matrix H V (x) is symmetric, its eigenvalues λ i are real numbers, and there is an orthogonal set of corresponding eigenvectors v i . us, for an eigenvector v i , the function t↦V(x + tv i ) has a local strict minimum at t � 0 if λ i > 0 and a local strict maximum if has a root of multiplicity 2 or higher, that is, is implies that we can expect for all points x in a connected component of C that the number of negative, positive, and zero eigenvalues of the Hessian of V is constant and correspond to the dimensions of the unstable, stable, and centre manifolds, respectively.
In our approach to gain information about the qualitative behaviour of the dynamical system defined by (1), we use a numerical method to compute an approximate complete Lyapunov function v: D ⟶ R on a compact set D. Our hypothesis is that the stability analysis, described above for a true complete Lyapunov function, also holds sufficiently accurately for an approximate one for us to gain information on the connected components of the chain recurrent set. is hypothesis will be tested Complexity in the examples. e numerical method we use, described in Section 2, delivers an explicit formula for v in terms of radial basis functions, the smoothness of which can be chosen arbitrarily high.
us, for x ∈ D and sufficiently small h ∈ R n , assuming v suitably approximates a true complete Lyapunov function V as above, we have for a point x in the chain recurrent set that ∇v(x) ≈ 0 and Now, we expect an eigenvector v i of H v (x) for an eigenvalue λ i ≈ 0 to be tangent to the connected component of the chain recurrent set, that is, a set along which v ′ is constant.
For an eigenvector v i of H v (x) to an eigenvalue λ i < 0 considerably smaller than zero, we expect t↦v(x + tv i ) to have a strict maximum and the dynamics to repel nearby solution trajectories in direction ± v i because v should be decreasing along solution trajectories.
Similarly, for an eigenvector v i of H v (x) to an eigenvalue λ i > 0 considerably larger than zero, we expect t↦v(x + tv i ) to have a strict minimum and the dynamics to attract nearby solution trajectories in direction ± v i . In all our examples, we used the criteria and from the results, it is clear that this classification is not sensitive to the exact value 3 because the values of the eigenvalues are clearly separated. e directions of the eigenvectors corresponding to positive/negative eigenvalues thus show the directions of the local stable/unstable manifold, respectively, while the directions of the eigenvectors corresponding to zero eigenvectors are tangential to the chain recurrent set. Now we will describe how we can use the eigenvalues to improve the determination of the chain recurrent set. In [10], the chain recurrent set of an ODE was localised by computing complete Lyapunov functions as in the last section and then using the criterion ‖∇v(x)‖ ≈ 0 for points x in the chain recurrent set. In practice, this was implemented by using a finite set of evaluation points X eval : � hZ n ∩ D, whose density is determined through the parameter h > 0. We fix a parameter c + > 0, for example, c + � 0.6, and define the approximation X 0 ⊂ X eval to the chain recurrent set C through x ∈ X 0 , if and only if us, C v,c + is an approximation of the chain recurrent set using the gradient of v and the threshold value c + . is method to approximate the chain recurrent set will in general either overestimate or underestimate the chain recurrent set. An underestimation is more problematic than an overestimation because the chain transitive components are then not necessarily connected. Note that the set X eval is discrete, and when we say that the components are not connected, we mean that there are points in X eval \X 0 between points in X 0 that belong to the same chain transitive component. Furthermore, for our evaluation grid X eval : � hZ n ∩ D, the classification of the connected components of X 0 boils down to x, y belonging to the same In our new approach, we first use the gradient of v to get a rough estimate X 0 : � C v,λ + of the chain recurrent set using the condition ‖∇v(x)‖ ≤ λ + , where we choose the parameter λ + so large that we get a considerable overestimation. Indeed, in all our examples, we set λ + � 14. en, we use a clustering algorithm to subdivide X 0 into connected components X 0 k , k � 1, 2, . . . , K. Subsequently, we fix a small c + ; in the examples, we used the values 0.6, 0.7, and 0.8 and compute the set C v,c + . Furthermore, we define Note that each X 0 k is an overestimation of a connected component of the chain recurrent set and C k v,c + ⊂ X 0 k an underestimation. at is, all or most points in C k v,c + are truly in the k th component, and all or most points in the k th component are also in X 0 k . We use C k v,c + to classify the stability properties of the k th component by analysing the Hessian matrix H v (x) of v at the points x ∈ C k v,c + . As we consider an approximation v of a true complete Lyapunov function, these over-and underestimations are also approximate. e classification is done by counting the negative, positive, and zero eigenvalues of the Hessian matrices. Finally, we obtain an improved approximation X 0 k of the k th component of the chain recurrent set, as the set of points x ∈ X 0 k , where either x ∈ C k v,c + or its Hessian matrix H v (x) is of the same type as at the points in C k v,c + . Note that although the analysis of maxima and minima of V using the gradient and the Hessian, described at the beginning of this section, is elementary, it is not clear how well it works in practice for our numerically computed approximation v to a complete Lyapunov function, in particular because our constructed function is the solution to a discretisation of an ill-posed PDE that does not have a solution. We identify the chain recurrent set as the area where the numerical solution v fails badly at fulfilling the ill-posed PDE in between the collocation points, and it is not a priori clear that the second-order derivatives of v carry enough information to classify the connected components of the chain recurrent set. However, in the next section, we investigate several systems with our new method and will see that it works very well in practice to estimate the connected components of the chain recurrent set and classify their stability properties.

Examples
As discussed in the introduction, we first estimate the chain recurrent set of a system given by an ODE (1), and then we compute the Hessian of a complete Lyapunov function approximation at points in our approximation of the chain recurrent set. From the eigenvalues and the corresponding eigenvectors, we seek to obtain information about attractive, repelling, and neutral directions.
Here, we apply our novel algorithm to various systems. First, an approximation v to a complete Lyapunov function is computed as described in Section 2. en, an overestimation X 0 : � x ∈ X eval |‖v(x)‖ ≤ λ + with λ + � 14 is computed and subdivided into its connected components X 0 k . en, the underestimations where c + ∈ [0.6, 0.8], depending on the example, the connected components of the chain recurrent set are computed, and the eigenvalues of the Hessian at the points in C k v,c + are used to classify the stability behaviour of the k th component. For the visual presentation, it is convenient to plot the eigendirections in red (unstable, negative eigenvalue), green (stable, positive eigenvalue), and black (neutral, zero eigenvalue), indicating the local stability of the connected component of the chain recurrent set. In a final step, improved estimations X 0 k , , of the components of the chain recurrent set are computed as described in Section 3.
In all our computations, we use the Wendland function ψ � ψ 6,4 (cr) with c � 0.5, resulting in 2 � (1/c) for the support radius of the Wendland function, and when we used the normalised ODE in (4), we set δ 2 � 10 − 8 . Furthermore, we choose the fineness parameter α > 0 for the density of the collocation such that the condition number κ 2 (A) of the collocation matrix A is ca. 10 12 . We always started with r(x) � −1 in the PDE (2); subsequent iterations were performed as described in Section 2 with m � 10 and ρ � 0.5 in formula (5). e best threshold parameter c + to estimate the chain recurrent set varies between 0.6 and 0.8 for the different examples. is is consistent with our findings in [20]. We solve the linear systems Aβ � α using (dense matrix) Cholesky factorisation. For larger systems in higher dimensions, one should use faster methods, taking advantage of the sparsity of A. For our examples, however, this is not necessary.
Since our method detects the connected components of the chain recurrent set, we categorize our examples using the same criterion.

One Periodic Orbit, One
Equilibrium.
e first system we considered was is system has one unstable equilibrium at (0, 0) and an attractive periodic orbit at Ω � (x, y) ∈ R 2 |x 2 + y 2 � 1 .
We set the fineness parameter α � 0.04 for the hexagonal collocation grid, and we used 10 iterations (Section 2). e computations were performed using f directly, rather than the normalised version in formula (4). e domain used was D � [−1.2, 1.2] 2 ⊂ R 2 , and to evaluate the function, we used the regular grid X eval � hZ ∩ D with parameter h � 0.008. Figure 1 shows the computed complete Lyapunov function v for system (18), together with its orbital derivative v ′ , the norm of the gradient ‖∇v‖, and the underestimation C v,0.8 of the chain recurrent set. e behaviour of the system is clearly represented: the repeller at the origin is a local maximum of the complete Lyapunov function, while the attractive periodic orbit is a local minimum. Recall that a complete Lyapunov function is decreasing along solution trajectories. e Figure 1(b) figure shows the orbital derivative: at most of the points, its value is close to −1 as a consequence of the collocation method. e points on the periodic orbit and around the critical point have values for the orbital derivative close to zero. e norm of the gradient is close to zero as well on a similar set; note the overshoot close to the area where it is zero. Figure 2 shows the eigenvalues and eigenvectors of the Hessian at the points of C v,0.8 . e upper figures show the eigenvalues, which are positive (green) and zero (black) on the periodic orbit. is indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). e lower figure shows the corresponding eigendirections (line segments parallel to the eigenvectors): as expected, the eigendirections corresponding to the zero eigenvalues are tangential to the chain recurrent set, while the ones corresponding to the positive (attracting) one are perpendicular and indicate the direction of attraction.
At the equilibrium point, we see two negative (red) eigenvalues, indicating a 2-dimensional unstable manifold. e lower figure shows the corresponding eigendirections in the same colours as the eigenvalues; recall that they are perpendicular because the Hessian is symmetric. Note that, at the equilibrium, several different directions can be seen because the approximation of the chain recurrent set covers a slightly larger set than just the equilibrium and thus consists of several points with corresponding eigendirections.
In Figure 3, we plot the overestimation X 0 � C v,14 of the chain recurrent set Figure 3(a) and the underestimation C v,0.8 Figure 3(b) and our refined approximation X 0 Figure 3(c), where we keep all points of the underestimation and add points of the overestimation if they have the same number of negative, positive, and zero eigenvalues as the underestimation (Section 3). For comparison, we zoom into a smaller area and compare the overestimation (black) with the refined approximation (green) in the Figure 3(d). e refined approximation is clearly superior to both the overestimated and the underestimated set. Finally, in Figure 4, we plot the eigenvalues of the Hessian at the points in X 0 close to the origin. At all points, both are negative, demonstrating that there is a 2-dimensional unstable manifold at the origin. We observe that the minimal value of the negative eigenvalues is obtained at the equilibrium.

Two Periodic Orbits, One
Equilibrium.
For our computations, we used the normalisation (6), the domain D � [−1.2, 1.2] 2 ⊂ R 2 , and the parameter values α � 0.04 and h � 0.008. We used 10 iterations (Section 2). Figure 5 shows the computed complete Lyapunov function v for system (19), together with its orbital derivative v ′ , the norm of its gradient ‖∇v‖, and the underestimation C v,0.8 of the chain recurrent set. We clearly see the qualitative behaviour of the system: the local minimum at the origin is an attracting equilibrium, Ω 2 (circle with radius 1/2, local maximum) is a repelling periodic orbit, and Ω 1 (circle with radius 1, local minimum) is an attracting periodic orbit.
While one can determine the properties of the chain recurrent set by checking local minima and maxima visually, we now use the eigenvalues and eigenvectors of the Hessian to determine these properties automatically. is indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). ey are negative (red) and zero (black) on the inner periodic orbit, indicating a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is repelling (the remaining one is negative). Finally, the equilibrium has two positive (green) eigenvalues, indicating stable behaviour (2-dimensional stable manifold). e lower figure shows the corresponding eigendirections: as expected, the eigendirections corresponding to the zero eigenvalues are tangential to the periodic orbits, while the others are perpendicular. For the outer orbit, they correspond to the attracting (positive) direction, while for the inner orbit, they show the repelling (negative) direction. At the equilibrium point, we see the corresponding attracting (green) eigendirections.
In Figure 7(a), we plot the overestimation X 0 � C v,14 of the chain recurrent set (black) and our refined approximation X 0 (green), computed as described in Section 3.
Finally, in Figure 7(b), we plot the eigenvalues of the Hessian at the points in X 0 close to the origin. At all points, both are positive, demonstrating that there is a 2-dimensional stable manifold. We observe that the maximal eigenvalue is obtained at the equilibrium.

Van der Pol.
Next, we consider the classical van der Pol system It has an unstable equilibrium at the origin and an attractive periodic orbit.
For our computations, we used the normalisation (6), the domain D � [−4, 4] 2 ⊂ R 2 , and the parameter values α � 0.044, and h � 0.008. We used 10 iterations (Section 2). Figure 8 shows the computed complete Lyapunov function v for the van der Pol system (20), together with its orbital derivative v ′ and the norm of its gradient ‖∇v‖, and the underestimation C v,0.6 of the chain recurrent set. is clearly shows the behaviour of the system: the local maximum at the origin corresponds to a 2-dimensional unstable manifold, while the periodic orbit is a local minimum and thus attracting. In Figure 10(a), we plot the overestimation X 0 � C v,14 of the chain recurrent set (black) and our refined approximation X 0 (green), computed as described in Section 3. Note that there are small gaps in the refined approximation close to (−1.5, 0.5) and (1.5, −1), that is, exactly where the magnitude of the positive eigenvalue is varying drastically, as shown in Figure 9. Finally, in Figure 10(b), we plot the eigenvalues of the Hessian at the  points in X 0 close to the origin. At all points, both are negative, demonstrating that there is a two-dimensional unstable manifold at the origin.
For a three-dimensional system, we cannot plot the complete Lyapunov function v as a function of its arguments, and hence we only plot the eigendirections of the Hessian on the underestimation C v,0.7 of the chain recurrent set. In particular, we set c + � 0.7. Figure 11 shows the connected components of the chain recurrent set C v,0.7 close to z � −1 in the xy-plane Figure 11(a) and the xz-plane Figure 11(b). e periodic orbit Ω −1 has a tangential eigendirection corresponding to a zero eigenvalue (black) and a 2-dimensional unstable manifold (red) in the z-direction and towards the equilibrium (0, 0, −1). e equilibrium (0, 0, −1) has a stable 2-dimensional manifold (green) in the xy-plane Figure 11(c) and a 1-dimensional unstable manifold (red)   Figure 13 shows the connected components of the chain recurrent set close to z � 0 in the xy-plane (Figure 13(a)) and the xz-plane (Figure 13(b)). e periodic orbit Ω 0 has a tangential eigendirection corresponding to zero eigenvalues (black), a 1-dimensional unstable manifold (red) towards the equilibrium (0, 0, 0), and a 1-dimensional stable manifold (green) in z-direction. e equilibrium (0, 0, 0) has a stable 3-dimensional manifold (green). Figure 14 shows the overestimation X 0 � C v,14 of the chain recurrent set (black) close to z � 1. Note the much    Complexity finer scale of the z-axis; the four different layers are connected. Figure 14(a) clearly shows parts of the overestimation (the triangles) that do not correspond to a component of the chain recurrent set. e underestimation C v,0.7 is disjoint from all these connected components of X 0 , and therefore, these components are removed in our improved approximation X 0 of the chain recurrent set. We only show the estimate close to the orbit at z � 1; the results close to z � 0 and z � −1 are similar and are omitted.

Analysis of the Results.
As can be seen from the examples, the new algorithm to localise the chain recurrent set and its connected components works very well and    compares favourably to earlier approaches in [7][8][9][10]. e only minor artefact we observed were the small gaps in the refined approximation for the van der Pol system, where the magnitude of the positive eigenvalue is drastically changing; this should and will be thoroughly analysed in future work. As an added bonus, the new method also computes algorithmically the dimensions of the stable, unstable, and centre manifold of the connected components of the chain recurrent set. While these stability properties of the components can be visually observed in two and three-dimensional systems, this becomes impossible for higher-dimensional systems.
However, the method described in this study can easily compute this information from the eigenvalues in any dimension. us, our new method is essential when computing complete Lyapunov functions in higher dimensions.
Finally, the new algorithm for computing a refined estimate of the connected components of the chain recurrent set and analysing their stability properties does not add to the computational complexity of the previous algorithm, analysed in [19], as it only requires the computation of the Hessian matrix on the set X 0 , which is usually much smaller than the evaluation grid X eval .  Figure 11: System (21). z � −1, periodic orbit Ω −1 and equilibrium at (0, 0, −1) in the xy-plane (a) and in the xz-plane (b). Ω −1 has a twodimensional unstable manifold (in z-direction and towards the equilibrium). Equilibrium (0, 0, −1) in the xy-plane (c) and in the xz-plane (d). Two-dimensional stable manifold (in xy-plane) and one-dimensional unstable manifold (in z-direction).  Figure 12: System (21). z � 1, periodic orbit Ω 1 and equilibrium at (0, 0, 1) in the xy-plane (a) and in the xz-plane (b). Ω 1 has a twodimensional unstable manifold (in z-direction and towards the equilibrium). Equilibrium (0, 0, 1) in the xy-plane (c) and in the xz-plane (d). Two-dimensional stable manifold (in xy-plane) and one-dimensional unstable manifold (in z-direction).  Figure 13: System (21). z � 0, periodic orbit Ω 0 and equilibrium at (0, 0, 0) in the xy-plane (a) and in the xz-plane (b). Ω 0 has a onedimensional unstable manifold (towards the equilibrium) and one-dimensional stable manifold (in z-direction). Equilibrium (0, 0, 0) in the xy-plane (c) and in the xz-plane (d) with a three-dimensional stable manifold.

Conclusions
We have extended an existing construction method for computing complete Lyapunov functions for nonlinear systems given by ODEs. Our new method considerably improves the approximation of the chain recurrent set of the ODE. Furthermore, it algorithmically subdivides it into connected components and analyses their stability properties. is enables the determination of the dimension of the components and the local directions of the stable and unstable manifolds. ese properties are derived from the eigenvalues and eigenvectors of the Hessian of the computed complete Lyapunov function. An implementation of our novel method will be included in a future freely distributed version of [21].
Data Availability e code used to perform all the computations is LyapXool-V2, and is available at https://github.com/LyapXool. e most general data to reproduce the figures is available at https://www.dropbox.com/sh/gl0l2h3aqxuagj5/AACUsd-BIjlFUGN6eZsmA9qQza?dl=0. 16 Complexity