GRID REFINEMENT IN THE CONSTRUCTION OF LYAPUNOV FUNCTIONS USING RADIAL BASIS FUNCTIONS

Lyapunov functions are a main tool to determine the domain of attraction of equilibria in dynamical systems. Recently, several methods have been presented to construct a Lyapunov function for a given system. In this paper, we improve the construction method for Lyapunov functions using Radial Basis Functions. We combine this method with a new grid refinement algorithm based on Voronoi diagrams. Starting with a coarse grid and applying the refinement algorithm, we thus manage to reduce the number of data points needed to construct Lyapunov functions. Finally, we give numerical examples to illustrate our algorithms.


1.
Introduction. The determination of the domain of attraction of an equilibrium is an important task in the analysis and derivation of dynamical systems, arising in many practical applications. Since it is difficult to determine the exact domain of attraction analytically, researchers have been seeking numerical algorithms to determine subsets of the domain of attraction. Most of these methods for computing domains of attraction are based on Lyapunov functions, which are functions that decrease along trajectories of the dynamical system. Sublevel sets of Lyapunov functions are positively invariant subsets of the domain of attraction. The construction of such Lyapunov functions, however, is very challenging.
In the last decades, several numerical methods to construct Lyapunov functions have been developed, for a review see [9]. These methods include the SOS (sums of squares) method, which is applicable for polynomial vector fields, introduced in [23] and available as a MATLAB tool box [22]. It constructs a polynomial Lyapunov function by semidefinite optimization.
The CPA method constructs a CPA (continuous piecewise affine) Lyapunov function using linear optimization [13,14]. A simplicial complex is fixed and the space of CPA functions which are affine on each simplex is considered. This space can be parameterized by the values on the vertices. The conditions of a Lyapunov function are transformed into a set of finitely many linear inequalities at the vertices, which include error estimates ensuring that the CPA Lyapunov function has negative orbital derivative inside each simplex. These linear inequalities are used as constraints of a linear programming problem, which can be solved by standard methods. While in the original method an arbitrarily small neighborhood of the equilibrium had to be cut out, a revised method can construct a CPA Lyapunov function also near the equilibrium by using a fan-like triangulation near the equilibrium [8].
A different method deals with Zubov's equation and computes a solution of this partial differential equation (PDE) [3]; the corresponding generalized Zubov equation is a Hamilton-Jacobi-Bellmann equation. This equation has a viscosity solution which can be approximated using standard techniques after regularisation at the equilibrium, for example one can use piecewise affine approximating functions and adaptive grid techniques [11].
The cell mapping approach [15] or set oriented methods [4] divide the phase space into cells and compute the dynamics between these cells; they have also been used to construct Lyapunov functions [12].
The RBF (Radial Basis Function) method, a special case of mesh-free collocation, considers a particular Lyapunov function, satisfying a linear PDE and solves it using mesh-free collocation [6,10]. For this method, a set of scattered grid points is used to find an approximation to the solution of the linear PDE. It is computed by solving a linear system of equations, for more details see Section 3.
In this paper, we extend the Radial Basis Function method to construct Lyapunov functions which will in turn determine subsets of the domain of attraction through sublevel sets of these Lyapunov functions. Our new contribution is to combine this construction method with a grid refinement algorithm in order to reduce effort and computation time.
The RBF method is very well suited for a refinement algorithm as the grid points can be scattered and do not require a special structure. Usually, if too few grid points are used, the constructed function has some areas with positive orbital derivative. Instead of refining the grid uniformly, the idea of the refinement algorithm is to refine the grid only where we need to do so, i.e., in areas with positive orbital derivatives. This is particularly useful when going to higher dimensions to save computation time both when constructing the function and evaluating the function, e.g. to determine sublevel sets.
The outline of the paper will be as follows: Section 2 gives an introduction to dynamical systems and Lyapunov functions. In Section 3 we explain the steps of the RBF construction method in detail. In Section 4 we present our refinement algorithm and apply it to numerical examples. Finally, Section 5 concludes the paper.
2. Dynamical systems and Lyapunov functions. Consider the following autonomous system of differential equations, which describes a dynamical systeṁ We denote by S t ξ := x(t) the solution of the initial value problemẋ = f (x), x(0) = ξ ∈ R d , t ≥ 0 and we assume that the solution exists for all positive times. Note that S t defines a dynamical system on Then, x(t) = x 0 is a constant solution of (1) for all t ≥ 0. The stability of an equilibrium point is determined by the behaviour of solutions in a neighborhood of an equilibrium. An equilibrium is called stable, if all solutions starting near the equilibrium point stay near the equilibrium for all future times. Moreover, it is called asymptotically stable if it is stable and the solutions starting near the equilibrium converge to it as time tends to infinity. It is called exponentially stable, if the convergence to the equilibrium is exponential. For an asymptotically stable equilibrium x 0 we are interested in finding the largest set of initial states from which the trajectories of solutions converge to the equilibrium as time tends to infinity. This set is called the domain of attraction.
Definition 1 (Domain of attraction). The domain of attraction of an asymptotically stable equilibrium x 0 is defined by Remark 1. The domain of attraction A(x 0 ) of an asymptotically stable equilibrium x 0 is non-empty and open.
2.1. Lyapunov Functions. The method of Lyapunov functions enables us to determine subsets of the domain of attraction of an asymptotically stable equilibrium through sublevel sets. A function V ∈ C 1 (R d , R) is called a Lyapunov function for the equilibrium x 0 if it has a local minimum at x 0 and a negative orbital derivative in a neighborhood of x 0 .
where ·, · denotes the standard scalar product in R d .

