Rev. Anal. Num´er. Th´eor. Approx., vol. 31 (2002) no. 1, pp. 89–102 ictp.acad.ro/jnaat
ON THE EXPANSION SCHEMES IN TRAJECTORY REVERSING METHOD
S¸TEFAN M ˇARUS¸TER∗
Abstract. 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.
MSC 2000. 65P40.
Keywords. nonlinear dynamical systems, stability regions, trajectory reversing method, expansion schemes.
1. 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 se- quence 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 integra- tion maps the points of Γ0 along the trajectories of the system into a new
∗University of the West, Timi¸soara, Romania, e-mail: [email protected].
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, Pjk, j = 1,2, ..., p on the boundary surface Γk and moving each these points by back- ward integration along the trajectories with the same step h(that is from the time tk to tk+1 = tk+h); the results are the points Pjk+1 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 Lya- punov 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) Deter- mining the critical level value of a given Lyapunov functionV; (B) Estimating the stability region via the function V; (C) Expanding the current estimate;
this step is performed via the following expansion schemes: the functionV(x) is replaced by eitherV(x+df(x)) orV(x+d/2(f(x+df(x)) +f(x)), d >0.
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 me- thod, 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.
2. EXPANSION SCHEMES
Consider a dynamical system which is described by the following nonlinear autonomous system of differential equations
(1) x˙ =f(x), f :D⊆Rn→Rn.
Suppose thatf(0) = 0 and that the null solutionx(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 V0 : Rn → R, sufficiently smooth, such that V0−1(0) = Γ0. That is Γ0 is the boundary of the zero level set ofV0. 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 steplengthh. 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 V0 such that V0−1(0) = Γ0, deter- mine V1 such thatV1−1(0) = Γ1 or at least such that V1−1(0) ≈Γ1. 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 Γt1, close to Γ1.
The transformation of Γ0 into Γt1 is given by
(2) X=x−hf(x),
which is just a step of backward numerical integration via Euler method.
Let h., .i , k.k denote the usual inner product and the corresponding (Eu- clidean) norm onRnrespectively. 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 setD0 ⊂D;
(b)kf0(u)−f0(v)k ≤kku−vk, ∀u, v∈D0.
These conditions and the boundedness ofD0 ensure that bothf andf0 are bounded on D0; let m, M be these boundaries, that is
kf(x)k ≤m, kf0(x)k ≤M, ∀x∈D0.
Suppose now that h < 1/M. Then, using perturbation lemma, it results that I−hf0(x) is invertible and
[I−hf0(x)]−1≤ 1
1−hM, ∀x∈D0. Define the functionϕ:D0→Rn by
ϕ(x) =x−h[I−hf0(x)]−1f(x).
Theorem1. Letxbe a given point inD0 and letXbe given by(2). Suppose that h satisfies the condition
0< h≤min2M1 , 1
M+√ km , and that S(X, r)⊂D0, where r≤4mh.
Then
(3) kx−ϕ(X)k ≤ 8km3 2h3.
Proof. Consider the functionF :S(X, r)→Rn given by F(¯x) = ¯x−hf(¯x)−X.
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
x¯1 = ¯x0−[I−hf0(¯x0)]−1(¯x0−hf(¯x0)−X) =ϕ(X).
Now, apply Mysovskii theorem [12] in order to estimate the error, that is the quantity k¯x1−xk.
We have
kF0(u)−F0(v)k ≤γku−vk, kF0(u)−1k ≤β, ∀u, v∈D0,
where γ = hk and β = 1/(1−hm), which are the first two conditions of Mysovskii. Also,kF0(¯x0)−1F(¯x0)k ≤η=hm/(1−hM), therefore the constant α from Mysovskii theorem is
α= 12γβη= km2 1−hMh ,
andα <1/2. Further, becauseP∞j=0α2j−1<2 andh/(1+hM)≤2h, it results r=ηP∞j=0α2j−1 <4mhand the condition S(X, r)⊂D0 is also satisfied.
Therefore, the Mysovskii theorem can be applied. We have ε1= η(1−αα 2) < 4α3η = k31−hMh < 2kh3 ,
kx¯1−x¯0k=hk[I−hf0(X)]−1f(X)k ≤h1−hMm ≤2mh.
Finally, it obtains
kx¯1−xk ≤¯ ε1k¯x1−x¯0k2 ≤ 2km3 2h3. Note that 2km2/3 depends, generally, of the magnitude ofD0 and the qual- ity of the approximation (3) depends of the size of x. For instance, consider the function
f(x) =
−x2 x1−x2+x21x2
,
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 D0 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
V0(ϕ(X))≈V0(x) = 0, ∀x∈Γ0,
which means thatX ∈Γt1 ⇔V0(ϕ(X))≈0 and we can takesV1 =V0◦ϕ.
Therefore, it results the following expansion schemes.
1. This scheme is just the above recurrence formula, that is
(4) Vk+1(x) =Vk(ϕ(x)).
2. Forh sufficiently small, the functionϕmay be approximated by ϕ(x)≈ x+hf(x) and (4) becomes
(5) Vk+1(x) =Vk(x+hf(x)).
3. If V is a real function defined on Rn, sufficiently smooth, and if h is sufficiently small, we can write
V(x+hf(x))≈V(x) +hhV0(x), f(x)i, whereV0 is the gradient of V. The expansion scheme is (6) Vk+1 =Vk+hhVk0, fi.
Remark1. IfThis the operator defined byTh(V) =V+hhV0, fi, then the third expansion scheme may be written asVk+1 =Th(Vk).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
x˙ =f(x, y), y˙ =g(x, y),
and we search for the function V :R2→Rin the explicit form V(x, y) =y− v(x). Note that, if we takes the operatorThalso in explicit form,Th(V(x, y)) = y −th(v(x)), then Th(k)(V) = y −ht(k)h (v), where the superscript indicate the iteration number, and the sequence {vk} will be generates by th. This operator results by a simple computation, taking into account that, in this case, V0(x, y) = (−v(x),1)T. It obtains th(v) = v−hF(v) where F(v(x)) = g(x, v(x))−v0(x)f(x, v(x)).The expansion scheme is
(7) vk+1 =vk−hF(vk).
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 V0 is a Lyapunov function then for a finite number of iteration Vk are also Lyapunov functions and the critical level sets of Vk are strictly increasing estimates and remain inside of the true stability region.
In the same manner (using certain formulas for reverse displacement, per- forming 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) vk+1 =vk◦ϕ−hg(ϕ, vk◦ϕ),
where the function ϕis defined byϕ(x) =x+h[1 +hf0(x, v(x))]f(x, v(x)).
3. 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 Rn. For example, we can use the ”distance” between two surfaces as is defined in [1].
The sequence of the functionsVk, given by any expansion schemes (4), (5) or (6), generally, do not have a punctual convergence. But the sequence{Vk−1(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, ifx∈Γ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}.
An example for the scheme 3. Consider again the system from example (2). The first five terms of the sequenceVk given by the expansion scheme 3,
starting withV0(x, y) =x2+y2−0.25 andh= 0.2, are V1(x, y) =0.6x2+ 0.6y2+ 0.8x3y−0.25,
V2(x, y) =0.36x2+ 0.36y2+ 0.64x3y+ 0.96x4y2−0.25,
V3(x, y) =0.216x2+ 0.216y2+ 0.46x3y+ 0.567x4y2+ 1.536x5y3−0.25, V4(x, y) =0.1296x2+ 0.1296y2+ 0.256x3y+ 0.384x4y2+ 3.072x6y4−0.25, V5(x, y) =0.07776x2+ 0.07776y2+ 0.15488x3y+ 0.2304x4y2−3.072x6y4
+ 0.6144x5y3+ 7.3728x7y5−0.25
It seems that this sequence of functions does not have a punctual convergen- ce; indeed, for instance, ifx=y= 1 then the corresponding numerical sequen- ce is 1.75, 1.75, 2.07, 2.71,3.721, 5.206, . . .; ifx=y= 2 then the correspond- ing numerical sequence is 7.75,17.35, 74.31, 438.214, 3175, 27230, . . .But the sequence of curves V0−1(0), V1−1(0), V2−1(0), V3−1(0), V4−1(0), V5−1(0), . . . seems to converge to a limit curve which approximates the boundary of stability re- gion. In the fig. 1 the initial curve and the second, the forth and the fifth curves are drawn.
Fig. 1. Expansion scheme 3.
Convergence analysis for the scheme3. LetL2(I) be the Hilbert space of the square summable real functions on the interval I, endowed with the usual inner product and corresponding norm and letL2d(I) the subset of L2(I) consisting of derivable functions. LetY be a bounded, convex and closed subset ofL2d(I).
Theorem 2. Suppose that the the operator th mapsY into itself,th :Y → Y, and that the function F, defined in scheme 4,satisfies the condition (9) F(u)−F(v) =δ(u−v), ∀u, v∈Y,
where δ is a real function such that 0 < d ≤ δ(x) ≤ D. Also, suppose that h≤1/2.
Then the operator th has a fixed point v∗.
Proof. Using (9) and the boundedness ofδ, we have
hF(u)−F(v), u−vi=hδ(x)(u−v), u−vi ≥dku−vk2, kF(u)−F(v)k2=kδ(x)(u−v)k2 ≤D2ku−vk2. Ifh <2d/D2 then
2hF(u)−F(v), u−vi> hkF(u)−F(v)k2. So, it results
kth(u)−th(v)k2=ku−v−h(F(u)−F(v)k2
=|u−vk2−2hhF(u)−F(v), u−vi+h2kF(u)−F(v)k2
<ku−vk2,
andth is nonexpansive. Using the fixed point theorem of Browder-G¨abel-Kirk (see, for example, [13, pp. 62]), it follows that th 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 + 2x2y, g(x, y) = −y. Let the interval from theorem 2 be I = (−∞,−ε]∪[ε,∞) and let Y be the set {ax, a < 1, x∈ I}. It can shown that Y ⊆L2d(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
F(u)−F(v) = 2(a+b−1)(u−v),
andd=D= 2(a+b−1),which is the condition (9). The conditionh≤2d/D2, which is also required by theorem 2, involvesh≤1/(a+b−1) which is satisfied for anyu, v because h≤1/2. Finally, if v∈Y then th(v)∈Y, because
th(v) =v−hF(v) = −2ha2+(1+2h)ax ,
and 0<−2ha2+(1+2h)a <1. This means thatth :Y →Y and all conditions of theorem 2 are satisfied.
It results that the operator th 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 {vk} is equivalent to the convergence of the numerical sequence {ak}, given by ak+1 = −2ha2k+ (1 + 2h)ak, which converges to the value one for all 0< a0 <1.
4. 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-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).
LethV0, fi:Rn→R be the function defined by hV0, fi(x) =hV0(x), f(x)i.
Suppose that V is a function which satisfies the condition (10) V(0)−1 =hV0, fi(0)−1.
This means that Th(V)(0)−1 = V(0)−1, that is the function V is invariant to the transformation Th and the property of previous section, ensures that V(0)−1 will approximates the boundary of the stability region. The computa- tion ofV 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:
V(a, x) = X
α1+...+αi n≤p
ai xα11...xαnn
where a = (a1, a2, ...) are the coefficients of polynomial and x = (x1, ..., xn) are independent variables. Generally, such a polynomial function cannot be a solution of (11), that is,generally, it can’t finda such that (11) be satisfied.
Therefore, we determine the polynomial V such that the condition (11) will be best satisfied. For example, it can determineaand a set of pointsXj, j= 1, ..., mas the solution of the following constraint optimization problem
min
a,Xj
n
X
j=1
hV0(a, Xj), f(Xj)i2, V(Xj) = 0, j = 1, ..., m.
In a concrete implementation, the set of points Xj, j = 1, ..., m may be chosen as follows. Let the point Xj be of the formXj = (x01,j, ..., x0n−1,j, xn,j), where the components x01,j, ..., x0n−1,j are given and xn,j is unknown. The fixed components must be chosen such that the pointXj = (x01,j, ..., x0n−1,j,0) belongs to the stability region. Thus, the constraint optimization problem involves as scalar unknowns the components ofa andxn,1, ..., xn,m.
Algorithm 3.Step 0. (Initializations) a = (a01, a02,. . .), xi,j = x0i,j, i = 1, n, j= 1, m;
Step 1. Computexn,j, j= 1, ..., mfrom conditions:
V(a, Xj) = 0, whereXj = (x01,j, ..., x0n−1,j, xn,j);
Step 2. Computea= (a11, a22, ...) such that X
j
hV0(a, Xj), f(Xj)i2 =minimum;
Step 3. Go to Step1.
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 pos- sibility 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,aandb, 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.
Example 1. This is a famous example studied in [7], [4], [3], [11]
x˙ =−x+ 2x2y, y˙=−y.
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) : xy = 1})
The polynomial of degree 2 in two variables, for a=−4, b= 4 is
V(x, y) =−5·10−5x2+ 0.999999xy+ 8·10−9y2−1·10−7x+ 8·10−9y−1, and, with the computer round of error,V(0,0)−1= Γ.
Example 2. The Van der Pol equations (studied in many papers) x˙ =−y,
y˙ =x−y+x2y.
The boundary of stability region is an unstable closed orbit. In the first two experiment (Fig. 2), two polynomials of degree two were computed for a=−1, b= 1 anda=−1.3, b= 1.3 respectively. The polynomials are
(a) V(x, y) = 0.918x20.731xy+ 0.385y2−1,
(b) V(x, y) = 0.5367x2−0.3677xy+ 0.226y2−1.
The next two experiments presents two polynomials of degree four computed for a = −1, b = 1 and a = −1.95, b = 1.95 respectively (Fig. 3). The polynomials are
(a)V(x, y) =−0.0263x4+ 0.121x3y+ 0.006x2y2+ 0.002xy3+ 0.0038y4 + 0.4254x2−0.5375xy+ 0.2339y2−1.,
(b)V(x, y) =−0.019x4+ 0.109x3y+ 0.006x2y2+ 0.0022xy3+ 0.0026y4 + 0.322x2−0.45xy+ 0.195y2−1.
In the Fig. 4 is drawn the graph of a polynomial of degree two for a =
−1.95, b= 1.95.
Fig. 2. The estimation of the boundary for example 2 by polynomials of degree two.
Fig. 3. The estimation of the boundary for example 2 by polynomials of degree four.
Example 3. This is an example studied in [2], [5]
x˙ =−2x+xy, y˙ =−y+xy.
Note that the boundary contains a saddle point of coordinate (1,2). The computed polynomials are
(a) V(x, y) = 0.1106x2y+ 0.144xy+ 0.2y−1,
(b) V(x, y) = 0.0048x4y2−0.0433x3y2+ 0.1051x2y+ 0.3127xy+ 0.1591y−1,
and the estimates are drawn in Fig. 5.
Fig. 4. The estimation of the boundary for example 2 by polynomials of degree two.
Fig. 5. The estimation of the boundary for example 3 by polynomials of degree three and six, respectively.
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 pointsXj are far to the stable point and they pack near the boundary.
REFERENCES
[1] Bogacki, P., Weistein, S. and Xu, Y., Distances between oriented curves in geo- metric modeling, Adv. Comp. Math.,7, pp. 593–621, 1997.
[2] Chiang, H., Hirsch, M.W.andWu, F.F., Stability regions of nonlinear autonomous systems, IEEE Trans. Aut. Control,33, pp. 16–27, 1988.
[3] Chiang, H. and Thorp, J. S., Stability regions of nonlinear dynamical systems: a costructive methodology, IEEE Trans. Aut. Control,34, pp. 1229–1241, 1989.
[4] Genesio, R., Tartaglia, M.,andVicino, A.,On the estimation of asymptotic stabil- ity regions: state of the art and new proposals, IEEE, Trans. on Aut. Control,AC-30, no. 8, pp. 747–755, 1985.
[5] Genesio, R. andVicino, A., New techniques for constructing asimptotic stability re- gions for nonlinear systems, IEEE Trans. Circuits Syst.,CAS-31, pp. 574–581, 1984.
[6] Guttalu, R. andFlashner, H., A numerical method for computing domains of at- traction for dynamical systems, Intern. J. Numer. Meth. Engrg.,26, pp. 875–890, 1988.
[7] Hahn, W., Theory and Applications of Lyapunov’s Direct Method, Prentice Hall, En- glewood Cliffs, NJ, 1963.
[8] Hsu, C. S., Yee and Cheng, H. C., Determination of global regions of asymptotic stability for difference dynamical systems, J. Appl. Mech.,44, pp. 147–153, 1977.
[9] Loccufier, M. and Noldus, E., A new trajectory reversing method for estimating stability regions of autonomous nonlinear systems, Nonlinear Dynamics, 21, pp. 265–
288, 2000.
[10] Mˇarus¸ter, S¸t., Experiments on the regions of asymptotic stability, An. Univ.
Timisoara, ser. Sti. Math.,XXVI, fasc. 3, pp. 53–66, 1988.
[11] Michel, A. N., Sarabudla, N. R. andMiller, R. K.,Stability analysis of complex dynamical systems some computational methods, Cir. Syst. Sign. Process, 1, pp. 171–
202, 1982.
[12] Ortega, J. M.andRheinboldt, W.C.,Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York, 1970.
[13] Rus, I.,Principles of the Fixed Point Theory, Ed. Dacia, Cluj-Napoca, 1980 (in Roma- nian).
[14] Stacey, A. J.andStonier, R. J.,Analitic estimatesfor the boundary of the region of asymptotic attraction, Dynam. Control,8, no. 2, pp. 177-189, 1998.
Received by the editors: January 18, 2002.
