# Eigenpairs for the Analysis of Complete Lyapunov Functions

Authors: Carlos Argáez; Peter Giesl; Sigurdur Freyr Hafstein
Journal: Complexity
 (2022)

---

## Abstract

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 flow, 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.


---

## Body

## 1. Introduction

We consider dynamical systems, whose dynamics are given by an autonomous system of ordinary differential equations (ODEs),(1)x.=fx,x∈ℝn,where f:ℝn⟶ℝn is a C1 vector field. Complete Lyapunov functions fully classify the qualitative behaviour of the system. A complete Lyapunov function is a scalar-valued function V:ℝn⟶ℝ, which is nonincreasing along solutions of (1). Moreover, it is strictly decreasing on a maximal set (with respect to set inclusion) [1–3]. If V is differentiable, then the decrease property can be expressed by V′x≤0, where V′x:=∇Vx⋅fx denotes the orbital derivative, that is, the derivative along solutions of (1). The chain recurrent set C is defined as the subset of ℝ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∈ℝn|∇Vx=0.In this study, we advance the method from [7–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.The 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. Then, we apply our new methodology to numerous systems from the literature in Section 4 and analyse the results before concluding in Section 5.


## 2. 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. The method uses collocation with radial basis functions to find a norm-minimal generalised interpolant v to the PDE(2)V′x=∑i=1n∂V∂xixfix=rx,wherer:ℝn⟶−∞,0,in a reproducing kernel Hilbert space ℋ. Typically, one sets rx=−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. The chosen radial basis function determines the kernel and thus the reproducing kernel Hilbert space. ℋ is a reproducing kernel Hilbert space if there exists a (unique) kernel Φ:D×D⟶ℝ such that Φ⋅,x∈ℋ for all x∈D and gx=g,Φ⋅,xℋ for all x∈D and all g∈ℋ, where ⋅,⋅,ℋ denotes the inner product of the Hilbert space ℋ. Often, the kernel is defined by an appropriately chosen radial basis function ψ:ℝ0+⟶ℝ 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. The 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.The procedure to compute a complete Lyapunov function approximationv:D⟶ℝ, D⊂ℝn compact, is as follows. First, fix a set of collocation points X=x1,x2,…,xN⊂D, where fxi≠0 and xi≠xj if i≠j. Note that a triangulation of the phase space is not needed. The quality of the approximation will depend on the so-called fill distance of X, that is, hX:=maxx∈Dminxi∈Xx−xi. 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. The 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,kcx, where ψl,k is a Wendland function with smoothness parameter k∈ℕ, l is fixed as l=⌊n/2⌋+k+1, and c>0 fixes the radius of the support of ψ as c−1. This results in a Hilbert space ℋ that consists of the functions of the Sobolev space W2k+n+1/2ℝn , and its norm is equivalent to the usual norm of W2k+n+1/2ℝn. By the Sobolev embedding theorems, W2k+n+1/2ℝn⊂Ckℝ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 ψir:=1/rd/drψi−1r, i=1,2,3, in optimised form for evaluation.For continuous linear operatorsLi:ℋ⟶ℝ, in our case Liv:=v′xi, one can show that(3)vx=∑i=1NβiLiyψ0y−x,where the βi are chosen such that(4)Liv=rxiholds for all i=1,2,…,N, and is the unique norm-minimal function in ℋ that satisfies the above conditions. The superscript y in Liy denotes the application of the operator Li to the function y↦ψ0y−x. The coefficient vector β=βi∈ℝN is determined by solving the system of equations Aβ=α, where the entries of A are given by aij=LixLjyψ0x−y and the entries of α by αi=rxi. Note that A is symmetric and positive definite; hence, Aβ=α always has a solution, even if the PDE (2) is ill-posed.Evaluating theLiyψ0 delivers the explicit formula(5)vx=∑k=1Nβkxk−xTfxkψ1x−xkfor 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 computev in some of the numerical examples. First, substituting the dynamical system (1) by(6)x.=f^x,wheref^x=fxδ2+fx2,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 withrx=−1 in (2), one can use the numerical solution v0:=v to this problem to generate a new right-hand side r of (2) and solve again for v1:=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 vj, we first set r˜xi=1/Yxi∑y∈Yxivj′y− for i=1,2,…,N and then normalise the values through rxi=N/∑l=1Nr˜xlr˜xi. Here, x−:=minx,0, and Yxi is a set of points aligned along the flow of the ODE system at the collocation point xi. In particular, we use(7)Yxi=xi±ρ⋅k⋅α⋅f^xjm|k∈1,…,m,where the parameter ρ∈0,1 determines how far the evaluation points will be placed away from xi and m∈ℕ 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.


### 2.1. Formulas

In this section, we compute the first and second-order derivatives ofv in formula (3) and the formulas for the Wendland functions:(8)∂v∂xix=∑k=1Nβk∂∂xixk−xTfxkψ1x−xk=∑k=1Nβk−fixkψ1x−xk+x−xkixk−xTfxkψ2x−xk(9)∂2v∂xi∂xjx=∑k=1Nβk∂∂xj−fixkψ1x−xk+x−xkixk−xTfxkψ2x−xk=∑k=1Nβk−x−xkjfixkψ2x−xk+δijxk−xTfxkψ2x−xk−x−xkifjxkψ2x−xk+x−xkix−xkjxk−xTfxkψ3x−xkwhere δij denotes the Kronecker delta and yi:=yi is the ith component of the vector y. These formulas can be used to compute the gradient ∇vx and the Hessian Hvx of v at arbitrary points x∈ℝn. The formula for the orbital derivative of v with respect to the ODE (1) is given by(10)v′x=∑k=1Nβk−ψ1x−xkfx,fxk+ψ2x−xkx−xk,fx⋅xk−x,fxk.The Wendland functionψr=ψ6,4cr was used in all the examples. The formulas for ψ0=ψ, ψ1, ψ2, and ψ3 defined above are(11)ψ0r=1−cr10429cr4+450cr3+210cr2+50cr+5,ψ1r=−26c21−cr95+45cr+159cr2+231cr3,ψ2r=3432c41−cr81+8cr+21cr2,ψ3r=−102960c61−cr71+7cr,for 0≤cr≤1 and otherwise equal to zero.


## 2.1. Formulas

In this section, we compute the first and second-order derivatives ofv in formula (3) and the formulas for the Wendland functions:(8)∂v∂xix=∑k=1Nβk∂∂xixk−xTfxkψ1x−xk=∑k=1Nβk−fixkψ1x−xk+x−xkixk−xTfxkψ2x−xk(9)∂2v∂xi∂xjx=∑k=1Nβk∂∂xj−fixkψ1x−xk+x−xkixk−xTfxkψ2x−xk=∑k=1Nβk−x−xkjfixkψ2x−xk+δijxk−xTfxkψ2x−xk−x−xkifjxkψ2x−xk+x−xkix−xkjxk−xTfxkψ3x−xkwhere δij denotes the Kronecker delta and yi:=yi is the ith component of the vector y. These formulas can be used to compute the gradient ∇vx and the Hessian Hvx of v at arbitrary points x∈ℝn. The formula for the orbital derivative of v with respect to the ODE (1) is given by(10)v′x=∑k=1Nβk−ψ1x−xkfx,fxk+ψ2x−xkx−xk,fx⋅xk−x,fxk.The Wendland functionψr=ψ6,4cr was used in all the examples. The formulas for ψ0=ψ, ψ1, ψ2, and ψ3 defined above are(11)ψ0r=1−cr10429cr4+450cr3+210cr2+50cr+5,ψ1r=−26c21−cr95+45cr+159cr2+231cr3,ψ2r=3432c41−cr81+8cr+21cr2,ψ3r=−102960c61−cr71+7cr,for 0≤cr≤1 and otherwise equal to zero.


## 3. 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 setC. 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(12)Vx+h=Vx+∇Vx⋅h+12hTHVxh+Oh3,where HVx is the Hessian matrix of V at x and Oh/h3 is bounded as h⟶0. Furthermore, if x∈C, we have ∇Vx=0 and(13)Vx+h−Vx=12hTHVxh+Oh3.Since the Hessian matrixHVx is symmetric, its eigenvalues λi are real numbers, and there is an orthogonal set of corresponding eigenvectors vi. Thus, for an eigenvector vi, the function t↦Vx+tvi has a local strict minimum at t=0 if λi>0 and a local strict maximum if λi<0. If t↦Vx+tvi is constant in a neighbourhood of t=0, we necessarily have λi=0. On the other hand, if λi=0, then t↦Vx+tvi−Vx has a root of multiplicity 2 or higher, that is, t↦Vx+tvi is “flat” at t=0. This 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⟶ℝ 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. This hypothesis will be tested in the examples. The 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. Thus, for x∈D and sufficiently small h∈ℝ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 ∇vx≈0 and(14)vx+h−vx≈12hTHvxh.Now, we expect an eigenvectorvi of Hvx 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 eigenvectorvi of Hvx to an eigenvalue λi<0 considerably smaller than zero, we expect t↦vx+tvi to have a strict maximum and the dynamics to repel nearby solution trajectories in direction ±vi because v should be decreasing along solution trajectories.Similarly, for an eigenvectorvi of Hvx to an eigenvalue λi>0 considerably larger than zero, we expect t↦vx+tvi to have a strict minimum and the dynamics to attract nearby solution trajectories in direction ±vi. In all our examples, we used the criteria(15)λ≈0ifλ≤3,λis positive ifλ>3,λis negative ifλ<−3,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.The 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 ∇vx≈0 for points x in the chain recurrent set. In practice, this was implemented by using a finite set of evaluation points Xeval:=hℤn∩D, whose density is determined through the parameter h>0. We fix a parameter γ+>0, for example, γ+=0.6, and define the approximation X˜0⊂Xeval to the chain recurrent set C through x∈X˜0, if and only if ∇vx≤γ+. It is convenient to define the set Cv,γ+⊂Xeval through(16)x∈Cv,γ+:⇔x∈Xeval,∇vx≤γ+.Thus,Cv,γ+ is an approximation of the chain recurrent set using the gradient of v and the threshold value γ+.This 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 setXeval is discrete, and when we say that the components are not connected, we mean that there are points in Xeval\X˜0 between points in X˜0 that belong to the same chain transitive component. Furthermore, for our evaluation grid Xeval:=hℤn∩D, the classification of the connected components of X˜0 boils down to x,y belonging to the same component X˜k0, if and only if there exist points x1,x2,…,xj∈X˜k0 with x1=x, xj=y, and xi+1−xi≤hn for i=1,2,…,j−1.In our new approach, we first use the gradient ofv to get a rough estimate X˜0:=Cv,λ+ of the chain recurrent set using the condition ∇vx≤λ+, where we choose the parameter λ+ so large that we get a considerable overestimation. Indeed, in all our examples, we set λ+=14. Then, we use a clustering algorithm to subdivide X˜0 into connected components X˜k0, k=1,2,…,K. Subsequently, we fix a small γ+; in the examples, we used the values 0.6, 0.7, and 0.8 and compute the set Cv,γ+. Furthermore, we define(17)Cv,γ+k≔Cv,γ+∩X˜k0.Note that eachX˜k0 is an overestimation of a connected component of the chain recurrent set and Cv,γ+k⊂X˜k0 an underestimation. That is, all or most points in Cv,γ+k are truly in the kth component, and all or most points in the kth component are also in X˜k0.We useCv,γ+k to classify the stability properties of the kth component by analysing the Hessian matrix Hvx of v at the points x∈Cv,γ+k. As we consider an approximation v of a true complete Lyapunov function, these over- and underestimations are also approximate. The classification is done by counting the negative, positive, and zero eigenvalues of the Hessian matrices. Finally, we obtain an improved approximation Xk0 of the kth component of the chain recurrent set, as the set of points x∈X˜k0, where either x∈Cv,γ+k or its Hessian matrix Hvx is of the same type as at the points in Cv,γ+k.Note that although the analysis of maxima and minima ofV 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.


## 4. 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 approximationv to a complete Lyapunov function is computed as described in Section 2. Then, an overestimation X˜0:=x∈Xeval|vx≤λ+ with λ+=14 is computed and subdivided into its connected components X˜k0. Then, the underestimations Cv,γ+k=x∈X˜k0|vx≤γ+, where γ+∈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 Cv,γ+k are used to classify the stability behaviour of the kth 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 Xk0, Cv,γ+k⊂Xk0⊂X˜k0, 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,4cr 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 κ2A of the collocation matrix A is ca. 1012. We always started with rx=−1 in the PDE (2); subsequent iterations were performed as described in Section 2 with m=10 and ρ=0.5 in formula (5). The best threshold parameter γ+ to estimate the chain recurrent set varies between 0.6 and 0.8 for the different examples. This 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.


### 4.1. One Periodic Orbit, One Equilibrium

The first system we considered was(18)x˙=−xx2+y2−1−y,y˙=−yx2+y2−1+x.This system has one unstable equilibrium at0,0 and an attractive periodic orbit at Ω=x,y∈ℝ2|x2+y2=1.We set the fineness parameterα=0.04 for the hexagonal collocation grid, and we used 10 iterations (Section 2). The computations were performed using f directly, rather than the normalised version in formula (4). The domain used was D=−1.2,1.22⊂ℝ2 , and to evaluate the function, we used the regular grid Xeval=hℤ∩D with parameter h=0.008.Figure1 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 Cv,0.8 of the chain recurrent set.Figure 1
System (18). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.8.
(a)(b)(c)(d)The 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. The Figure1(b) figure shows the orbital derivative: at most of the points, its value is close to −1 as a consequence of the collocation method. The points on the periodic orbit and around the critical point have values for the orbital derivative close to zero. The 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.Figure2 shows the eigenvalues and eigenvectors of the Hessian at the points of Cv,0.8. The upper figures show the eigenvalues, which are positive (green) and zero (black) on the periodic orbit. This indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). The 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.Figure 2
System (18). (a) Eigenvalues of the Hessian at points of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected component of the chain recurrent set. The periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, while the equilibrium has a 2-dimensional unstable manifold.
(a)(b)At the equilibrium point, we see two negative (red) eigenvalues, indicating a 2-dimensional unstable manifold. The 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 Figure3, we plot the overestimation X˜0=Cv,14 of the chain recurrent set Figure 3(a) and the underestimation Cv,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). The refined approximation is clearly superior to both the overestimated and the underestimated set.Figure 3
System (18). (a) The overestimation of the chain recurrent set X˜0=Cv,14. (b) The underestimation of the chain recurrent set Cv,0.8. (c) The refined approximation with our new method X˜0 of the chain recurrent set. (d) Comparison between the overestimated chain recurrent set (in black) and the reduced chain recurrent set (green) in a zoomed area.
(a)(b)(c)(d)Finally, in Figure4, 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.Figure 4
System (18). Eigenvalues of the Hessian at the points in X˜0 close to the origin. Since both are negative, v must have a maximum at the critical point 0,0.


