DETERMINATION OF THE BASIN OF ATTRACTION OF A PERIODIC ORBIT IN TWO DIMENSIONS USING MESHLESS COLLOCATION

. A contraction metric for an autonomous ordinary diﬀerential equation is a Riemannian metric such that the distance between adjacent solutions contracts over time. A contraction metric can be used to determine the basin of attraction of a periodic orbit without requiring information about its position or stability. Moreover, it is robust to small perturbations of the system. In two-dimensional systems, a contraction metric can be characterised by a scalar-valued function. In [9], the function was constructed as solution of a ﬁrst-order linear Partial Diﬀerential Equation (PDE), and numerically constructed using meshless collocation. However, information about the periodic orbit was required, which needed to be approximated. In this paper, we overcome this requirement by studying a second-order PDE, which does not require any information about the periodic orbit. We show that the second-order PDE has a solution, which deﬁnes a contraction metric. We use meshless collocation to approximate the solution and prove error estimates. In particular, we show that the approximation itself is a contraction metric, if the collocation points are dense enough. The method is applied to two examples.

1. Introduction. Ordinary differential equations model numerous real-world situations and arise from Physics, Biology, Chemistry, Economics and many other areas. Often, the long-time behaviour of such systems is determined by periodic orbits, and it is an important, but non-trivial task to prove the existence and uniqueness of periodic orbits, as well as to analyse their stability. The basin of attraction of an asymptotically stable periodic orbit, consisting of all initial conditions such that the corresponding solutions converge to the periodic orbit, is of particular interest.
Among the classical analytical methods to determine the basin of attraction of a periodic orbit is the method of Lyapunov functions. These, however, require the exact position of the periodic orbit, which in general is not known analytically. A contraction metric is an analytical method to determine the basin of attraction, which does not require any information about the periodic orbit. Moreover, contraction metrics are robust with respect to small perturbations of the system.
A contraction metric is a Riemannian metric, such that the distance of adjacent solutions with respect to the metric decreases over time. The Riemannian metric can be expressed by a matrix-valued function M : R n → R n×n , where M (x) is a symmetric, positive definite matrix, such that v, w M := v T M (x)w defines a point-dependent scalar product for v, w ∈ R n . A contraction metric to determine the basin of attraction of a periodic orbit only assumes contraction of adjacent solutions in the (n−1)-dimensional space perpendicular to f (x), the right-hand side of the autonomous ODEẋ = f (x) under consideration, where x ∈ R n . Contraction metrics can also be used to study the stability and basin of attraction of equilibria, where the contraction is required with respect to any adjacent solution, but in this paper we will focus on the case of a periodic orbit.
In the literature, contraction metrics appear under different names, see [12, Section 2.10] for an overview. In [18], the distance between solutions with respect to a Riemannian metric was studied. For the case of a periodic orbit with the Euclidean metric as contraction metric see [4], and for a general Riemannian contraction metric see [15,14,22], and also [17,3]. Closely related areas include incremental stability, see [19,1], and Finsler-Lyapunov functions [6].
Converse theorems, showing the existence of a contraction metric for exponentially stable periodic orbits have been considered in [7]. It turns out that in 2-dimensional systems, there exists a contraction metric of the specific form M (x) = e 2W (x) I, where W is a scalar-valued function and I the 2 × 2 identity matrix. This is the reason why we focus on 2-dimensional systems in this paper.
A method to explicitly construct a contraction metric for periodic orbits in autonomous ODEs using meshless collocation, in particular Radial Basis Functions, was proposed for 2-dimensional systems in [9]. The contraction metric is of the specific form M (x) = e 2W (x) I, where W is a scalar-valued function; W is characterised as the solution of a first-order PDE and approximated by meshless collocation. The approximation itself, if sufficiently close, is a contraction metric; the paper also includes error estimates. However, the PDE involves quantities that need to be numerically approximated.
For three-and higher-dimensional systems [10], the contraction metric is constructed on and near the periodic orbit by first numerically approximating the periodic orbit and the solution to its first variation equation, namelyẏ = Df (x(t))y, where x(t) is a solution ofẋ = f (x), in this case a periodic solution. Finally, the metric is extended to the remaining basin of attraction by using a Lyapunov function; note that outside a neighborhood of the periodic orbit there exists a contraction metric of the specific form M (x) = e 2W (x) I in any dimension.
In [11], a contraction metric for a time-periodic ODE with periodic orbit was constructed as a matrix-valued, continuous and piecewise affine (CPA) function. The contraction conditions become the constraints of a semidefinite optimisation problem.
The Sums-of-Squares (SOS) method was used to construct a contraction metric to show global stability of an equilibrium of a polynomial system [2], as well as to construct a contraction metric for a periodic orbit of an autonomous system [20].
In this paper, we will improve the approach in [9]. For a 2-dimensional system, we characterise the function W , defining a contraction metric through M (x) = e 2W (x) I, by a second-order PDE which does not involve unknown quantities. We then apply meshless collocation to approximate W and show that the approximation is a contraction metric if the collocation points are sufficiently dense. Parts of this paper are based on the second author's PhD thesis [21].
The outline of the paper is the following: In Section 2 we introduce a contraction metric, state the implications and a converse theorem, and discuss its characterisation as solution of a first-and a second-order PDE. In Section 3, we recall meshless collocation, and apply the general theory to our specific second-order differential operator. The main result of the paper is given in Section 4 where we prove a constructive converse theorem, i.e. we show that the method constructs a contraction metric on any compact subset of the basin of attraction, if the collocation points are dense enough. The paper ends with examples in Section 5 and conclusions. An appendix contains explicit formulas for the computations using meshless collocation.
2. Contraction metric. We consider the autonomous ordinary differential equationẋ where x ∈ R 2 and f ∈ C σ (R 2 , R 2 ), with σ ≥ 1 although further restrictions on σ are specified later. We will state many results in R n for general n, but restrict ourselves to n = 2 where necessary. The solution of the initial value problem with x(t) = ξ defines a dynamical system with flow operator S t ξ = x(t), as we will assume that solutions exist for all t ≥ 0. We are interested in showing the existence and uniqueness of an exponentially stable periodic orbit, and in determining its basin of attraction.
Definition 2.1. A periodic orbit Ω of a dynamical system with flow operator S t is called exponentially stable if it is orbitally stable and there are δ, µ > 0 such that dist(x, Ω) ≤ δ implies dist(S t x, Ω)e µt t→∞ −→ 0. Its basin of attraction A(Ω) is defined by Here, we define dist(x, K) = min y∈K x − y for x ∈ R n and a compact set K ⊂ R n , where · denotes the Euclidean norm in R n .
The orbital derivative, the derivative along a solution of (1), is defined below.
We will use a contraction metric, which is a Riemannian metric with respect to which the distance between adjacent solutions decreases. The distance v is required to decrease only in directions perpendicular to f (x), if we aim to find a periodic orbit. If the metric contracts in all directions, then the attractor is an equilibrium. Definition 2.3. A contraction metric (for a periodic orbit) forẋ = f (x), f ∈ C 1 (R n , R n ) in K ⊂ R n , where K does not contain any equilibrium, is a matrixvalued function M ∈ C 1 (K, R n×n ) such that M (x) is symmetric and positive definite for all x ∈ K and is negative for all x ∈ K. Here, M (x) denotes the matrix with entries ∇M ij (x), f (x) , the orbital derivative, and v, w = For an exponentially stable periodic orbit in a two-dimensional system, there always exists a contraction metric of the form M (x) = e 2W (x) I, where W is a scalar-valued function and I the 2 × 2 identity matrix. Then, L M (x) = L W (x), where the latter function is defined in (2).
We cite [9, Theorem 1.1], which shows that if a contraction metric in the sense above exists in a positively invariant, compact, connected set K, then there exists a unique exponentially stable periodic orbit Ω in K. Furthermore K is a subset of the periodic orbit's basin of attraction A(Ω). Note that L in (3) Theorem 2.4. Let W ∈ C 1 (R 2 , R) and let ∅ = K ⊂ R 2 be a compact, connected and positively invariant set, which contains no equilibrium. Moreover assume Then there exists one and only one periodic orbit Ω ⊂ K. The periodic orbit is exponentially stable, and the Floquet exponent different from 0 is less than or equal to −ν. Moreover, K ⊂ A(Ω) holds.

