COMPUTING CONTINUOUS AND PIECEWISE AFFINE LYAPUNOV FUNCTIONS FOR NONLINEAR SYSTEMS

. We present a numerical technique for the computation of a Lyapunov function for nonlinear systems with an asymptotically stable equilibrium point. The proposed approach constructs a partition of the state space, called a triangulation, and then computes values at the vertices of the triangulation using a Lyapunov function from a classical converse Lyapunov theorem due to Yoshizawa. A simple interpolation of the vertex values then yields a Continu- ous and Piecewise Aﬃne (CPA) function. Veriﬁcation that the obtained CPA function is a Lyapunov function is shown to be equivalent to veriﬁcation of several simple inequalities. Numerical examples are presented demonstrating diﬀerent aspects of the proposed method.


1.
Introduction. Lyapunov's Second or Direct Method [21] (see also [14,26,30]) has proved to be one of the most useful tools for demonstrating stability properties. This is largely due to the fact that if one has a Lyapunov function at hand there is no need to explicitly generate system solutions in order to determine stability. Unfortunately, this frequently trades the difficult problem of generating system solutions for the equally difficult problem of constructing a Lyapunov function.
Converse Lyapunov theorems provide existence results for Lyapunov functions; i.e., assuming a particular stability property holds then there exists an appropriate Lyapunov function [23,19,29,30,18]. However, such results are largely not constructive in nature and, in fact, depend explicitly on the solution trajectories of 228 SIGURDUR F. HAFSTEIN, CHRISTOPHER M. KELLETT AND HUIJUAN LI the system under study. As a consequence, various approaches have been proposed for the numerical construction of Lyapunov functions such as collocation methods [5,15], graph theoretic methods [2,16], semidefinite optimization for sum-of-squares polynomials (known as the SOS method) [24,25], and linear programming to generate continuous and piecewise affine (CPA) Lyapunov functions [22,1,9,12]. (See [11] for a survey of computational methods.) This latter approach, sometimes called the CPA method, is the starting point for this paper. In the CPA method, a domain of the state space is partitioned into simplices (called a triangulation) and a linear program is constructed to obtain numerical values at each simplex vertex. This linear program is constructed in such a way that the convex interpolation of these values yields a Lyapunov function that is CPA; that is, a CPA Lyapunov function. However, a shortcoming of this approach is that the linear program can be quite large with the number of variables being at least the number of vertices in the triangulation and the number of constraints being at least the number of simplices in the triangulation times the dimension of the state space. Consequently, solving the linear program can be quite slow.
In this paper we consider systems described by ordinary differential equationṡ where we assume f : R n → R n is twice continuously differentiable (i.e., f ∈ C 2 (R n , R n )), f (0) = 0, and denote the solution to (1) by φ : R ≥0 × R n → R n . It is well known that the asymptotic stability of the equilibrium at the origin is equivalent to the existence of a strict Lyapunov function for the system [18]. As an alternate approach to constructing a CPA Lyapunov function, we compute simplex vertex values by numerically approximating a Lyapunov function from the converse Lyapunov theorem demonstrated by Yoshizawa [29,30]. Verification that the resulting CPA function is in fact a CPA Lyapunov function can be done by checking straightforward linear inequalities similar to those that comprise the constraints in the linear programming approach. Several similar approaches for different kinds of systems and/or using different kind of numerical approximation of a Lyapunov function have surfaced in the literature recently, see e.g. [20,10,4,3]. We emphasize that the approach proposed in this paper only requires asymptotic stability of the equilibrium at the origin, not exponential stability. While the construction of Yoshizawa requires the solutions of (1) for every initial condition, only the solution over a finite time horizon is required. Furthermore, in our approach, this finite time solution is not required for every initial condition in the considered region, but only at the vertices of the triangulation. It is satisfaction of the aforementioned linear inequalities that is the crucial step in demonstrating a CPA Lyapunov function rather than constructing a numerical approximation of the construction of Yoshizawa. In practice, numerically approximating the construction of Yoshizawa provides a principled guess for the vertex values of the triangulation.
The benefit of this approach in constructing CPA Lyapunov functions over the linear programming approach is three-fold; (i) in all examples so far considered, a significant speed-up in computation time is achieved; (ii) it is usually possible to obtain a CPA Lyapunov function on a larger domain; and (iii) the algorithm always provides some information in the form of regions where the constructed CPA function is in fact a CPA Lyapunov function. With respect to this last point, when the linear programming method fails, it returns no information. While it is difficult to directly compare the computational burden of the linear programming approach and the approach proposed herein, in Section 5 both techniques are applied to four examples where our approach yields both shorter computation times and a larger domain of definition.
The paper is organized as follows: in Section 2 we describe the construction of CPA functions on a given triangulation and the linear inequalities used to verify if a given CPA function is, in fact, a Lyapunov function. In Section 3 we describe the Lyapunov function construction due to Yoshizawa and describe the form of the stability estimates required. In Section 4 we propose an algorithm for constructing CPA functions and verifying that they are CPA Lyapunov functions. In Section 5 we present several numerical examples and in Section 6 we provide some concluding remarks.
2. Continuous and piecewise affine Lyapunov functions. In the sequel, we will define continuous and piecewise affine (CPA) functions on suitable triangulations. For a set Ω ⊂ R n , we denote the interior of Ω by Ω • , the closure of Ω by Ω, the boundary of Ω by ∂Ω, and the complement of Ω by Ω C . For a vector x ∈ R n , we denote the 2-norm by |x|, the 1-norm by |x| 1 , and the infinity-norm by |x| ∞ . We denote the 2-norm of matrices by · . We denote the positive real numbers by R >0 and the nonnegative real numbers by R ≥0 . Given ε ∈ R >0 we define B ε := {x ∈ R n : |x| < ε}. We denote the closed convex hull of an ordered set of points x i ∈ R n , i = 0, 1, . . . , N by co(x 0 , x 1 , . . . , x N ). If the vectors x 0 , x 1 , . . . , x N ∈ R n are affinely independent, i.e. the vectors x 1 − x 0 , x 2 − x 0 , . . . , x N − x 0 are linearly independent, then the set S = co(x 0 , x 1 , . . . , x N ) is called an N -simplex. For any nonempty subset {x i0 , x i1 , . . . , x i K } of the vertices of S the K-simplex co(x i0 , x i1 , . . . , x i K ) is said to be a face of S. We make use of the standard function classes K, K ∞ , and KL (see [14,17] Definition 2.1. A finite collection T = {S 1 , S 2 , . . . , S N } of n-simplices in R n is called a suitable triangulation if i) S ν , S µ ∈ T intersect in a common face or not at all. ii) With D T := ∪ ν S ν , D • T is a simply connected neighborhood of the origin. iii) If 0 ∈ S ν , then 0 is a vertex of S ν . Remark 1. Property i), often called shape regularity in the theory of finite element methods, is needed so that we can parameterize every continuous function, affine on every simplex, by specifying its values at the vertices, cf. Remark 2. Property ii) ensures that D T is a natural domain for a Lyapunov function and, without Property iii), a function affine on each of the simplices could not have a local minimum at the origin.
In what follows, we will define simplices by fixing an ordered set of vertices and considering the closed convex hull of those vertices. For a given suitable triangulation, T , and with D T := ∪ S∈T S, we denote the set of all continuous functions f : D T → R that are affine on every simplex S ∈ T by CPA[T ].