### 4.2. Two Periodic Orbits, One Equilibrium

The second system we consider is(19)x˙=−xx2+y2−1x2+y2−14−y,y˙=−yx2+y2−1x2+y2−14+x.This system has a stable equilibrium at0,0, an attractive periodic orbit Ω1=x,y∈ℝ2|x2+y2=1, and a repelling periodic orbit Ω2=x,y∈ℝ2|x2+y2=1/4.For our computations, we used the normalisation (6), the domain D=−1.2,1.22⊂ℝ2, and the parameter values α=0.04 and h=0.008. We used 10 iterations (Section 2).Figure5 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 Cv,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.Figure 5
System (19). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.8.
(a)(b)(c)(d)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.Figure6 shows the eigenvalues and eigenvectors of the Hessian at points of the underestimation Cv,0.8 of the chain recurrent set. The upper figures show the eigenvalues: the eigenvalues are positive (green) and zero (black) on the outer periodic orbit. This indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). They 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).Figure 6
System (19). (a) Eigenvalues of the Hessian at points of the underestimation Cv,0.8 of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected components of the chain recurrent set. The outer periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, and the inner periodic orbit has a 1-dimensional neutral direction and a 1-dimensional unstable one, while the equilibrium has a 2-dimensional stable manifold.
(a)(b)The 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 Figure7(a), we plot the overestimation X˜0=Cv,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.Figure 7
System (19). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. Our refined approximation X˜0 of the chain recurrent set is plotted in green. It includes the underestimation Cv,0.8 and points from the overestimation X˜0 at the inner orbit if one eigenvalue is negative and one zero. At the outer orbit, it includes points from X˜0 if one eigenvalue is negative and one is zero. (b) Eigenvalues of the points in X˜0 close to the critical point at the origin. All are positive, so the equilibrium has a 2-dimensional stable manifold.
(a)(b)


