ON THE EXPANSION SCHEMES IN TRAJECTORY REVERSING METHOD

. The paper deals with certain expansion schemes in trajectory rever-sing method for estimating asymptotic stability region of nonlinear dynamical systems. The asymptotic behavior of the sequence of estimates is investigated. Some numerical examples are given.


INTRODUCTION
The trajectory reversing method it seems to be one of the most powerful method for estimating the asymptotic stability region of autonomous nonlinear dynamical systems.The first papers concerning this method were due by Genesio, Tartaglia and Vicino [4], [5] and also by Hsu [8].There are two main ways in which a concrete implementation of this idea may be done.
First.The boundary of stability region is synthetized from a number of system trajectories obtained by backward integration of the differential system which describe the dynamical system, starting from the equilibrium points.These trajectories, starting in a neighborhood of an asymptotic stable point, tend to the boundary of stability region, while the trajectories, starting near an equilibrium point on the boundary, remain related to the boundary and give essential information about it [4], [5].In [2] such a procedure is based on topological properties of the equilibrium points and closed orbits on the stability boundary.Several necessary and sufficient conditions are given to determine whether an equilibrium point or a closed orbit is on the stability boundary.
Second.The stability region (or its boundary) is approximated by a sequence of estimates, consisting of certain domains (or surfaces) around of the stable equilibrium point.Starting from an initial estimate Ω 0 , inside of the true stability region, one performs a backward integration and obtains a new estimate Ω 1 .If Γ 0 denote the boundary surface of Ω 0 , the backward integration maps the points of Γ 0 along the trajectories of the system into a new surface Γ 1 which bounds the new estimate Ω 1 .The sequence {Ω k } should satisfies the following properties: 1. Ω k ⊂ Ω k+1 , that is {Ω k } should be a strictly monotonically increasing sequence; 2. Every estimate Ω k should belongs to the true stability region.
In [6] the backward integration is performed by choosing p points, P k j , j = 1, 2, ..., p on the boundary surface Γ k and moving each these points by backward integration along the trajectories with the same step h (that is from the time t k to t k+1 = t k + h); the results are the points P k+1 j which define the boundary of the new estimate.The sequence {Γ k } is proved to converge to the boundary surface of the true stability region.
In [14] the possibility of yielding the successive estimates in analytic form is studied.Starting with an initial parametrically defined surface within the true stability region, it use Euler method to produce a sequence of parametrically defined surfaces which approximates the required boundary.
A constructive methodology was proposed in [3].It starts with a given Lyapunov function and yields a sequence of Lyapunov functions which are then used to estimate the stability region.The sequence is shown to satisfies the conditions 1, 2. The methodology proceeds in three main steps: (A) Determining the critical level value of a given Lyapunov function V ; (B) Estimating the stability region via the function V ; (C) Expanding the current estimate; this step is performed via the following expansion schemes: the function V (x) is replaced by either Steps (B) and (C) are then reapplied iteratively.The first expansion scheme is related to the backward Euler numerical procedure.This idea also appear in our paper [10], experiments 5 and 6, pp.62-64, fig.2.3-2.7.
Loccufier and Noldus [9] recently proposed a new trajectory reversing method, a combination of Lyapunov techniques, trajectory reversing and some topological properties of the stability boundary.This method provides an accurate estimation of the true stability region for a wide classes of high order nonlinear dynamical system.
In this paper the expansion schemes based on Euler method are studied and developed.We try to answer to the following question: What is the asymptotical behavior of successive estimates produced by such expansion schemes?In section 2 the expansion schemes are constructed and motivated.The particular case of second order system and explicit form of Lyapunov function is considered in section 3. A convergence theorem is then proved.Section 4 contains an algorithm of trajectory reversing type and a number of illustrative examples.