Remark 2.
The orbital derivative is the derivative along solutions: with the chain rule we have The following theorem shows how Lyapunov functions are used to find subsets of the domain of attraction; note that the requirement of the local minimum at x 0 is a consequence of the assumptions. The theorem states that sublevel sets of a Lyapunov function are positively invariant subsets of the domain of attraction, see e.g. [6,Theorem 2.24].
, V is decreasing along solutions in K \{x 0 }. Then K ⊂ A(x 0 ), K is positively invariant and V is called a Lyapunov function.

2.2.
Existence of Lyapunov functions. There are several results on converse theorems, ensuring the existence of Lyapunov functions in different contexts, e.g. already in 1949 [21]; for a review see [19]. These converse theorems, however, use the explicit solution of (1) and are thus often not useful to construct a Lyapunov function explicitly.
We will assume that x 0 is an exponentially stable equilibrium and we will consider two classes of Lyapunov functions such that V (x) < 0 holds for all x ∈ A(x 0 ) \ {x 0 }. We will characterise them by equations for V (x); these are linear first-order Partial Differential Equations (PDE) for V . Both functions have the same degree of smoothness as f , i.e. they are C σ functions. In this paper, we will extend the RBF method to construct Lyapunov functions by a refinement algorithm. In the next section we will recall the main ingredients of this method before we introduce the refinement in Section 4.
3. Construction of Lyapunov functions using Radial Basis Functions. The two Lyapunov functions T and Q, introduced in the previous section, satisfy partial differential equations for their orbital derivative. Therefore, their solutions can be approximated using the Radial Basis Functions collocation method. Note that the approximation itself will be a Lyapunov function.
3.1. Numerical solutions for PDE's using the Radial Basis Functions collocation method. Meshless collocation based on Radial Basis Function is an effective tool to solve linear PDE's. It has outstanding properties, such as approximating arbitrarily scattered data in multidimensional space as well as providing high order of accuracy which have made it a preferable method for the numerical solutions of partial differential equations. For a general introduction to Meshless collocation, in particular Radial Basis Functions, see [2,28]. For the application of RBF to the construction of Lyapunov functions, see [6], where details for the following overview of the method can be found, as well as [10].
A Radial Basis Function is a real-valued function whose value depends only on the distance from the origin i.e., Ψ(x) = ψ( x ), or alternatively from another point called centre, i.e., Ψ(x, x j ) = ψ( x − x j ), where · denotes the Euclidean norm in R d . There is a one-to-one correspondence between a Radial Basis Function, or more generally kernel, and its Reproducing Kernel Hilbert Space (RKHS). The approximate solution of the PDE will be a norm-minimal interpolant in the RKHS; in our brief overview, however, we do not discuss this relation further, the interested reader is referred to [10].
In the following section we will introduce a family of compactly supported Radial Basis Functions that enables us to approximate functions with certain smoothness i.e., not necessarily C ∞ . Note that the corresponding RKHS is a Sobolev space with equivalent norm.
The Wendland functions, introduced by Wendland [27], are compactly supported Radial Basis Functions, which are polynomials on their support. for r ∈ R + 0 . Here we set x + = x for x ≥ 0 and x + = 0 for x < 0. If we fix the parameter l := d 2 + k + 1 depending on the space dimension d and the parameter k, then the function Ψ(x) = ψ l,k (c x ) with c > 0 is a C 2k function with compact support. For dimensions d = 2 or d = 3, we give some Wendland functions in the following table.
k ψ l,k 1 ψ 3,1 (cr) = (1 − cr) 4 + (4cr + 1) 2 ψ 4,2 (cr) = (1 − cr) 6 + (35c 2 r 2 + 18cr + 3) 3 ψ 5,3 (cr) = (1 − cr) 8 + (32c 3 r 3 + 25c 2 r 2 + 8cr + 1) 4 ψ 6,4 (cr) = (1 − cr) 10 + (429c 4 r 4 + 450c 3 r 3 + 210c 2 r 2 + 50cr + 5) Table 1. The Wendland functions ψ 3,1 (cr), ψ 4,2 (cr), ψ 5,3 (cr) and ψ 6,4 (cr) with l = k + 2 and a scaling parameter c > 0. Note that these functions are the Wendland functions of Definition 3 up to a constant. Now let us consider a general linear partial differential equation of the form where L is a linear differential operator of the form In our case, the differential operator L will be given by the orbital derivative of a function V with respect to system (1), namely The operator L in (8) is a first order differential operator of the form (7) with where the superscript x denotes the application of the operator with respect to the variable x. The approximant v : R d → R of V will be given by where Ψ(x) is the Radial Basis Function. The coefficients β k are determined by claiming that the interpolation condition is satisfied for all grid points x j ∈ X N , or in other words that the PDE is satisfied at all points x j ∈ X N . This will lead to a linear system for β, i.e., Aβ = α. If the points x j are pairwise distinct and no equilibria, the symmetric matrix A is positive definite, so in particular non-singular. Hence, the system has a unique solution β.
The interpolation matrix entries of A = (a jk ) j,k=1,...,N are given by and the right-hand side α = (α j ) j=1,...,N is given by which are chosen to be one of our Lyapunov functions V ( Finally, we calculate the approximant v(x) and its orbital derivative v (x), using the following formulas, by evaluating and taking the orbital derivative of (9).
Remark 3 (The value ofc). Changing the value ofc has the effect of multiplying the solution β of the linear system Aβ = α by a positive constant. As a consequence also the value of the approximant v and its orbital derivative v will be multiplied by the same positive constant. This means that the areas of the phase space, where v is positive, are independent of the value ofc. Indeed, this follows from Aβ = α = −c = −(1, 1, . . . , 1) Tc , since the interpolation matrix A is independent of the value ofc and from formulas (10) and (11).
The following error estimate was given in [10,Corollay 4.11]. Note, that W τ 2 (Ω) denotes the usual Sobolev space on Ω ⊂ R d .
The reconstruction v of the Lyapunov function V with respect to the operator (8) and a set K ⊂ Ω : GRID REFINEMENT IN THE CONSTRUCTION OF LYAPUNOV FUNCTIONS   7 where h := sup x∈Ω min xj ∈X N x − x j denotes the fill distance, i.e. the maximum distance which a point in Ω can have from any point in X N .
Remark 4. The set K in the theorem above can be any compact subset of A(x 0 ) as r can be chosen so large that the sublevel set Ω of V satisfies Ω ⊃ K. A similar statement for the function satisfying V (x) = −c follows also from [10, Corollay 4.11], but with no data sites on the boundary. Note, however, that in this case the function satisfying V (x) = −c is not unique up to a constant, the error estimates on the orbital derivative, however, still hold.
The error estimate (14) implies that the approximant v of the Lyapunov function V is actually a Lyapunov function, i.e. satisfies v (x) < 0, if the grid is dense enough. Let us make this more precise: by choosing the fill distance h so small that Remark 5. In both cases, the approximation may fail to have negative orbital derivative near the equilibrium x 0 ; in the first case the error estimate requires x − x 0 2 > , and in the second case the function V is not defined in x 0 . If the equilibrium is exponentially stable, one can use the Lyapunov function of the linearised system, the so-called local Lyapunov function, to deal with this small neighborhood of x 0 , for details see [6], or a modified method, see [7]. In this paper, we will not deal with this local problem in more detail.
The error estimate uses the fill distance as a measure; hence, we are led to choose a uniformly fine grid. In examples, however, it turns out that an approximation with negative orbital derivative can be achieved with fewer points using a nonuniform grid. Moreover, the goal is not necessarily to have a good approximation of V , but to construct a function with negative orbital derivative. For example, when approximating the solution of V (x) = − x − x 0 2 , a larger error is permissable for points far away from the equilibrium.

3.2.
Steps of the construction method. The construction method is based on considering the Lyapunov function satisfying the PDE's stated in Section 2.2. We will explain this construction method in an example.

Example 1. Consider the simple linear system
The system has one asymptotically stable equilibrium x 0 = (0, 0). Now we will go through the steps of the method.
Here, we fix a regular grid where h > 0 is the distance between points in x-and y-direction and will be specified below.

NAJLA MOHAMMED AND PETER GIESL
• Use the RBF method to approximate the Lyapunov function V = T which satisfies T (x) = −1 by the approximant v and then calculate its orbital derivatives v using (10) and (11). • From the ansatz of the Lyapunov function V , we know that the orbital derivative is negative at every point in our grid, i.e., v (x i ) < 0 for all x i ∈ X N ⊂ R 2 , but there may be points in [−1, 1] 2 , where the orbital derivative is positive. The error estimate tells us that if h is small enough, the orbital derivative will be negative except for a small neighborhood of the equilibrium. Hence, we start with a certain h, and check the sign of the orbital derivative at the points between the grid points. If we have points where v (x) > 0, we need a finer grid. Therefore, we choose h smaller and smaller until we have v (x) < 0 for all points in the desired area. In this example, we used h = 1/6, resulting in  During the calculations of the Lyapunov function v we may find some points with positive orbital derivatives. This occurs because either the grid we have chosen is not fine enough, or the grid points do not all lie in the domain of attraction. We cannot exclude the second case, so we use a finer grid to calculate the Lyapunov function v and see if the problem will be solved. However, instead of using a regular, finer grid, the question is whether the refinement can be done more efficiently by using an irregular grid with fewer points. Our goal is to have a fine grid but to avoid expensive computations which are caused by refining the whole area rather than only the parts where we need to add more points. Therefore, we combine this construction method with a refinement algorithm.
4. Grid Refinement Algorithm. As the RBF approximation method is meshfree, i.e. there is no special connectivity between points, the procedure of adding points will be easy and convenient. The proposed algorithm is recursive and uses Voronoi diagrams. In each step, given a grid, we generate a Voronoi diagram for our grid points, then we consider the Voronoi vertices of each cell of this diagram as possible points to be added to the grid. Finally, we run a test on each Voronoi vertex and decide whether we add the point to the grid or not.
We have used Voronoi vertices as new possible points for our grid since they are equidistant to three or more previous grid points, and thus lie "in between" the previous grid points. A brief introduction to Voronoi diagram as well as the strategy of the refinement algorithm will be given in the following two sections, respectively.

Voronoi diagram.
A Voronoi diagram is a geometric structure that divides a d-dimensional space into cells based on the distance between sets of points in the space [24]. Many algorithms have been proposed for computing Voronoi diagrams. The fundamental and most popular ones include: The Divide and Conquer algorithm, and Fortunes's Sweep Line algorithm [1]. For the purpose of explaining the structure of Voronoi diagrams, we will explain a very simple but less efficient algorithm using perpendicular hyperplanes.
Let S = {s 1 , s 2 , . . . , s n } ⊂ R d be a set of n arbitrarily distributed and distinct sites (points) in R d . The perpendicular bisector algorithm works as follows: for each pair of sites in S we construct a perpendicular hyperplane to the line segment joining these sites. At the end of this process, we will have intersections of finitely many hyperplanes which build up cells, with a convex polygon structure, known as Voronoi regions. The boundaries of each region are called Voronoi regions and the intersections of multiple Voronoi edges are called Voronoi vertices. For more details see [1,20,17]. Mathematically, the Voronoi region of a point s i in S is defined by where · denotes the Euclidean distance. This means that for every point x ∈ R d within a Voronoi region V i the Euclidean distance of x to the site s i , which is also inside the region, is smaller than the Euclidean distance of x to any other site s j . Remark 6. As the Voronoi region is the intersection of finitely many hyperplanes, each Voronoi region is a convex polygon.
Delaunay triangulation is the dual structure of Voronoi diagram. The relation between them is simply that the sites of the Voronoi diagram are the vertices of the Delaunay triangulation and vice versa.
Voronoi diagrams and Delaunay triangulations have enormous applications in different scientific fields, especially in mesh generation and nodes insertion procedures. Sibson [26] developed an interpolation method based on Voronoi diagrams, the method is known as natural neighbour interpolation method. Moreover, there is a refinement algorithm called Ruppert's Delaunay refinement algorithm [25]. Some recent works include [18], presenting a refinement procedure to develop the gradient smoothing method using Delaunay triangulation for the adaptive analysis in solid mechanics, and [29], where a Voronoi neighbour criterion is used to construct the adaptive radial point interpolation method. Voronoi diagrams have also been used in kernel-based adaptive particle methods for numerical flow simulation [16], and for a thinning algorithm in multistep interpolation with Radial Basis Functions [5].

4.2.
The refinement algorithm strategy. We will now describe the refinement algorithm in detail.

Fix a compact neighbourhood K ⊂ R d of the equilibrium x 0 and a Radial
Basis Function. Let n = 1 and start with an initial set of grid points X n = {x 1 , x 2 , . . . , x Nn } ⊂ K, not containing any equilibrium. 2. Calculate a Lyapunov function v n using the RBF method with the grid X n . 3. Generate Voronoi vertices, Y n = {y 1 , y 2 , . . . , y Mn } ⊂ R d , for the grid points X n . Exclude points in Y n which are equilibria or which lie outside K. 4. Run a test on each vertex in the set Y n and check whether v n (y j ) < 0 (y j ∈ Y − n ) or v n (y j ) ≥ 0 (y j ∈ Y + n ), where j = 1, . . . , M n and Y n = Y − n ∪ Y + n . 5. Define new grid X n+1 = X n ∪ Y + n . 6. n → n + 1, repeat the steps 2. to 5. until Y + n = ∅. Example 2. Consider the system in Example 1. We will explain the refinement algorithm starting with a coarse equidistant grid X N1 = {(x, y) ∈ R 2 | x, y ∈ {−h, 0, h}} \ {(0, 0)} with h = 1 and the set K = [−1, 1] 2 , i.e. N 1 = 8.
• The first step of the algorithm n = 1: We start with an initial grid X 1 = {x i } 8 i=1 and calculate the Lyapunov function v 1 on the grid X 1 . Figure 4 (a) shows the level set v 1 (x, y) = 0 with the area where v 1 (x, y) > 0 in green. We can also see the Voronoi diagram for our grid points and the circled Voronoi vertices which will all be added to the set X 1 , since they are all located in the green area. In this case, we have Y + 1 = Y 1 and Y − 1 = ∅, i.e. all four new points (the equilibrium is not part of Y 1 ) have positive orbital derivative and will be added to the grid.  • The second step of the algorithm n = 2: Our new set of grid points is now Figure 4 (b) shows the level set v 2 (x, y) = 0. Again, we generate a Voronoi diagram for the set X 2 and determine the Voronoi vertices to be added to the existing grid points, see Figure 4 (b).
• The third step of the algorithm n = 3: The new set of grid points is X 3 = X 2 ∪ {y j } 12 j=1 . After going though the same steps again we show in Figure 5 (a) the grid and level set v 3 (x, y) = 0.
• After the fourth step, the set Y + 4 = ∅ and the algorithm terminates as all Voronoi vertices have negative orbital derivative. Figure 5 (b) shows the final set of grid points X 4 = X 3 ∪ {y j } 12 j=1 and no more areas with positive orbital derivative (green areas). Thus, we have found a Lyapunov function and will now be able to determine sublevel sets, which are subsets of the domain of attraction. We plot the graph of the Lyapunov function in Figure 6 (a) and its sublevel sets in Figure 6    After successfully testing the refinement algorithm on a linear system, we will now examine it on a nonlinear one.    For this example, we approximate the Lyapunov function satisfying V (x) = − x 2 . We have used the Wendland function ψ 6,4 with c=1 and started with the grid X 1 = {(x, y) ∈ R 2 | x, y ∈ {0, ±h, . . . , ±1}} \ {(0, 0)} with h = 0.2, i.e. N 1 = 24 points. Figure 7 (a) shows the region where v 1 (x, y) > 0 (green) and Figure 7 (b) shows the Voronoi diagram and vertices Y 1 . The red points marked with o (14 points) have positive orbital derivative and form the set Y + 1 , these are the points that will be added to the previous grid. Figure 8 (a) shows the level sets of the approximation v 2 , using the refined grid X 2 = X 1 ∪ Y + 1 ; the orbital derivative of v 2 for all points of X 2 is negative by construction.
After four refinement steps, the algorithm terminates and the final function satisfies v 4 (x, y) < 0 everywhere; no green areas (where v 4 (x, y) > 0) occur in Figure 8 (b). To show that v 4 is negative, we have checked that the sign of v 4 (x, y) is negative on the grid X check = {(x, y) | x, y ∈ {0, ±h check , ±2 h check , . . . , ±1}} \ [−0.1, 0.1] 2 with h check = 10 −3 -note that we have excluded a small neighborhood of the equilibrium because of the discussion in Remark 5. Figure 9 (a) shows the graph of the constructed Lyapunov function v 4 (x, y), and Figure 9 (b) shows some of its sublevel sets. To construct a Lyapunov function with a regular grid which has negative orbital derivative on X check , we need a grid of N = 360 points.
The points in the green area (red o) form the set Y + 1 and will be added to our existing set X 1 of grid points; there are two Voronoi vertices lying in the white set (blue *), they form the set Y − 1 and will not be added to the grid. x(x 2 + y 2 − 1) − y(z 2 + 1), y = y(x 2 + y 2 − 1) + x(z 2 + 1), z = 10z(z 2 − 1).
The system has an exponentially stable equilibrium at (0, 0, 0) and its domain of attraction is given by In [6,Example 6.4], an RBF approximation with 137 points resulted in a large area near the equilibrium with positive orbital derivative; larger than the set [−0.2, 0.2] 3 , which is later excluded in our example. With a modified algorithm, using the Taylor polynomial at the equilibrium, this was overcome in [6, Example 6.4].
As generally the domain of attraction is not known, we have chosen to use a grid in the set K = [−0.9, 0.9] 3 , which is not a subset of the domain of attraction, but also does not include other invariant sets. This is a more realistic, but also more challenging test case for the method.
We choose the Wendland function ψ 6,4 with c = 0.6. We have started with a regular grid X N = {(x, y, z) ∈ R 3 | x, y, z ∈ {0, ±h, ±2h, . . . , ±0.9}}\{(0, 0, 0)} and we have approximated the Lyapunov function V = Q satisfying Q = − x 2 . In order to check the sign of the orbital derivative of the constructed Lyapunov    To construct a Lyapunov function v with a regular grid which has negative orbital derivative on X check we need N = 2196 points (h = 0.15).   Now, applying our refinement algorithm with an initial set N 1 = 342 grid points, gives us a final set of grid points N 4 = 458, thus reducing again the number of points needed by a factor 4. The constructed Lyapunov function v has negative orbital derivative on the set X check . Figure 10 (a) shows the initial grid X 1 , and Figure 10 (b) the grid X 4 after the refinement algorithm. Figure 11 Table 2 we have listed the value of h for the initial grid and the corresponding number of points in the initial grid N initial as well as the number of points in the grid after the refinement algorithm has terminated. Moreover, it shows the running time for solving the linear systems in all refinement steps (time 1) and for calculating and plotting the orbital derivative of the constructed Lyapunov function for the final set of grid points (time 2), note that this time is proportional to the number N f inal of grid points in the final grid. For all calculations we have used a standard laptop with an Intel(R) Core (TM) i5-3550 CPU @ 3.30 GHz processor.  and for calculating and plotting the orbital derivative of the constructed Lyapunov function for the final set of grid points (time 2). "Patches" means that after termination of the refinement algorithm there are still areas with positive orbital derivative remaining. The shortest successful construction of a Lyapunov function is achieved with an initial grid of 24 points.
For small (16) and larger (48, 80, 120) numbers in our starting grid, the refinement algorithm stopped, but there were patches of areas where v (x, y) > 0, see, for example, Figure 12. On the other hand, when we started with 24, 36 or 168 grid points, we achieved good results with no patches remaining at the end, see, for example, Figure 13. There were also some cases, for 224 and 288 starting points, where the refinement algorithm did not add any more points but we still had small patches. Finally, with 360 points, we reached the case where the refinement did not add any points and at the same time we did   The reason for these patches is that the last refinement step placed all Voronoi vertices in areas where v (x, y) < 0, thus missing the areas where refinement would be necessary. Thus, a further test to ensure that the orbital derivative is negative is necessary; therefore, we have checked the sign of v on the grid X check = {(x, y) ∈ R 2 | x, y ∈ {0, ±h check , . . . , ±1}} \ [−0.1, 0.1] 2 with h check = 10 −3 . Moreover, future work will include an improvement of the refinement algorithm to add points to the grid in these areas.
The refinement algorithm solves (small) linear systems in each refinement step, while for the grid with 360 points only one (larger) linear system needs to be solved (time 1). Comparing the four successful constructions of a Lyapunov function (24,36,168 and 360 initial points), the starting grids with 24 and 36 have a shorter time 1 than the grid with 360 points, where no refinement step is necessary, so here solving several small systems is shorter than solving a large one. Starting with 168 points, however, takes longer. However, as the number of points in the final grid is proportional to the time to calculate and plot the orbital derivative (time 2), which takes much longer than solving the linear systems, all smaller starting grids (24,36,168) [6]. This table shows the initial sets of regular grid points and the final number of points after refinement. Moreover, it shows the running time for solving the linear systems in all refinement steps (time 1) and for calculating and plotting the orbital derivative of the constructed Lyapunov function for the final set of grid points (time 2). The shortest successful construction of a Lyapunov function is achieved with an initial grid of 168 points, without any refinement step (time 1); the shortest construction and plot (time 1+ time 2), however, is achieved with an initial grid of 36 points.
A(0, 0) = {(x, y) ∈ R 2 | x 2 + y 2 < 1}, so K actually is not a subset of the domain of attraction. As usually the domain of attraction is not known in advance, this is a realistic situation, and even in this situation, the refinement algorithm works well, see  Table 3. Although in this example, the shortest time to solve the linear system(s) is achieved with the grid of 168 points, without any refinement, again the sum of time 1 and time 2 (solving linear system(s) and plotting the orbital derivative) is the longest for the initial grid of 168 points. 4.3.3. Third example. The influence of the starting grid on the refinement algorithm was also investigated in the 3-dimensional system introduced in Example 4. In Table 4, we have presented the values of h and the corresponding number of regular starting grid points, the final number of grid points after the refinement process and the running time in each case.  Again, when we started with a coarse grid (124 initial points), the refinement algorithm stopped but with some remaining patches where v (x, y, z) > 0, see Figure  14 (a). Moreover, starting with a finer grid of 728 or 1330 initial points, the refinement algorithm did not add any more points but we still have patches, see Figure  14 (b) and (c). However, starting with a regular grid of 342 points, and after 4 refinement steps, we ended up with a desirable result with no patches remaining at the end. A similar result was achieved with a regular fine grid of 2196 points, except for a small neighborhood of the origin, see Figure 15. Recall that, for this example, we used X check = {(x, y, z) ∈ R 3 | x, y, z ∈ {0, ±h check , . . . , ±0.9}}\[−0.2, 0.2] 3 with h check = 10 −2 , to check the sign of v is negative everywhere for 342 and 2196 points in the initial grid.
Here, there is a considerable difference between the two successful constructions of more than a factor four, both in the times to solve the linear systems (time 1) and in the time to plot the orbital derivative (time 2).
Summarising, starting with a grid which is too coarse or too fine may end the refinement algorithm with a function that has still areas where the orbital derivative is positive. A moderately coarse grid is the most promising starting point. In most cases, the time to calculate a Lyapunov function was shorter by using the refinement algorithm, even if it involves solving a system of linear equations in each refinement step (time 1). An even greater advantage of the refinement algorithm, however, becomes apparent when using the calculated Lyapunov function further, e.g., to calculate and plot its orbital derivative (time 2), which is proportional to the number of points in the final grid; here, the reduction in the number of grid points pays off considerably. 5. Conclusion and outlook. In this paper we have introduced a refinement algorithm for the Radial Basis Function method to construct Lyapunov functions. The recursive refinement of the grid is done by considering the Voronoi vertices of the previous grid as possible new grid points, which are added if the orbital derivative of the previous approximation is non-negative. The algorithm terminates, if no points are added.
Compared to using a regular grid, the refinement algorithm is able to reduce the required number of grid points by nearly a factor of 4 in some examples and a factor of 2 in others, depending on the system we are solving and which Lyapunov function we approximate (i.e., T or Q ). This reduces the computational effort and time considerably when evaluating the constructed Lyapunov function. It even succeeded in constructing a Lyapunov function if the grid was placed in a set, which was not a subset of the domain of attraction.
When the refinement algorithm terminates, since all Voronoi vertices have already negative orbital derivative, there may still be points in between, where the orbital derivative is positive. Hence, the values of v (x) have been checked to be negative on a very fine grid. Future work will improve this by providing a reliable verification of the negativity of the orbital derivative using error estimates. This can then be used for a modified refinement algorithm, where points in these areas of positive orbital derivative will be added, and the refinement continues until a Lyapunov function is found. Further work on this refinement algorithm could also include the derivation of optimal initial grid configurations. This paper has established the first refinement algorithm for the construction of Lyapunov functions using Radial Basis Functions. It is a considerable improvement from the regular grid, that was used until now, and provides a systematic way to reduce the number of grid points and tackle larger problems in higher dimensions.
(a) The refinement with 124 initial grid points.
(b) The refinement with 1330 initial grid points.
(c) The refinement with 728 initial grid points (4 little green patches). Figure 14. (a) The level set v (x, y, z) = 0 of the orbital derivative after the last refinement step started with 124 points and ended up with 180 points, (b) and (c) the level set v (x, y, z) = 0 of the orbital derivative calculated with 728 and 1330 points, where the refinement did not add more points. In all cases we can see some patches remaining at the end (green), where v (x, y, z) > 0.