### 4.3. Van der Pol

Next, we consider the classical van der Pol system(20)x˙=y,y˙=1−x2y−x.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,42⊂ℝ2, and the parameter values α=0.044, and h=0.008. We used 10 iterations (Section 2).Figure8 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 Cv,0.6 of the chain recurrent set. This 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.Figure 8
System (20). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.6.
(a)(b)(c)(d)Figure9 shows the eigenvalues and eigenvectors of the Hessian at points of the chain recurrent set. The top figures show the eigenvalues, which are positive (green) and zero (black) on the periodic orbit and both negative at the equilibrium (red). The lower figure shows the corresponding eigendirections: the eigenvectors corresponding to the zero eigenvalues are tangential to the periodic orbit, while the ones corresponding to the positive (attracting) one are perpendicular to it. At the equilibrium point, we see the corresponding repelling eigendirections.Figure 9
System (20) : (a) Eigenvalues of the Hessian at points of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected component of the chain recurrent set. The periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, while the equilibrium has a 2-dimensional unstable manifold.
(a)(b)In Figure10(a), we plot the overestimation X˜0=Cv,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.Figure 10
System (20). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. Our refined approximation X0 of the chain recurrent set is plotted in green. We include points if there is one positive and one zero eigenvalues. (b) Eigenvalues of the points in X˜0 close to the critical point at the origin. Both are negative, so v must have a maximum, that is, the critical point is unstable with a 2-dimensional unstable manifold.
(a)(b)