First-order PDE.
While Theorem 2.4 shows the implications of a function W with the above properties, the following theorem is concerned with the existence of such a function W and its characterisation by a first-order PDE.
Further results, so-called converse theorems, on the existence of general contraction metrics are available, but we are particularly interested in the characterisation by a linear PDE, which can be used to approximate its solution by meshless collocation.
The following theorem and its proof are from [9]. Note that (4) and (2) imply Let Ω be an exponentially stable periodic orbit with period T and Floquet exponents 0 and −ν < 0.
In [9], equation (4) was used to approximate its solution W by meshless collocation, and thus constructing a contraction metric. However, for this approach the constant µ needs to be numerically approximated, which was achieved by numerically approximating the periodic orbit and its non-trivial Floquet exponent.
However, the main advantage of using contraction metrics is that no information about the periodic orbit is required; hence, we will overcome the requirement in this paper.

2.2.
Second-order PDE. By rearranging equation (4) and taking the orbital derivative on both sides we can eliminate the term µ from the collocation equation. This means that the periodic orbit will not need to be numerically approximated to create the right-hand side of the collocation equation. The equation, however, will be a second-order PDE, and thus the collocation will require the solution W , and thus f , to be sufficiently smooth.
Let Ω be an exponentially stable periodic orbit.
Then there exists a function W ∈ C σ−2 (A(Ω), R) such that where L is given by (3) and D is the following linear second-order differential operator Conversely, a solution W of (5) satisfies equation where µ is defined as in Theorem 2.5.
Proof. From Theorem 2.5 we have Rearranging (8), noting that f (x) = 0 for all x ∈ A(Ω), and taking the orbital derivative gives as µ is constant. Then the quotient rule and multiplying through by f (x) 2 gives which shows (5). Now assume that W satisfies (5). By rearranging the equation, we obtain for all x ∈ A(Ω).
Integration along the solution S τ x from 0 to t gives which we can write as with .