Remark 2.
A function V ∈ CPA[T ] is uniquely determined by its values at the vertices of the simplices of T . To see this, let S ν = co(x 0 , x 1 , . . . , x n ) ∈ T . Every point x ∈ S ν can be written uniquely as a convex combination of its vertices, Additionally, V has a representation on S ν as V (x) = w T ν (x − x 0 ) + a ν for some w ν ∈ R n and some a ν ∈ R. In what follows, for V ∈ CPA[T ] and x ∈ S ν we denote Then, as shown in [9,Remark 9], ∇V ν is linear in the values of V at the vertices For a locally Lipschitz function V : R n → R ≥0 , the upper Dini derivate at x ∈ R n in the direction w ∈ R n is defined by Our subsequent results will be valid on a domain D ⊂ R n minus a fixed arbitrarily small neighborhood of the origin. We define a CPA[T ] Lyapunov function that accounts for this. Definition 2.2. Let T be a suitable triangulation and let V ∈ CPA[T ] be a positive definite function. Let ε ∈ R >0 be such that If there is a constant α * ∈ R >0 such that The implication of a CPA[T ] Lyapunov function for (1) on D T \ B ε is slightly weaker than asymptotic stability of B ε as we make precise in the following theorem. This is frequently referred to as "practical stability" in the literature or as "ultimate boundedness" in [14]. By a slight abuse of notation, for a set Ω ⊂ R n , we denote the reachable set of (1) from Ω at time t ∈ R ≥0 by φ(t, Ω) := ∪ x∈Ω φ(t, x).  . Let S ν = co(x ν 0 , x ν 1 , . . . , x ν n ) ∈ T and let µ ν ∈ R ≥0 satisfy max i,j,k=1,2,...,n x∈Sν For each S ν , for i = 0, 1, . . . , n define the constants Then, for every S ν such that the inequalities The usefulness of Theorem 2.4 is that it reduces the verification that a given function V ∈ CPA[T ] is a Lyapunov function for (1) to the verification of a finite number of inequalities (7). In the linear programming approach used in [1,9,12,22], the linear inequalities are used as constraints in a linear program and, hence, a solution necessarily satisfies (7). By contrast, in this paper, we propose fixing the vertex values by a computational procedure described in the next section, and subsequently verifying the inequalities (7).
We now turn to the question of the existence of a CPA[T ] Lyapunov function. As we will demonstrate in Theorem 2.7, if a CPA[T ] function approximates a twice continuously differentiable Lyapunov function, and if the triangulation T is sufficiently fine, then the CPA[T ] function is in fact a CPA[T ] Lyapunov function. To do this, we require the following definitions. We additionally need that the simplices in the triangulation T are not too close to being degenerate; that is, no n-simplex should be close to being of dimension n − 1. This property can be quantified as follows: For an n-simplex S ν := co(x 0 , x 1 , . . . , x n ) ∈ T define its shape-matrix as X ν by writing the vectors Note that, as we have defined simplices based on an ordered set of vertices, the shape matrix thus defined is unique. The degeneracy of the simplex S ν is quantified by is the spectral norm of the inverse of X ν (see part (ii) in the proof of [1,Theorem 4.6]). To see why this quantity captures a "distance-to-degeneracy" of the n-simplex S ν , observe that degeneracy corresponds to the presence of linearly dependent rows in X ν , resulting in X ν being singular. If rows are nearly linearly dependent, possibly as a result of vertices being close to each other, then the spectral norm of X −1 ν will be large. Of course, we may wish to use very small simplices in order to reduce the error between a given Lyapunov function and its CPA approximation, and hence a reasonable measure of distance-todegeneracy should also scale the spectral norm of the inverse of X ν by the diameter of the simplex, leading to the quantity diam(S ν ) X −1 ν . Definition 2.6. Given a neighborhood of the origin D ⊂ R n , a locally Lipschitz Note that, as the above definition is essentially local in nature (unless the neighborhood D is the entire space), the functions α, α 1 ∈ K can be replaced by positive definite functions α, α 1 : R ≥0 → R ≥0 which, locally, can always be bounded from below by functions of class-K.
While the above definition with W locally Lipschitz is sufficient to conclude local asymptotic stability of the origin, the following result requires a twice continuously differentiable Lyapunov function in order to obtain certain numerical estimates in the proof.
Then for every R ∈ R >0 there exists a δ R ∈ (0, ε) so that, for any suitable triangulation T satisfying Further, for any large enough R ∈ R >0 there are such suitable triangulations.
Proof. For any R ∈ R >0 so large that (14) can be fulfilled one can construct a suitable triangulation satisfying (12), (13), and (14). For example, it is not difficult to see that with R = 2 one can construct such a suitable triangulation satisfying (12), (13), and (14) as in [9,Definition 13]. Indeed, one can take any δ R between zero and ε that is smaller than inf{|x − y| : x ∈ C, y ∈ D C } and the triangulation T C K,b in [9, Definition 13] with K = 0 and b = δ R / √ n. In summary, this triangulation starts from integer grid points that are then scaled down by the constant b. Simplices that do not intersect the interior of C are then discarded. For the rest of the proof assume that we have such a triangulation T and that δ R is so small that (12) and (13) are fulfilled. We first derive some inequalities and then we fix δ R > 0 so small that the CPA[T ] approximation V to W on D T is a CPA Lyapunov function for (1) on For an arbitrary but fixed S ν = co( and define A := max z∈D, i,j=1,2,...,n Let X ν be as in (8) and define χ := max ν X −1 ν . Following the proof of part (iii) of [1, Theorem 4.6] we can show that Define and observe that since and, using the definitions (15) and (8) Since W ∈ C 2 (D), ∇W (x) is bounded on the compact set D and we can define holds uniformly in ν. Let ∇V ν,i denote the i th component of ∇V ν . We then see that |∇V ν,i | ≤ G and hence |∇V ν | 1 ≤ nG. Define and let E i,ν ∈ R ≥0 be defined by (6) with µ ν = µ * . Then, from (13) and (6) we have Using (10), (16), and (17) we calculate Now, fix δ R ∈ (0, ε) so that 2δ R nA 1 2 χn Since δ R < and S ν ∩ B C = ∅ we have |x i | ≥ − δ R > 0. Since α ∈ K, and with the bounds (21) and (20), the linear constraints are satisfied for all vertices x i of S ν . Further, because V is defined as interpolated values of W , we have by (11) Since W is positive definite, so is V . Consequently, Corollary 1 proves the theorem.
Theorem 2.7 implies that it is always possible to find a triangulation that admits a CPA Lyapunov function approximating a twice continuously differentiable Lyapunov function.
We note that the assumption of twice differentiability is required in proving the bound (16), which, with (18), can be seen to be a bound on the difference between the slope of the CPA approximation on S ν and the gradient of the Lyapunov function W at each vertex of S ν ; i.e., a bound on |∇V ν − ∇W (x i )|. Since the right-hand side of (16) goes to zero as the diameter of the simplex S ν goes to zero, we see that ∇V ν being close to ∇W (x i ) for all vertices defining the simplex requires at least continuity of ∇W (x). In fact, as can be seen from the definition of the constant A, what is additionally required is that the second derivative of W needs to exist and be bounded inside each simplex.
3. Yoshizawa construction of Lyapunov functions. We now turn to the question of how to define the vertex values of each simplex in order to obtain a CPA Lyapunov function. We propose using a numerical approximation of a construction initially proposed by Yoshizawa in proving a converse Lyapunov theorem [29, Theorem 1].
Let the open set D ⊂ R n be such that D is forward invariant for (1) and the origin is contained in D. Suppose (1) is KL-stable on D; i.e., there exists β ∈ KL so that |φ(t, x)| ≤ β(|x|, t), ∀x ∈ D, t ∈ R ≥0 . (23) It was shown in [28, Proposition 1] that KL-stability is equivalent to (local) asymptotic stability of the origin for (1) where D is contained in the basin of attraction. See also [14,Definition 2.9] where asymptotic stability is defined in terms of a bound of class-KL. When D = R n , KL-stability is equivalent to global asymptotic stability of the origin for (1). We will refer to the function β ∈ KL of (23) as a stability estimate.
Definition 3.2. Given a stability estimate β ∈ KL, let α 1 , α 2 ∈ K ∞ come from Lemma 3.1 with λ = 2. We call the function V : R n → R ≥0 defined by a Yoshizawa function for (1).  (24) is continuous on D and locally Lipschitz on D\{0} and satisfies and the decrease condition for all x ∈ D and all t ∈ R ≥0 . Furthermore, with T : D\{0} → R ≥0 defined by for all x ∈ D\{0}, we have Observe that, for any x ∈ D\{0}, taking the maximum over any interval [0, T ] where T ≥ T (x) will not change the value of the Yoshizawa function. Furthermore, since the Yoshizawa function is locally Lipschitz, (25) and (26) imply that the Yoshizawa function is a Lyapunov function when the system under study is KLstable. Furthermore, the argument for (26) follows that of [28, Section 5.
Similarly, that the Yoshizawa function is locally Lipschitz on D\{0} is shown in [28, Section 5.1.2]. In [28,Claim 2] it is shown that for T : D\{0} → R ≥0 given by the Yoshizawa function satisfies By using the upper and lower bounds (25) we see that  (7) at each vertex of the triangulation. 5. If necessary, refine the triangulation and repeat steps 2-4. Computationally, steps 1 and 5 are also a feature of the linear programming approach [22] to computing CPA Lyapunov functions, though it may be necessary to refine the triangulation in order for the linear program to be feasible. By contrast, the calculations proposed in Algorithm 1 can be carried out for any triangulation. Assuming a triangulation that admits a feasible solution to the linear program of [22], the difference between the linear programming approach [22] and the approach proposed in Algorithm 1 lies in steps 2-4. Steps 3 and 4 are computationally straightforward.
Step 2 requires some discussion.
In computing the Yoshizawa function (28), we require a stability estimate β ∈ KL, functions α 1 , α 2 ∈ K ∞ from Lemma 3.1, and, at each vertex of the triangulation, a solution to (1) over the finite time window [0, T (x)]. We first address issues with the solution and finite time window and then comment on the stability estimate and K ∞ functions.
As a closed form solution of (1) is generally not available, we will resort to numerical integration in order to obtain an approximate solution φ(t, x) for use in the calculation of V (x) given by (28). For this approach to be numerically tractable, it is important that the time horizon T (x) given by (27) is not too large.
For an exponential stability estimate β ∈ KL bounded as β(s, t) ≤ α(s)e −µt for some µ ∈ R >0 , α ∈ K ∞ satisfying α(s) ≥ s, and all (s, t) ∈ R 2 ≥0 , the functions Note that the usual definition of exponential stability via a KL-estimate requires a KL function of the form β(s, t) = kse −λt for k > 1 and λ > 0. Hence, while (34) decays exponentially fast in time, the transient overshoot also grows exponentially, rather than linearly, with the size of the initial condition.
Remark 4 (Stability Estimates). There are two difficulties we encounter in trying to calculate (24). The first difficulty lies with finding a stability estimate β ∈ KL or even with verifying that a particular stability estimate such as (34) holds for a particular system (1). There seems to be little that can be done to circumvent this problem. However, in practice, since we only compute the Yoshizawa function on a compact domain containing the origin, a global stability estimate is not required.
The second difficulty is that Sontag's lemma on KL-estimates is not constructive and, to the best of the authors' knowledge, given an arbitrary β ∈ KL, there are currently no constructive techniques for finding α 1 , α 2 ∈ K ∞ . In what follows we used a multi-threaded four-step Adams-Bashforth solver to perform the numerical integration. In the first three examples we explicitly derive stability estimates for the Yoshizawa construction and use these parameters in the computation of the Lyapunov functions. In the fourth example, which is much more demanding, we omit this and found after some experimenting that α 1 (s) = s 2 and T = 3 deliver a solution. Similar CPA Lyapunov functions were constructed using a simple Euler integration scheme. The results thus obtained are similar to those presented and are thus omitted. 5.1. Triangulation. The triangulation in the subsequent examples involves a fanlike region around the origin. Such a triangulation was first proposed in [6] for two dimensional systems and generalized in [8] to arbitrary dimensions. As a complete enumeration of edges in the triangulation is quite involved, we refer to [6,8] for a detailed description of the construction of such triangulations. Here, we provide an abbreviated description that captures the essence of the construction in two dimensions.

SIGURDUR F. HAFSTEIN, CHRISTOPHER M. KELLETT AND HUIJUAN LI
Fix positive integers K, k ∈ Z >0 , with K > k, and define the preliminary vertex set by (35) Simplex edges are then defined as shown in Figure 1 for k = 2 and K = 4. x 2 x 1 Figure 1. Initial triangulation before scaling.
To obtain the triangulation used for the examples that follow, this initial triangulation of Figure 1 is then scaled by the mapping F : R n → R n defined by For the case where K = 4 and k = 2, this yields the triangulation shown in Figure 2. This scaling results in vertices lying on concentric circles or spheres. Other scalings can certainly be used. Note that in defining the initial vertex set we fixed a large square, characterized by K, and excised a smaller square, characterized by k, from the interior. However, for systems with complicated dynamics or for equilibria with complicated regions of attraction, particularly in higher dimensions, it can be useful to define more complicated regions. Using the two-dimensional case as an example, one straightforward modification to achieve this is to allow different constants defining the initial vertex region; i.e., choose K i , k i ∈ Z >0 , i = 1, 2, 3, 4, and then define the preliminary vertex set using the regions [−K 1 , Computation times for each of the three following examples using this triangulation and the Adams-Bashforth integration scheme are summarized in Table 1 for both the approach proposed herein and the linear programming approach of [22]. By way of comparison, an alternate triangulation is used for Example 4 while still using the Adams-Bashforth solver.

5.2.
Example 1 -First order nonlinear system. Consider the systeṁ which has solution The origin is thus asymptotically stable, but not exponentially stable. We observe that the norm of the solution is in fact a KL function and, consequently, immediately provides a stability estimate. We can verify that the functions With α 1 = α 2 , we see that, from (27), T (x) = T = 1.
We define a triangulation with K = 200 and k = 19, where modifying the procedure of Section 5.1 for a one-dimensional system is straightforward. We calculate the values at the simplex vertices by approximating (28) and an affine interpolation of these values on each simplex then yields a CPA function. We verify that the inequalities (7)  We note that, for any p ∈ Z ≥1 and c ∈ R >0 , a Lyapunov function for (37) is given by (38) Figure 3 shows the CPA-Lyapunov function V 2 (x) for system (37) as well as the known Lyapunov function (38) with p = 2, c = 0.01.
We observe that the origin is globally exponentially stable as the eigenvalues are at −1 ± i. By solving the Lyapunov equation A T P + P A = −Id, a Lyapunov function is given by By explicitly calculating the solutions of (39) we see that the system satisfies the stability estimate From (32), with α(s) = 7s and µ = 1 we see that α 1 (s) = s 2 , α 2 (s) = 49s 2 , and T (x) = T = 4.892 computed via (27). According to the above proposed procedure, we define a triangulation as described in Section 5.1 with K = 90 and k = 20 and the scaling F of (36). The values at the simplex vertices are given by approximating the solution of (39) by numerical integration over the time window [0, 4.892] and computing the value of the Yoshizawa function (28). This then defines a continuous and piecewise affine function V 1 (x) for (39) as shown in Figure 4. It is straightforward to numerically verify that the inequalities (7) are satisfied for all simplex vertices where S ν ∩ B C 0.048 = ∅ and hence V 1 (x), is a CPA-Lyapunov function on D T \B 0.048 = B 0.972 \B 0.048 .
The function V (x) given by (40) has a similar though slightly different shape. Level curves for V 1 (x) are shown in Figure 5 and level curves for V (x) are shown in Figure 6 for comparison. While V 1 is a CPA approximation of the true Yoshizawa function, the corners in Figure 5 indicate that the Yoshizawa function is not continuously differentiable. The implicit function theorem delivers the existence of a continuously differentiable function, if the maximum occurs at a unique time t 1 . The nondifferentiability appears to be due to the fact that the time at which the maximum occurs in (24) is not a continuous function of x and, in fact, on either side of the corners (i.e., approximately either side of the line x 1 = 0) the maximum is attained at two different times t 1 , t 2 ∈ [0, T ], t 1 = t 2 .

5.4.
Example 3 -Second order nonlinear system. Consider the two-dimensional nonlinear system given bẏ This system has the unit circle as a periodic orbit and the origin as a locally exponentially stable equilibrium. On any compact subset of the unit ball, the simple quadratic is a known Lyapunov function. Fix R ∈ (0, 1). Then, for any initial conditions satisfying and, from (32), we can calculate and, from (27), T (x) = T = 1. Note that (44) indicates that the origin is a locally exponentially stable equilibrium, though it is clearly not exponentially stable in the large on its domain of attraction.
For this example,with R = 0.94478 and using the numerical procedure previously outlined, a CPA-Lyapunov function V 3 (x) of system (42) was computed and is shown in Figure 7. The triangulation is defined as in Section 5.1 with K = 90, k = 10, and the mapping F of (36), yielding a region that is coincident with the region on which the stability estimate is valid; i.e., on B √ R . The inequalities (7) hold for all simplices such that S ν ∩ B 0.012 = ∅ and so  The extreme shape of V 3 (x) as shown in Figure 7 demonstrates the exponential decrease obtained by the Yoshizawa function. In particular, near the periodic orbit, but in the domain of attraction of the origin, the convergence rate is very slow and hence the Lyapunov function must be very steep. As the trajectory nears the origin, its convergence rate increases and hence the Lyapunov function need not be as steep. Furthermore, the stability estimate (44) becomes increasingly conservative as trajectories move away from the boundary of the domain B √ R , contributing further to the very flat Lyapunov function near the origin. 5.5. Example 4 -Third order nonlinear system. Finally, consider the third order systemẋ The Jacobian at the origin has eigenvalues −1 ± i √ 2 and −3 and the origin is thus a locally exponentially stable equilibrium.
Fix the scaling in the Yoshizawa function (28) Figure 8. Note that the level sets are not, in fact, spheres and the level set shown in Figure 8 is squashed or flattened in the region x 2 , x 3 > 0 and x 1 < 0. By way of comparison, we also applied the linear programming approach proposed in [22], with the largest level set obtained shown in Figure 9. The triangulation used in this computation is given by 0.01(±i 2 , ±j 2 , ±k 2 ) for i, j, k = 0, 1, . . . , 9. Note that the obtained level set is on a domain with a radius less than half the size of that obtained via the proposed approach. Despite the fact that the triangulation used in Figure 9 used 44 times fewer grid points (i.e., 6,859), the computation took almost 60 minutes on the same PC using the state-of-the-art Gurobi Linear Program solver.
In order for the linear program to have a solution, there cannot be any constraint violations anywhere on the computational domain. In particular, for the considered example it was necessary to define a "quadratic" triangulation, so that one has smaller simplices closer to the origin, in order to obtain a feasible solution. However, using Algorithm 1, one can define a domain as large as desired and then, by checking the linear inequalities (7) at each vertex, determine a region where the orbital derivative is negative. This significantly simplifies the initial setup of the computational problem as choosing a computational domain larger than the basin of attraction does not lead to an infeasible problem necessitating a refinement of the computational domain. 5.6. Computation times. Table 1 summarizes the computation times required for the four presented examples. The verification of the inequalities (7) (denoted by 'VT' in Table 1) requires negligible time in comparison to the computation of the values of the Lyapunov function (denoted by 'CT' in Table 1). Also note that Figure 9. From the LP approach of [22], the largest level set containing the origin. relatively more time is needed to construct the simplicial complex (denoted 'TC' in Table 1) and for solution of the LP problem (denoted 'LP' in Table 1). This can be considerably simplified if the intention is only to use the Yoshizawa construction. For the triangulations referenced in the table, using the Yoshizawa construction leads to the necessary inequalities (7) being satisfied at every vertex; i.e., for the described triangulation, we successfully computed a CPA Lyapunov function on that triangulation. i4790K@4.6GHz with 32GB memory). 'TC' denotes the time needed to create the triangulation and set up the LP problem, 'CT' the time required to compute the Yoshizawa function at all vertices, 'VT' the time required for verification of the linear inequalities (7), and 'LP' the time needed to solve the LP problem with Gurobi using default parameters.
-No solution after 4 hours. † -Not appropriate here, c.f. text to Example 4.

6.
Conclusions. In this paper we have presented a novel technique, summarized in Algorithm 1, for the numerical construction of Lyapunov functions given a stability estimate in the form of a KL-bound on the norm of system trajectories. For a suitable triangulation of the state space, at each simplex vertex we calculate the value of a Lyapunov function construction due to Yoshizawa [29,30]. From these values, we then define a CPA function on the domain minus an arbitrarily small neighborhood of the origin. We can verify that the CPA function thus defined is a Lyapunov function (Corollary 1) by checking a simple linear inequality (7) at each vertex of the triangulation.
It is important to note that any CPA function that satisfies the inequalities (7) is a CPA Lyapunov function. Theorem 2.7 guarantees that a CPA function that approximates a twice continuously differentiable Lyapunov function is, in fact, a CPA Lyapunov function. In this sense, the method proposed in Algorithm 1 can be seen as a way to make a "principled guess" for a CPA function that is likely to satisfy (7), despite possible crude approximations made in the process of computing the Yoshizawa function.
We observe that in the numerical examples of Section 5 there is a significant improvement in both decreasing computation time and increasing computation domain when using Algorithm 1 over the linear programming approach of [22]. Further reductions in computation time can be made by moving to a parallel computation architecture based on the observation that Steps 2 and Steps 4 of Algorithm 1 can be done for each vertex independent of every other vertex.