### 4.4. Three-Dimensional Case

Finally, we consider the three-dimensional system(21)x˙=xx2+y2−1−xz2+1,y˙=yx2+y2−1+xz2+1,z˙=10zz2−1.This system has three equilibria at0,0,±1 and 0,0,0 as well as three periodic orbits Ωz0=x,y,z∈ℝ3|x2+y2=1,z=z0 for z0∈0,±1.For our computations, we used the normalisation (6), the domain D=−1.2,1.23⊂ℝ3, and the parameter values α=0.07 and h=0.015; note that this is a considerably coarser grid than in the two-dimensional examples. We used 5 iterations (Section 2).For a three-dimensional system, we cannot plot the complete Lyapunov functionv as a function of its arguments, and hence we only plot the eigendirections of the Hessian on the underestimation Cv,0.7 of the chain recurrent set. In particular, we set γ+=0.7.Figure11 shows the connected components of the chain recurrent set Cv,0.7 close to z=−1 in the xy-plane Figure 11(a) and the xz-plane Figure 11(b). The 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. The 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) in z-direction Figure 11(d). The connected components of the chain recurrent set at z=1 show a similar behaviour, as shown in Figure12.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 two-dimensional 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).
(a)(b)(c)(d)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 two-dimensional 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).
(a)(b)(c)(d)Figure13 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)). The 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. The equilibrium 0,0,0 has a stable 3-dimensional manifold (green).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 one-dimensional 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.
(a)(b)(c)(d)Figure14 shows the overestimation X˜0=Cv,14 of the chain recurrent set (black) close to z=1. Note the much 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. The underestimation Cv,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.Figure 14
System (21). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. The green points represent the refined estimate X˜0 using the classification rule: two eigenvalues must be negative and one must be zero. (b) The results in the xy-plane, z=1.
(a)(b)