PETER GIESL AND JAMES MCMICHEN
We show the equation (7) first for x ∈ Ω. For x = p ∈ Ω we have (9) = c p T F + νT by definition of F and using 1 g., see [9, Part I of the proof of Theorem 2.1]. This shows c p = − ν F = −µ, and thus (9) shows (7) for all points on the periodic orbit Ω.
Now let x ∈ A(Ω)\Ω. In particular, for every p ∈ Ω there is a sequence t n n→∞ −→ ∞ with lim n→∞ S tn x = p. Then (9) with t = t n shows, letting n → ∞, that as all functions are continuous. As W also satisfies this equation with −µ instead of c x , as shown above, we can conclude c x = −µ, which shows (7).
We now give an explicit form for the differential operator defined in (6).
Proof. We have In the next section we will approximate the solution of (5) using meshless collocation.
3.1. General method. In this section we recall the main facts from meshless collocation. We first introduce Reproducing Kernel Hilbert Spaces, and then formulate the generalised interpolation problem, corresponding to the PDE that we seek to solve. We state an error estimate that will be used in the next section. We follow [13, Section 2], see also [24].
Let O ⊂ R n be a bounded set with Lipschitz continuous boundary. A Reproducing Kernel Hilbert Space (RKHS) is a Hilbert space H of functions g : O → R with inner product ·, · H such that the following properties hold with a kernel Φ : O × O → R: We want to solve the following generalized interpolation problem: Given N linearly independent functionals λ 1 , . . . , λ N ∈ H * , where H * denotes the dual of H, and N numbers r i ∈ R, i = 1, . . . , N , find the norm-minimal interpolant w ∈ H satisfying It is well known that there is a unique norm-minimal interpolant, which is a linear combination of the Riesz representers of the functionals, and the coefficients can be determined by solving a system of N linear equations. If H is a RKHS, then we have the following formula for the norm-minimal interpolant w: where the superscript y in λ y j denotes the application of the functional with respect to y. Note that α ∈ R N is the solution of the linear system Aα = r, where r = (r i ) i=1,...,N ∈ R N and A = (a ij ) i,j=1,...,N ∈ R N ×N is given by The matrix A is positive definite, since the functionals are assumed to be linearly independent.
In the following we consider the Sobolev space W τ 2 (R n ) with τ > n/2, which is a RKHS. While the reproducing kernel is rather complicated, there is a reproducing kernel, defined by a Wendland function, see Definition 3.1, which generates the same Hilbert space, but with a different, yet equivalent norm. We will still call such a kernel a reproducing kernel of the RKHS. Definition 3.1 (see [23]). Let k ∈ N 0 and l ∈ N. We define x + = x for x ≥ 0 and x + = 0 for x < 0.
From now on, we choose a Wendland function with smoothness degree k ≥ 3 and choose Φ to be the kernel of the RKHS H = W τ 2 (R n ) as above, with τ = k + n+1 2 . Now let us apply the method to solve a PDE of the form where L is a differential operator of the form c β : O → R and m ∈ N, where r(x) is a given function. We assume that τ = k + n+1 2 > m + n/2 and that (12) has a solution W ∈ H = W τ 2 (R n ). We call a point x ∈ R n a singular point of L if c β (x) = 0 for all |β| ≤ m.
We choose collocation points X N = {x 1 , . . . , x N } and define the functionals λ j = δ xj • L to find the norm-minimal interpolant with the data r j = r(x j ).
The following proposition shows that the functionals are linearly independent if the collocation points x j are no singular points of L, see [13, Proposition 3.3].
Proposition 3. Let Φ be a reproducing kernel of W τ 2 (R n ) with τ > m + n/2 and L be as above. Let X N = {x 1 , . . . , x N } be pairwise distinct points, which are not singular points of L. Then, the functionals λ j = δ xj • L are linearly independent over W τ 2 (R n ).
Finally, we cite an error estimate, see [13,Corollary 3.6], for the kernel given by a Wendland function with smoothness index k. The error estimate is in terms of the fill distance, which measures how dense the collocation points X N lie in O.  where C is a constant independent of W and h = sup x∈O min xj ∈X N x − x j . Moreover, we have 3.2. Approximating W . The method above was used in [9] to construct a Riemannian metric by using the first-order differential operator D 1 w(x) := W (x) and approximating the PDE (4). We, however, will use the second-order PDE (5) with differential operator L = D, see (6), so m = 2 and right-hand side is given by Let us now restrict ourselves to the case n = 2. Let O be a bounded domain with Lipschitz continuous boundary and O ⊂ A(Ω). We will use meshless collocation with Wendland functions as discussed in Section 3.1. Let k ≥ 3 and σ ≥ k + 3.5, which in particular means σ ≥ 6.5. We check the assumptions of Theorem 3.2: note that , which can be seen from (10). The solution of (5) satisfies W ∈ C σ−2 (O) ⊂ W Recall that ψ 0 (r) = φ l,k (cr), where c > 0 and φ l,k denotes the Wendland function with l = k + 2, see Proposition 3.1. The vector α ∈ R N is the solution of the linear system Thus, we consider a δ-neighborhood Ω δ of the periodic orbit, such that the time that all solutions enter Ω δ is finite, which enables us to prove the estimate for the other points. Then for all > 0 there exists a h( ) > 0 such that for all sets X N ⊂ O of N pairwise distinct collocation points with fill distance h ≤ h( ) we have where w denotes the norm-minimal interpolant and W the solution of (5).
Hence, the right-hand side of (19) is smaller than f 2 m 0 . Thus using the definition of 0 . Thus, Now we fix a point on the periodic orbit p ∈ Ω and integrate along the solution S τ p from 0 to t with t ∈ [0, T ), Since W (x) = −µ f (x) − L(x), see (7), and t < T , we have , Integrating along the periodic orbit and using 1 T T 0 L(S t p)dt = −ν, e.g., see [9, Part I of the proof of Theorem 2.1], we get Equation (20) can also be rearranged to give Following a similar process yields This gives which shows, in particular, (17) for every x = p ∈ Ω.
3. The estimate on O \ Ω Fix an arbitrary z ∈ O \Ω and denote t 1 ≥ 0 as the first time that dist(S t1 z, Ω) = S t1 z − p ≤ δ with suitable p ∈ Ω; note that t 1 ≤ T 1 . We wish to show that the following estimate holds true Let D 1 be the functional denoting the first order orbital derivative. Define λ := (δ St 1 z • D 1 ) and γ := (δ p • D 1 ) and note that λ, γ ∈ H * . Then we have using (14).