EXPANSION SCHEMES
Consider a dynamical system which is described by the following nonlinear autonomous system of differential equations Suppose that f (0) = 0 and that the null solution x(t) ≡ 0 is asymptotically stable.Let Ω be the asymptotic stability region of the origin and let Γ be its boundary.Let Γ 0 be an initial estimate of Γ and suppose that there exists a function V 0 : R n → R, sufficiently smooth, such that V −1 0 (0) = Γ 0 .That is Γ 0 is the boundary of the zero level set of V 0 .Starting from Γ 0 it can obtains a new estimate Γ 1 by the standard trajectory reversing technique: let x be an arbitrary point of Γ 0 and let x(t) be the reversing trajectory of (1) which starts from x. Now, perform a reversing displacement along this trajectory with the steplength h and repeat this procedure for all trajectories starting from points on Γ 0 .This means that Γ 0 is shifted along the trajectories in reversing sense with the steplength h.The new position of Γ 0 will be the next estimate Γ 1 .Thus, we yield a sequence {Γ k } of estimates which approximates the boundary Γ of the true stability region.
The following problem arises: knowing V 0 such that For this last purpose, we consider the following slight modification of the procedure: the displacement are performed not just along the trajectories, but along the tangencies of the trajectories.It results the estimate Γ t 1 , close to Γ 1 .The transformation of Γ 0 into Γ t 1 is given by (2) which is just a step of backward numerical integration via Euler method.Let ., ., .denote the usual inner product and the corresponding (Euclidean) norm on R n respectively.Throughout this paper we will consider that f will satisfies the following two basic conditions: (a) f is F-differentiable on a convex and bounded set These conditions and the boundedness of D 0 ensure that both f and f are bounded on D 0 ; let m, M be these boundaries, that is Suppose now that h < 1/M .Then, using perturbation lemma, it results that I − hf (x) is invertible and Define the function ϕ : Theorem 1.Let x be a given point in D 0 and let X be given by (2).Suppose that h satisfies the condition km , and that S(X, r) ⊂ D 0 , where r ≤ 4mh.Then Clearly, x = x is a solution of the equation F (x) = 0. Perform one step with the Newton method starting from the point x0 = X.It results Now, apply Mysovskii theorem [12] in order to estimate the error, that is the quantity x1 − x .
We have where γ = hk and β = 1/(1 − hm), which are the first two conditions of Mysovskii.Also, and α < 1/2.Further, because ∞ j=0 α 2 j −1 < 2 and h/(1+hM ) ≤ 2h, it results r = η ∞ j=0 α 2 j −1 < 4mh and the condition S(X, r) ⊂ D 0 is also satisfied.Therefore, the Mysovskii theorem can be applied.We have Note that 2km 2 /3 depends, generally, of the magnitude of D 0 and the quality of the approximation (3) depends of the size of x.For instance, consider the function which is the right side of the Van der Pol equation.Let h = 0.01.If x = (1, 0.5) T then ϕ(X) = (0.99999999, 0.50000088) T , which is in accordance with (3), while if x = (5, 4) T then ϕ(X) = (4.99994720,400528037) T .Note also that if the stability region is bounded then D 0 may be chosen to covers this region and we can dispose of h that the approximation be suitable.
Based on theorem 1, various expansion schemes can be obtained.Since, from (4), x ≈ ϕ(X), we have Therefore, it results the following expansion schemes.1.This scheme is just the above recurrence formula, that is 2. For h sufficiently small, the function ϕ may be approximated by ϕ(x) ≈ x + hf (x) and ( 4) becomes 3. If V is a real function defined on R n , sufficiently smooth, and if h is sufficiently small, we can write where V is the gradient of V .The expansion scheme is ( 6) , then the third expansion scheme may be written as V k+1 = T h (V k ).It is remarkable the fact that this operator is defined by the linear part of the Taylor development of V (x + hf (x)).
4. This scheme is just the third scheme for the case of the explicit second order system.In this case, we will written the system as ẋ = f (x, y), ẏ = g(x, y), and we search for the function V : R 2 → R in the explicit form V (x, y) = y − v(x).Note that, if we takes the operator T h also in explicit form, , where the superscript indicate the iteration number, and the sequence {v k } will be generates by t h .This operator results by a simple computation, taking into account that, in this case, In the paper [10] the schemes ( 5) and ( 6) were considered in somewhat different form and some numerical experiments concerning the possibility of the estimation of stability regions by these schemes were made.The second expansion scheme (6) was also considered by Chiang and Thorp [3] who have shown that if V 0 is a Lyapunov function then for a finite number of iteration V k are also Lyapunov functions and the critical level sets of V k are strictly increasing estimates and remain inside of the true stability region.
In the same manner (using certain formulas for reverse displacement, performing one step with the Newton method and using Mysovskii theorem), it can obtain various other expansion schemes.For example, in the case of explicit second order systems, it can obtain the expansion scheme ( 8) where the function ϕ is defined by ϕ

