Determination of the area of exponential attraction in one-dimensional finite-time systems using meshless collocation

We consider a non-autonomous ordinary differential equation over a finite time interval \begin{document}$[T_1,T_2]$\end{document} . The area of exponential attraction consists of solutions such that the distance to adjacent solutions exponentially contracts from \begin{document}$T_1$\end{document} to \begin{document}$T_2$\end{document} . One can use a contraction metric to determine an area of exponential attraction and to provide a bound on the rate of attraction. In this paper, we will give the first method to algorithmically construct a contraction metric for finite-time systems in one spatial dimension. We will show the existence of a contraction metric, given by a function which satisfies a second-order partial differential equation with boundary conditions. We then use meshless collocation to approximately solve this equation, and show that the resulting approximation itself defines a contraction metric, if the collocation points are sufficiently dense. We give error estimates and apply the method to an example.


1.
Introduction. When analysing ordinary differential equations, classical stability concepts study the behaviour of solutions as time tends to infinity. In applications, however, one is often more interested in the development over a finite time interval.
In [16,17], weaker concepts of attraction in finite-time systems were introduced. Here, the notion of attraction does not require the distance between adjacent trajectories to decrease at every point in time.
The concept of finite-time areas of exponential attraction was introduced in [8], using a concept of attraction that requires the distance between adjacent trajectories to decrease from the start to the end time of the finite-time interval, and hence, allows for the distance between trajectories to increase for parts of their evolution. It is shown that the existence of a contraction metric is both a sufficient and necessary condition for a set to be an area of exponential attraction. A contraction metric is a Riemannian metric, with respect to which the distance decreases monotonously over the whole time interval.
In this paper, we limit the spatial dimension to one as then the contraction metric can be characterised by a scalar-valued function W which satisfies a first-order partial differential equation, involving the (unknown) rate of exponential attraction of a solution. By using the fact that this rate is constant along solutions, we show that W can be characterised by a second-order PDE with boundary conditions, involving only known quantities. This boundary value problem is then used to approximate the contraction metric using meshfree collocation. We show error estimates as well as how the approximated contraction metric can be used to determine a subset of the area of exponential attraction. Parts of this paper are based on the second author's PhD thesis [14].
Let us give an overview over the paper: Section 2 recalls the definition of an area of exponential attraction and discusses how it is characterised by a contraction metric given by a scalar-valued function, which satisfies a first-order PDE. In Section 3 we show how the scalar-valued function can be characterised by a second-order PDE boundary value problem and discuss how an approximate solution to the boundary value problem determines a subset of the area of exponential attraction. Section 4 recalls meshless collocation, which is then used to numerically approximate a solution of the boundary value problem. Finally, we apply the method to an example in Section 5 and make conclusions in Section 6.
2. Area of exponential attraction. In this section we recall some definitions concerning attraction in finite-time dynamical systems, in particular the area of exponential attraction. This section is included for the convenience of the reader and, apart from Theorem 2.5, does not contain new results.
We define a finite-time contraction metric and cite a result showing that the existence of a contraction metric is a sufficient condition for a connected set to be an area of exponential attraction. We then show a converse theorem, proving that a contraction metric exists and can be characterised by a first-order PDE. In the following Section 3, we will characterise the contraction metric by a secondorder PDE, which has the advantage that its right-hand side does not contain any unknown quantities.
2.1. Preliminaries. We consider a non-autonomous differential equation of the formẋ We assume that solutions to the initial value problem with x(T 1 ) = x 0 exist for all t ∈ I; they are unique, see e.g. [4,Theorem 1.2]. We define ϕ : I × I × R n → R n , where ϕ(t, t 0 , x 0 ) denotes the solution x(t) of the initial value problem of system (1) with initial value x(t 0 ) = x 0 . Let · denote the Euclidean norm in R n and denote by B η (x 0 ) = {x ∈ R n | x − x 0 < η} the open ball with centre x 0 . We define the Hausdorff distance between two sets A, B ⊂ R n by dist(A, B) = sup x∈A inf y∈B x − y . The theory also works for any other fixed norm in R n ; note that the definitions of finite-time attractivity and the area of exponential attraction depend on the choice of the norm, which is in contrast to the classical stability concepts.
Let us recall the following definitions from [8].
2. G is called a non-autonomous set if G(t) = ∅ for all t ∈ I; furthermore, G is called connected, compact and open if all the t-fibres are respectively connected, compact and open. 3. A non-autonomous set G is called positively invariant with respect to (1), if ϕ(t, τ, G(τ )) ⊂ G(t) and invariant if ϕ(t, τ, G(τ )) = G(t), for all t > τ and t, τ ∈ I.
The following definition from [8] compares the distance between trajectories at the end time T 2 to the start time T 1 and allows the distance between nearby trajectories to increase for part of their evolution. It is thus weaker than, e.g. the concept of finite-time attraction from [11], where solutions are required to reduce the Euclidean distance between nearby trajectories at all times. 1. Let µ : I → R n be a solution of (1). The solution µ is attractive on I if there exists η > 0 such that and the negative number is the rate of exponential attraction.
The area of exponential attraction differs from traditional concepts such as the domain of attraction as it is not related to a single solution. This is because in the classical case of asymptotic attraction we are usually interested in an invariant set that all other solutions will move towards. However, finite-time systems lack sufficient time for the long-term behaviour to emerge and as such individual solutions do not play a dominant role in their study.   ∂x (µ(T 1 )) < 1. The rate of exponential attraction is given by Proof. The first formula follows from [8,Proposition 2.5], restricted to the case where the phase-space is one-dimensional, i.e. n = 1. To show the second formula, consider the first variation equation associated with the solution µ(t), namelẏ Denote the (fundamental) solution of (3) as φ : I → R, then the solution will take the form φ(t) = exp t T1 f x (τ, µ(τ ))dτ . We use the fact that the solution of the first variation equation satisfies φ(t) = ∂Ft ∂x (µ(T 1 )) (see [12, Theorem 3.1 and

2.2.
Sufficiency. To give a sufficient condition for the determination of an area of exponential attraction, we consider a contraction metric, i.e. a metric with respect to which the distance between adjacent solutions decreases at all times. In general R n , a contraction metric is given by a matrix-valued function x) ∈ R n×n denotes the Jacobian with respect to x and M (t, x) ∈ R n×n denotes the matrix with entries M ij (t, x). As we are interested in the case of n = 1, our time-varying Riemannian metric will be a positive scalar-valued function and we can assume that M (t, x) = e 2W (t,x) , where W ∈ C 1 (I × R, R). In this case, a simple calculation shows The following theorem shows that a contraction metric can be used to determine an area of exponential attraction. The theorem follows from [8, Theorem 4.2, Corollary 4.3], adapted to the case n = 1, as all connected intervals in R are convex.
Theorem 2.4. Consider the differential equation (1) with n = 1, let G ⊂ I × R be a nonempty, connected, compact and invariant non-autonomous set and let W ∈ x) and assume that there exists a ν > 0 such that Then G is an area of exponential attraction. Furthermore, i.e. −ν is an upper bound on the rate of exponential attraction of all solutions in G.