PETER GIESL AND JAMES MCMICHEN
Then by following a similar methodology to the proof of Theorem 3.35 from [8] we can bound λ − γ H * by the following, see (11), Denoting r := S t1 z − p and using Taylor's Theorem there is anr ∈ [0, r] such that Furthermore due to the choice of δ made in (18) and r ≤ δ we have (25) By using (25), we get from (24) that This shows (23). Now we show that (17) holds true for z ∈ O \ Ω. From (20) and the definition of for all x ∈ O. Then, integrating along the orbit through z ∈ O \ Ω, we obtain This becomes, Rearranging and multiplying through by − f (z) gives In a similar way we can show We have using (22) and (23). Putting (29) into (28) shows (17).
Finally we show that w defines a contraction metric as L w is negative. Then the norm-minimal interpolant w of W with respect to the grid X N satisfies Proof. We apply Theorem 4.1 to W and w. By (17) and (7) we have for all Hence for all x ∈ O, we have

5.
Examples. We apply the method to two systems: the first is a system for which we can derive an analytical formula for L M (x). This allows us to calculate the error in the approximation. The second example is the Van-der-Pol system. For this section we introduce the notation x = (x, y) T .
The system has an equilibrium at the origin and an exponentially stable periodic orbit Ω := {(x, y) ∈ R 2 | x 2 + y 2 = 1}, A(Ω) = R 2 \ {0}. We have µ = 2 and therefore we know the true value of W which is given by We define O := {x ∈ R 2 | 0.5 < x < 1.5}. It is clear that O is positively invariant. The collocation points X N form a hexagonal grid as was shown is optimal in [16]. The fill distance is h = 0.1 The collocation points can be seen in Figure  1, left, along with the periodic orbit Ω and four sample trajectories. We have used the kernel given by ψ 0 (r) = φ 6,4 1 5.5 r . In Figure 1, right, we show some level sets of L w (x), which show that L w (x) < 0 for all x ∈ O. In Figure 2, the approximation error L w (x) − L W (x) is shown.