THE ASYMPTOTICAL BEHAVIOR
The successive estimates must be "close" to the boundary of the stability region; this means that we need a suitable topology in the space of surfaces from R n .For example, we can use the "distance" between two surfaces as is defined in [1].
The sequence of the functions V k , given by any expansion schemes (4), ( 5) or (6), generally, do not have a punctual convergence.But the sequence {V −1 k (0)} may be convergent to a surface Γ h which must have the following important property: The limit surface Γ h is invariant to the transformation (2) that is, if x ∈ Γ h then also X ∈ Γ h .Moreover, Γ h approximates arbitrarily well the boundary of the true stability region as h → 0.
This remarkable property will be pointed out for the scheme 3 by the next example; for the scheme (7) we will give a convergence theorem (theorem 2 in this section).First of all we will verify this property for the scheme (8) and for the nonlinear system considered in example 1 (Section 4).
By a simple computation it result that the function v * (x) = a/x, where a = (−1 + 2h + √ 1 + 4h)/4h is a fixed point of the iteration (8), that is Γ h = {(x, y) : xy = a}.The invariance property of Γ h to the transformation (2) can be also verified.Moreover, if h → 0 then a → 1 and Γ h tend to the true stability region of the system (Γ = {(x, y) : xy = 1}.
Then the operator t h has a fixed point v * .
Proof.Using ( 9) and the boundedness of δ, we have So, it results and t h is nonexpansive.Using the fixed point theorem of Browder-Gäbel-Kirk (see, for example, [13, pp.62]), it follows that t h has at least one fixed point v * .
Application.Consider again the system from example 2, that is the right side of the considered system is f (x, y) = −x + 2x 2 y, g(x, y) = −y.Let the interval from theorem 2 be I = (−∞, −ε] ∪ [ε, ∞) and let Y be the set { a x , a < 1, x ∈ I}.It can shown that Y ⊆ L 2 d (I) and that it is a bounded, convex and closed set.Let u, v ∈ Y given by u(x) = a/x, v(x) = b/x.By a simple computation, it results and d = D = 2(a+b−1), which is the condition (9).The condition h ≤ 2d/D 2 , which is also required by theorem 2, involves h ≤ 1/(a+b−1) which is satisfied for any x , and 0 < −2ha 2 +(1+2h)a < 1.This means that t h : Y → Y and all conditions of theorem 2 are satisfied.
It results that the operator t h has a fixed point in Y , v * (x) = 1/x; the graph of this function is just the boundary oh the true stability region of the system.Observe that this curve is invariant with respect to the mapping (9) and that in this particular case it does not depends of h.Remark 2. In fact, the convergence of the sequence {v k } is equivalent to the convergence of the numerical sequence {a k }, given by a k+1 = −2ha 2 k + (1 + 2h)a k , which converges to the value one for all 0 < a 0 < 1.

AN ALGORITHM AND ILLUSTRATIVE EXAMPLES
Based on the expansion schemes (1-5) some algorithms for estimating the boundary of the stability regions can be developed.The simplest are just the recurrence formulas (4)(5)(6)(7)(8).In [10] some of such algorithms are presented and numerical experiment concerning the efficiency and accuracy of the algorithms are given.
The algorithm we present uses the expansion scheme 3, formula (6).Let V , f : R n → R be the function defined by V , f (x) = V (x), f (x) .Suppose that V is a function which satisfies the condition (10) V (0 This means that T h (V )(0) −1 = V (0) −1 , that is the function V is invariant to the transformation T h and the property of previous section, ensures that V (0) −1 will approximates the boundary of the stability region.The computation of V which satisfies the condition (11) is the main idea of the algorithm.We search for a function V as a polynomial of degree p, in several variables: where a = (a 1 , a 2 , ...) are the coefficients of polynomial and x = (x 1 , ..., x n ) are independent variables.Generally, such a polynomial function cannot be a solution of (11), that is,generally, it can't find a such that (11) be satisfied.Therefore, we determine the polynomial V such that the condition (11) will be best satisfied.For example, it can determine a and a set of points X j , j = 1, ..., m as the solution of the following constraint optimization problem min a,X j n j=1 V (a, X j ), f (X j ) 2 , V (X j ) = 0, j = 1, ..., m.
In a concrete implementation, the set of points X j , j = 1, ..., m may be chosen as follows.Let the point X j be of the form X j = (x 0 1,j , ..., x 0 n−1,j , x n,j ), where the components x 0 1,j , ..., x 0 n−1,j are given and x n,j is unknown.The fixed components must be chosen such that the point X j = (x 0 1,j , ..., x 0 n−1,j , 0) belongs to the stability region.Thus, the constraint optimization problem involves as scalar unknowns the components of a and x n,1 , ..., x n,m .
This algorithm has been applied to some examples we have found in the literature; in the following a part of them is presented to illustrate the possibility of estimation the stability region.In each example we have used 20 equidistant points (m = 20) on the horizontal axis, bounded by two given numbers, a and b, inside of the true stability region.The boundary of the true stability regions was drawn by a continue lines, while the computed estimates, by dotted lines.
Note that the boundary of the stability region of this equation have two branches which runs to infinity (the boundary of stability region is Γ = {(x, y) : The polynomial of degree 2 in two variables, for a and, with the computer round of error, V (0, 0) −1 = Γ.Note that the boundary contains a saddle point of coordinate (1,2).The computed polynomials are (a) V (x, y) = 0.1106x 2 y + 0.144xy + 0.2y − 1, (b) V (x, y) = 0.0048x 4 y 2 − 0.0433x 3 y 2 + 0.1051x 2 y + 0.3127xy + 0.1591y − 1, and the estimates are drawn in Fig. 5.

CONCLUSIONS
The results of experiments are encouraging.If the boundary is given just by a polynomial, then our algorithm gives this polynomial (example 1).In other cases, the estimates have a suitable accuracy for a moderate degree of polynomials (Example 2, Fig. 4b and example 3, Fig. 5b).
The estimates do not depend essentially of the particular characteristics of the boundary; the boundaries of considered examples have totally different shapes (branches which run to infinity, closed orbit, saddle point).
The computed estimates are not always inside of the true stability regions (example 2, Fig. 4); this undesirable situation seems to appears if the given points X j are far to the stable point and they pack near the boundary.

Fig. 2 .
Fig. 2. The estimation of the boundary for example 2 by polynomials of degree two.

Fig. 3 .
Fig. 3.The estimation of the boundary for example 2 by polynomials of degree four.

Fig. 4 .
Fig. 4. The estimation of the boundary for example 2 by polynomials of degree two.

Fig. 5 .
Fig. 5.The estimation of the boundary for example 3 by polynomials of degree three and six, respectively.