### 4.5. 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–10]. The 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. Thus, 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 Xeval.


## 4.1. One Periodic Orbit, One Equilibrium

The first system we considered was(18)x˙=−xx2+y2−1−y,y˙=−yx2+y2−1+x.This system has one unstable equilibrium at0,0 and an attractive periodic orbit at Ω=x,y∈ℝ2|x2+y2=1.We set the fineness parameterα=0.04 for the hexagonal collocation grid, and we used 10 iterations (Section 2). The computations were performed using f directly, rather than the normalised version in formula (4). The domain used was D=−1.2,1.22⊂ℝ2 , and to evaluate the function, we used the regular grid Xeval=hℤ∩D with parameter h=0.008.Figure1 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 Cv,0.8 of the chain recurrent set.Figure 1
System (18). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.8.
(a)(b)(c)(d)The 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. The Figure1(b) figure shows the orbital derivative: at most of the points, its value is close to −1 as a consequence of the collocation method. The points on the periodic orbit and around the critical point have values for the orbital derivative close to zero. The 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.Figure2 shows the eigenvalues and eigenvectors of the Hessian at the points of Cv,0.8. The upper figures show the eigenvalues, which are positive (green) and zero (black) on the periodic orbit. This indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). The 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.Figure 2
System (18). (a) Eigenvalues of the Hessian at points of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected component of the chain recurrent set. The periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, while the equilibrium has a 2-dimensional unstable manifold.
(a)(b)At the equilibrium point, we see two negative (red) eigenvalues, indicating a 2-dimensional unstable manifold. The 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 Figure3, we plot the overestimation X˜0=Cv,14 of the chain recurrent set Figure 3(a) and the underestimation Cv,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). The refined approximation is clearly superior to both the overestimated and the underestimated set.Figure 3
System (18). (a) The overestimation of the chain recurrent set X˜0=Cv,14. (b) The underestimation of the chain recurrent set Cv,0.8. (c) The refined approximation with our new method X˜0 of the chain recurrent set. (d) Comparison between the overestimated chain recurrent set (in black) and the reduced chain recurrent set (green) in a zoomed area.
(a)(b)(c)(d)Finally, in Figure4, 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.Figure 4
System (18). Eigenvalues of the Hessian at the points in X˜0 close to the origin. Since both are negative, v must have a maximum at the critical point 0,0.