5.2.
Van-der-Pol equation. We study the Van-der-Pol equation, which was also considered in [9], We approximate W using the kernel given by ψ 0 (r) = φ 6,4 (r). We have used 3583 collocation points, forming a hexagonal grid with fill distance h = 0.1 √ 5 2 ; the collocation points and the level set L w (x) = 0 are shown in Figure 3, left. In Figure  3, right, we highlight the set {x ∈ R 2 | L w (x) > 0} (green); a positively invariant lies in the set {x ∈ R 2 | L w (x) < 0}, bounded by the L w (x) = 0 (magenta). The set K thus contains a unique exponentially stable periodic orbit (black) and is a subset of its basin of attraction.
set within the white set {x ∈ R 2 | L w (x) < 0} is a subset of the basin of attraction of Ω.
To determine a positively invariant set K we follow the method given in [9]. We approximate the solution of the equation V (x) = −1 using meshless collocation with the kernel given by ψ 0 (r) = φ 5,3 r 5.5 and 602 collocation points. In Figure 4, left, the collocation points and the level set v (x) = 0 (red) are shown and we see that there are small areas near the periodic orbit where v (x) > 0, as V (x) = −1 has no solution on the periodic orbit. However, the set K = {x ∈ R 2 | v(x) ≤ −0.5} (blue), see Figure 4, left, is positively invariant as v (x) < 0 holds for all x ∈ ∂K. Moreover, the set K is a compact and connected set that does not contain any equilibria and is fully contained within the set where L w (x) < 0, see Figure 4, right. Hence, by Theorem 2.4, the set K contains a unique exponentially stable periodic orbit Ω and K ⊂ A(Ω).
Note that we used more collocation points than the method in [9]; however, we used a uniform grid, while [9] used an adaptive one. Moreover, we did not rely on any quantities of the periodic orbit being known to calculate the approximation.
Conclusions. In this paper we have presented a method to construct a contraction metric for a general 2-dimensional autonomous ODE with an exponentially stable periodic orbit. The existence of a contraction metric implies existence and uniqueness of an exponentially stable periodic orbit, gives an estimate of its non-trivial Floquet exponent and determines a subset of its basin of attraction. A contraction metric does not require any information about the periodic orbit.
A contraction metric is a Riemannian metric, with respect to which adjacent solutions decrease. It can be described by a point-dependent scalar product v, w x := v T M (x)w for v, w ∈ R n , defined by a matrix-valued function M (x), where M (x) is symmetric and positive definite. Converse theorems on the existence of a contraction metric are available in any dimension; in two-dimensional systems, however, the contraction metric can be chosen to be of the special form M (x) = e 2W (x) I with a scalar-valued function W . In dimensions 3 and higher, there are counterexamples showing that no contraction metric of this particular type exists. Hence, the construction of a contraction metric in dimension 2 is easier to achieve, as we only need to construct a scalar-valued function rather than a matrix-valued function. Hence, the method of this paper is not directly transferable to construct a contraction metric in higher dimensions, which will be a subject for future work.
Previous work has shown that for two-dimensional systems W can be characterised as the solution of a first-order PDE and can be approximated using meshless collocation. However, the first-order PDE required information about the periodic orbit. This paper has improved the method by showing that the solution W of a second-order PDE defines a contraction metric through M (x) = e 2W (x) I, and the second-order PDE does not require any information about the periodic orbit. We have applied meshless collocation to approximate W by an approximation w. The approximation w itself defines a contraction metric, if the collocation points are dense enough, which was shown using an error estimate.
Appendix A. Explicit formulas. In this appendix we compute explicit formulas for the orbital derivative of w, where w is given by (15). We denote the operators of the first and second orbital derivative by D 1 W := W , D 2 W := W , respectively. Then w becomes We have D y 1 ψ 0 ( x − y ) = −ψ 1 ( x − y ) x − y, f (y) , and