2.3.
Necessity: First-order PDE. We will derive a converse theorem, i.e. given an area of exponential attraction, we prove the existence of a contraction metric in Theorem 2.5. This is a different result than the converse theorem [8, Theorem 5.1]: firstly, by restricting ourselves to n = 1 we can assume that M (t, x) is given by e 2W (t,x) and thus, secondly, characterise W as the solution of a partial differential equation rather than an inequality. This will later enable us to approximate W and thus M .
Theorem 2.5. Consider the differential equation (1) with n = 1 and f ∈ C σ (I × R, R) with σ ≥ 2. Let the compact non-autonomous set G ⊂ I × R be an area of exponential attraction with −ν < 0 being the maximal rate of exponential attraction of all solutions in G.
Then there exists a unique function W : where γ(x) is the exponential rate of attraction of the solution ϕ(t, T 1 , x) on I with initial condition x. Moreover, we have Proof. By integrating equation (5) along a solution ϕ(t, This defines the unique W , which satisfies equation (5) and the boundary condition W (T 1 , x) = 0 for all x ∈ G(T 1 ). The other boundary condition is also fulfilled as Note that the function W (t, x) is C σ−1 as all quantities on the right-hand side of (6) are. Since G is an area of exponential attraction with maximal rate −ν, we conclude the bound on L M (t, x).
3. Second-order PDE. Equation (5) in Theorem 2.5 gives an equation involving W , however it cannot yet be used to construct an approximation to W as γ is unknown. We seek to find an equation that only involves a linear differential operator applied to W and known values.
By taking another derivative, the equation becomes a second-order PDE, involving only known quantities. In this section, we will prove that such a function W exists and satisfies a second-order PDE with boundary conditions and, conversely, that a function W with these properties defines a contraction metric. In Section 4 we will approximate the solution to the second-order PDE boundary value problem.
3.1. Necessity. The following corollary to Theorem 2.5 shows that W satisfies a second-order PDE. Corollary 1. Let the conditions of Theorem 2.5 be satisfied with σ ≥ 3. Then the following statements are true: there exists a function W ∈ C σ−1 (G, R) with W (T 1 , x) = 0 for all x ∈ G(T 1 ) and W (T 2 , x) = 0 for all x ∈ G(T 2 ), (7) and Proof. The boundary conditions follow from Theorem 2.5. Equation (7) follows from (5) by taking the derivative with respect to t and noting that γ does not depend on t.
For the last statement we take the orbital derivative of (4) and get Then by using (7) we immediately see that (8) is true.
Corollary 1 provides us with a partial differential equation which can be numerically solved to approximate W ; this approximation can be used to approximate of the area of exponential attraction. Bringing together the results of Theorem 2.4 and Theorem 2.5, we show that the sublevel set {(t, x) | L M (t, x) < λ} with λ < 0 is an area of exponential attraction. The following theorem confirms that the function W that satisfies W = −f x on the interior of the domain and W = 0 on the temporal boundary is the same function as W from Theorem 2.5 and that it characterises an area of exponential attraction.
Let the compact non-autonomous set G ⊂ I × R be an area of exponential attraction with maximal rate of exponential attraction −ν < 0.
Conversely, consider a nonempty, connected, compact and invariant non-autonomous set G ⊂ I×R and assume the function W ∈ C σ−1 (G, R) solves the boundary value problem Then W is the function described in Theorem 2.5. Moreover, L M (t, T 1 , x) = 0 for all x ∈ G(T 1 ) and for all t ∈ I, where M (t, x) = e 2W (t,x) .
Fix λ < 0 and let be a connected non-autonomous set. Then G 1 is an area of exponential attraction and λ is an upper bound on the rate of exponential attraction of all solutions in G 1 .
Proof. The first statement follows from Corollary 1. Now we fix x 0 ∈ G(T 1 ) and let µ(t) = ϕ(t, T 1 , x 0 ). By integrating W (t, x) along µ(t) we obtain the following We integrate (10) again from T 1 to T 2 and obtain This gives using (2) of Proposition 1. Equations (10) and (11) show that the function W given in this corollary satisfies W (t, µ(t)) = −f x (t, µ(t)) + γ(x 0 ). It is therefore the (unique) function from Theorem 2.5. Furthermore, as L M (t, µ(t)) = γ(x 0 ) and γ(x 0 ) does not depend on t, we have L M (t, µ(t)) = 0. The last statement of the corollary follows immediately from Theorem 2.5 where it is shown that L M (t, µ(t)) = γ(x 0 ) is the rate of exponential attraction of the solution ϕ(t, T 1 , x 0 ).

Approximate solution.
We will later approximate the solution of the boundary value problem W by the approximation w and show error estimates on w (x)− W (x) in the interior as well as w(x) − W (x) on the temporal boundary. Hence, we seek to conclude that if these quantities are small, then also L m is close to L M , where m(t, x) = e 2w(t,x) and M (t, x) = e 2W (t,x) . The first step is the following lemma, see [7,Lemma 4.1].
Using Lemma 3.2 we now bound the error between L m , using our approximation, and the rate of exponential attraction, γ(µ(T 1 )) which is also equal to L M (t, µ(t)) for all t ∈ I. Theorem 3.3. Let G be a non-autonomous set and f ∈ C σ (I × R, R) with σ ≥ 3. Let W ∈ C 2 (G, R) satisfy the boundary value problem (9) and let w ∈ C 2 (G, R) be a function satisfying Let µ(t) be a solution to (1), satisfying (t, µ(t)) ∈ G for all t ∈ I, and denote by γ(µ(T 1 )) its rate of exponential attraction.
In the next corollary, we will show how to use a sufficiently good approximation w to W to determine a subset of the area of exponential attraction.
Corollary 2. Let G be a non-autonomous set and f ∈ C σ (I × R, R) with σ ≥ 3. Let W ∈ C 2 (G, R) satisfy the boundary value problem (9) and let w ∈ C 2 (G, R) be a function satisfying as well as the sublevel set H λ := {(t, x) ∈ G | L m (t, x) < λ} for λ ∈ R. Let ν > 0 such that H −ν+ * ⊂ G. Then H −ν− * is a subset of an area of exponential attraction of rate −ν.

Meshless collocation. Meshless collocation, in particular with Radial Basis
Functions, is a powerful method to solve linear PDEs, for more background see [3,20]. Section 4.1 is included for the convenience of the reader and recalls known results, in particular error estimates, following [9]. In Section 4.2 we will then apply the method to our particular case. This will provide us with estimates on Let Ω ⊂ R d be a bounded domain with a piecewise C τ ,s boundary ∂Ω (for more details see [21]), where τ ∈ R + , τ > d/2 and s = τ − τ ∈ [0, 1). Let ∅ = Γ ⊂ ∂Ω. We will later apply the general method to Ω ⊂ I × R ⊂ R 2 with Γ ⊂ {T 1 , T 2 } × R.
Let L be a linear differential operator of the form Lu(x) = |α|≤m c α (x)D α u(x). A singular point of the operator L is a point x ∈ R d such that c α (x) = 0 for all |α| ≤ m. Denote L 0 = id. To solve the boundary value problem we choose two sets of pairwise distinct collocation points, namely X 1 = {x 1 , . . . , x N } ⊂ Ω and X 2 = {y 1 , . . . , y M } ⊂ Γ ⊂ ∂Ω. We consider the Sobolev space W τ 2 (R d ), which is a reproducing kernel Hilbert space (RKHS), and we will use a kernel, given by a Wendland function. Definition 4.1 (see [19]). Let k ∈ N 0 and l ∈ N. For r ∈ R + 0 we define by recursion the Wendland function φ l,k , where x + = x for x ≥ 0 and x + = 0 for x < 0, through The following proposition follows from [6,Proposition 3.11] and the arguments in the proof.
Then Ψ ∈ C 2k (R d × R d , R) is a reproducing kernel of a RKHS which is norm- From now on, we choose a Wendland function with fixed k, called smoothness index, and consider Ψ and W τ (R d ) as above with τ = k + d+1 2 . We assume, moreover, that the boundary value problem (13) has a solution u ∈ W τ 2 (R d ). The mixed ansatz for the approximation s (which is the norm-minimal interpolant of the data in the RKHS) of the function u is given by where δ x0 f (x) = f (x 0 ) denotes Dirac's δ-distribution and the superscript y denotes the application of the operator with respect to the variable y.
The coefficient vector (α, β) ∈ R N +M is determined by the interpolation conditions which are equivalent to the linear system where The matrix A C C T B is positive definite as the points in X 1 are no singular points of the operator L, see [7, Lemma 2.6], and the points in X 2 are no singular points of L 0 . We then have the following error estimate, see [9,Corollary 3.12], where h X1,Ω = sup x∈Ω min xj ∈X1 x − x j and h X2,Γ = sup y∈Γ min yj ∈X2 y − y j denote the meshnorms of the collocation points.
Then, for sufficiently small mesh-norms, we have where C is a constant independent of u.

4.2.
Application to second-order boundary value problem. We will now apply the general theory to our specific problem. Fix k ∈ N, k ≥ 3. Let Ω ⊂ (T 1 , T 2 ) × R be a bounded, open non-autonomous set with piecewise C k+1,1/2 -boundary. Denote the temporal boundary of Ω by Γ := ∂Ω ∩ ({T 1 , T 2 } × R). We consider the boundary value problem Let us define our collocation points, the collocation points in Ω, namely X 1 := {x i = (t i , x i ) ∈ Ω for i = 1, ..., N }, and the collocation points on the temporal boundary X 2 := {ỹ j = (τ j , y j ) ∈ Γ for j = 1, ..., M }. The ansatz for the approximation is with the second-order differential operator Lw(t, x) = w (t, x), where denotes the orbital derivative andx = (t, x). For explicit formulas of the terms in equation (20) and the collocation matrix (16), see [7,Appendix].
We solve the linear system (16) to determine the coefficient vector (α, β), which can be used to compute w(t, x) and w (t, x). Then, by using the approximation of w (t, x) to determine level sets of L m (t, x), we can approximate areas of exponential attraction. We apply Theorem 4.2 to our case to obtain the following error estimates. (Ω) satisfy (19) and let w be the approximant to W as described in Section 4.1.
The above theorem provides us with bounds on W (t, x) − w (t, x) L∞(G) and W (t, x) − w(t, x) L∞(Γ) , as required by Theorem 3.3 and Corollary 2. Indeed, for a non-autonomous set G ⊂ I × R, define Ω := G • . As W, w ∈ C 2 (G, R), the first estimate in Lemma 4.3 also holds for W (t, x) − w (t, x) L∞(G) .

5.
Example. We consider a non-autonomous Bernoulli-type problem with x ∈ R on the time interval [0, T 2 ]. Note that in this example the distance between trajectories both increases and decreases along the time domain. The system iṡ The solution is of the form Since t 0 exp − 2 3 τ 3 + 2τ 2 − 3 2 τ dτ cannot be analytically solved we cannot determine an analytic form of L M .
This example was used to approximate the domain of attraction in [7] with T 2 = 2. The domain of attraction of a solution is defined below; note that by [8, Theorem 6.1] the area of exponential attraction is a subset of the domain of attraction of a solution in it.
and G µ is the maximal set such that this is true (with respect to set inclusion) and which contains the graph of µ.  Figure 1. The collocation points X 1 and X 2 as well as some numerically computed solutions of the system (21) with T 2 = 2.
We have used the kernel given by ψ 0 (r) = φ 6,4 (0.5r) with the Wendland function φ 6,4 and the following collocation points for the time interval with T 2 = 2 We show X 1 and X 2 as well as some numerically approximated solutions in Figure  1. Figure 2 shows the function L m (t, x) and Figure 3 shows some level sets of L m (t, x). Using Corollary 2, H −ν− * is a subset of an area of exponential attraction of rate −ν if H −ν+ * ⊂ G; hence, the sublevel set H 0 is an approximation of a subset of the area of exponential attraction. Since L m is a good approximation of the function L M in Theorem 3.1, see Theorem 3.3, the sublevel set H 0 is actually a good approximation of the area of exponential attraction; the intersection of ∂H 0 with G 0 are the points (0, ±0.1789). Note that H 0 is a (small) subset of the domain of attraction of the zero solution, see [7,Example 1], which is to be expected, cf. Theorem 6.1 [8].
When approximating the domain of attraction for a solution using the method given in [8], one cannot determine information about the domain of attraction in a small neighbourhood around the solution. However, the method presented here does not have that problem so can in particular be used to show that the neighbourhood around the solution lies in the domain of attraction for that solution. Now we study this example with different end times T 2 ∈ {0.4, 0.8, 1.2, 1.6, 2, 2.4}. Considering the system (21) linearised around x = 0, we observe that Hence, the linearised system is asymptotically stable except for T 2 = 3 2 . The area of exponential attraction is approximated by the zero level set of L m in Figure 4; as T 2 increases from 0 to 1.5, the area in x-direction shrinks, and it grows again for T 2 > 1.5.
6. Conclusion. In this paper we have considered a non-autonomous ordinary differential equation on a finite time interval in one spatial dimension. To determine the area of exponential attraction, consisting of exponentially attracting solutions, we have used a contraction metric. We have shown that the contraction metric can be characterised by a scalar-valued function, which satisfies a second-order PDE with boundary conditions, and we have solved this boundary value problem numerically using meshless collocation. We have derived error estimates and have applied the method to an example. The method in the form described above is only applicable to systems in one spatial dimension. The reason is that in these systems there exists a Riemannian metric of the form M (t, x) = e 2W (t,x) , given by a scalar-valued function W ; this is not true for higher-dimensional systems. This is similar to autonomous systems, where a Riemannian metric given by a scalar-valued function only exists in dimension up to 2 in general. For higher-dimensional systems, the Riemannian metric is a matrix-valued function M (t, x). To apply a similar methodology to this paper, we would need to find a PDE which is satisfied by M (t, x), involving only known quantities. Once such a PDE is found, one can again employ meshless collocation to solve matrix-valued PDEs, see [10].  Figure 4. Zero level sets of L m (t, x) for different values of T 2 , namely 0.4 (black), 0.8 (blue), 1.2 (red), 1.6 (green), 2 (magenta) and 2.4 (cyan). The 0-level set of L m is an approximation of the area of exponential attraction. The size of the area of exponential attraction in x-direction shrinks until T 2 = 1.5 and then grows again.