## 4.2. Two Periodic Orbits, One Equilibrium

The second system we consider is(19)x˙=−xx2+y2−1x2+y2−14−y,y˙=−yx2+y2−1x2+y2−14+x.This system has a stable equilibrium at0,0, an attractive periodic orbit Ω1=x,y∈ℝ2|x2+y2=1, and a repelling periodic orbit Ω2=x,y∈ℝ2|x2+y2=1/4.For our computations, we used the normalisation (6), the domain D=−1.2,1.22⊂ℝ2, and the parameter values α=0.04 and h=0.008. We used 10 iterations (Section 2).Figure5 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 Cv,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.Figure 5
System (19). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.8.
(a)(b)(c)(d)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.Figure6 shows the eigenvalues and eigenvectors of the Hessian at points of the underestimation Cv,0.8 of the chain recurrent set. The upper figures show the eigenvalues: the eigenvalues are positive (green) and zero (black) on the outer periodic orbit. This indicates a 1-dimensional connected component of the chain recurrent set (one zero eigenvalues) that is attracting (the remaining one is positive). They 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).Figure 6
System (19). (a) Eigenvalues of the Hessian at points of the underestimation Cv,0.8 of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected components of the chain recurrent set. The outer periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, and the inner periodic orbit has a 1-dimensional neutral direction and a 1-dimensional unstable one, while the equilibrium has a 2-dimensional stable manifold.
(a)(b)The 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 Figure7(a), we plot the overestimation X˜0=Cv,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.Figure 7
System (19). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. Our refined approximation X˜0 of the chain recurrent set is plotted in green. It includes the underestimation Cv,0.8 and points from the overestimation X˜0 at the inner orbit if one eigenvalue is negative and one zero. At the outer orbit, it includes points from X˜0 if one eigenvalue is negative and one is zero. (b) Eigenvalues of the points in X˜0 close to the critical point at the origin. All are positive, so the equilibrium has a 2-dimensional stable manifold.
(a)(b)


## 4.3. Van der Pol

Next, we consider the classical van der Pol system(20)x˙=y,y˙=1−x2y−x.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,42⊂ℝ2, and the parameter values α=0.044, and h=0.008. We used 10 iterations (Section 2).Figure8 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 Cv,0.6 of the chain recurrent set. This 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.Figure 8
System (20). (a) Complete Lyapunov function v. (b) Its orbital derivative v′. (c) Norm of its gradient ∇v. (d) Underestimation of the chain recurrent set Cv,0.6.
(a)(b)(c)(d)Figure9 shows the eigenvalues and eigenvectors of the Hessian at points of the chain recurrent set. The top figures show the eigenvalues, which are positive (green) and zero (black) on the periodic orbit and both negative at the equilibrium (red). The lower figure shows the corresponding eigendirections: the eigenvectors corresponding to the zero eigenvalues are tangential to the periodic orbit, while the ones corresponding to the positive (attracting) one are perpendicular to it. At the equilibrium point, we see the corresponding repelling eigendirections.Figure 9
System (20) : (a) Eigenvalues of the Hessian at points of the chain recurrent set: green (positive, attractive), black (zero, neutral), and red (negative, repelling). (b) Corresponding eigendirections, indicating attractive, repelling, and neutral directions to the connected component of the chain recurrent set. The periodic orbit has a 1-dimensional neutral direction and a 1-dimensional attractive one, while the equilibrium has a 2-dimensional unstable manifold.
(a)(b)In Figure10(a), we plot the overestimation X˜0=Cv,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.Figure 10
System (20). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. Our refined approximation X0 of the chain recurrent set is plotted in green. We include points if there is one positive and one zero eigenvalues. (b) Eigenvalues of the points in X˜0 close to the critical point at the origin. Both are negative, so v must have a maximum, that is, the critical point is unstable with a 2-dimensional unstable manifold.
(a)(b)


## 4.4. Three-Dimensional Case

Finally, we consider the three-dimensional system(21)x˙=xx2+y2−1−xz2+1,y˙=yx2+y2−1+xz2+1,z˙=10zz2−1.This system has three equilibria at0,0,±1 and 0,0,0 as well as three periodic orbits Ωz0=x,y,z∈ℝ3|x2+y2=1,z=z0 for z0∈0,±1.For our computations, we used the normalisation (6), the domain D=−1.2,1.23⊂ℝ3, and the parameter values α=0.07 and h=0.015; note that this is a considerably coarser grid than in the two-dimensional examples. We used 5 iterations (Section 2).For a three-dimensional system, we cannot plot the complete Lyapunov functionv as a function of its arguments, and hence we only plot the eigendirections of the Hessian on the underestimation Cv,0.7 of the chain recurrent set. In particular, we set γ+=0.7.Figure11 shows the connected components of the chain recurrent set Cv,0.7 close to z=−1 in the xy-plane Figure 11(a) and the xz-plane Figure 11(b). The 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. The 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) in z-direction Figure 11(d). The connected components of the chain recurrent set at z=1 show a similar behaviour, as shown in Figure12.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 two-dimensional 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).
(a)(b)(c)(d)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 two-dimensional 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).
(a)(b)(c)(d)Figure13 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)). The 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. The equilibrium 0,0,0 has a stable 3-dimensional manifold (green).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 one-dimensional 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.
(a)(b)(c)(d)Figure14 shows the overestimation X˜0=Cv,14 of the chain recurrent set (black) close to z=1. Note the much 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. The underestimation Cv,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.Figure 14
System (21). (a) The overestimation of the chain recurrent set X˜0=Cv,14 is plotted in black. The green points represent the refined estimate X˜0 using the classification rule: two eigenvalues must be negative and one must be zero. (b) The results in the xy-plane, z=1.
(a)(b)


## 4.5. 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–10]. The 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. Thus, 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 Xeval.


## 5. 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. This enables the determination of the dimension of the components and the local directions of the stable and unstable manifolds. These 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].



---

*Source: 3160052-2022-08-08.xml*