Ngoc Anh Nguyen
email: ngocanh.nguyen.rs@gmail.com

Sorin Olaru
email: sorin.olaru@centralesupelec.fr

A family of piecewise affine control Lyapunov functions

Keywords: Control Lyapunov function, Robust control, Convex lifting

This paper presents a novel method to construct a family of piecewise affine control Lyapunov functions. Unlike most of existing methods which require the contractivity of their domain of definition, the proposed control Lyapunov functions are defined over a so-called N -step controllable set, which is known not to be contractive. Accordingly, a robust control design procedure is presented which only requires solving a linear programming problem at each sampling time. The construction is finally illustrated via a numerical example.

Introduction

As a fundamental concept in control theory [START_REF] Lyapunov | Problème général de la stabilité du mouvement[END_REF], Lyapunov stability has been applied in intensive studies related to stability analysis as well as control design. For the design purpose, control Lyapunov functions are usually employed to synthesize controllers guaranteeing closed-loop stability in the sense of Lyapunov, see among the others [START_REF] Zubov | Methods of AM Lyapunov and their application[END_REF]; [START_REF] Khalil | Nonlinear systems[END_REF]. Such control Lyapunov functions are usually chosen a priori with special structural properties. More clearly, in the case of linear optimal control, suitable quadratic objective functions represent control Lyapunov candidates, see e.g. [START_REF] Anderson | Optimal control: linear quadratic methods[END_REF]; [START_REF] Chmielewski | On constrained infinite-time linear quadratic optimal control[END_REF]; [START_REF] Daafouz | Parameter dependent Lyapunov functions for discrete time systems with time varying parametric uncertainties[END_REF]. Moreover, model predictive control (MPC) usually employs finite/infinite horizon quadratic cost functions as control Lyapunov candidates, see for instance [START_REF] Kothare | Robust constrained model predictive control using linear matrix inequalities[END_REF]; [START_REF] Cuzzola | An improved approach for constrained robust model predictive control[END_REF]; [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF]. Extensive studies about control Lyapunov functions for nonlinear systems have been found in the literature, see among the others [START_REF] Primbs | Nonlinear optimal control: A control Lyapunov function and receding horizon perspective[END_REF]. In case the underlying system is subject to constraints, such a control Lyapunov function should be determined such that the recursive feasibility is ensured. This problem is closely related to the determination of the domain of attraction.

Piecewise linear control Lyapunov functions date back to the studies in [START_REF] Gutman | An algorithm to find maximal state constraint sets for discrete-time linear dynamical systems with bounded controls and states[END_REF]; [START_REF] Nguyen | Interpolation of polytopic control Lyapunov functions for discrete-time linear systems[END_REF] for the nominal case, and are subsequently extended for the robust case to cope with additive disturbances and/or polytopic uncertainty in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF]; [START_REF] Rakovic | Parameterized robust control invariant sets for linear systems: Theoretical advances and computational remarks[END_REF]; [START_REF] Nguyen | Implicit improved vertex control for uncertain, timevarying linear discrete-time systems with state and control constraints[END_REF], leading to simple design formulations as linear programming problems. However, these studies require that such control Lyapunov functions be defined over contractive sets to guarantee its strict decrease and the recursive feasibility.

This paper aims to present the construction of a more general family of control Lyapunov candidates in the context of constrained control, namely piecewise affine functions. These candidates are defined over a so-called N -step controllable set for a given positive integer N, which is obtained from an increasing sequence of N polytopes, and is known to be not necessarily (one step) contractive. Accordingly, we prove that the conditions of a Lyapunov function (the positivity and the strict decrease) are satisfied within this N -step controllable set with a suitable robust control algorithm. Note that the complexity of the proposed control Lyapunov functions increases as N becomes larger, leading to a more complex control algorithm over the ones using contractive set in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF]; Nguyen et al. (2017c), since the number of constraints in the proposed algorithm is larger than the ones in two latter references. However, since this method only requires solving a linear program at each sampling instant, these results can be used for constrained control systems with fast dynamics, e.g. vibration attenuation system, c.f. Gulan et al. (2017b,a).

Generalities and basic notions

Throughout the paper, R, R + , N, N >0 denote the field of real numbers, the set of nonnegative real numbers, the set of nonnegative integers and the positive integer set, respectively. The following index set is also defined, for ease of presentation, with respect to a given N ∈ N >0 : I N = {1, 2, . . . , N }. A polyhedron is defined as the intersection of finitely many closed halfspaces. A polytope is defined as a bounded polyhedron. Also, V(P ) is understood as the set of vertices of polytope P. The distance from a point x ∈ R d to a set S ⊂ R d denoted by ρ S (x) is defined as ρ S (x) := min y∈S (x -y) T (x -y).

Given two sets S 1 , S 2 ⊂ R d , we define the set S 1 \S 2 :=

x ∈ R d : x ∈ S 1 , x / ∈ S 2 . Also, the Minkowski sum of these two sets, denoted by S 1 ⊕ S 2 , is defined as follows:

S 1 ⊕ S 2 = {x 1 + x 2 : x 1 ∈ S 1 , x 2 ∈ S 2 } .
We use ∂S to denote the boundary of a compact set S. Also, int(S) represents the interior of a full-dimensional set S and conv(S) denotes its convex hull. A function α(•) : R + → R + is said to be of class K, if it is continuous, strictly increasing and α(0) = 0.

Problem settings

In this paper, we consider a linear time-varying system:

x k+1 = A(k)x k + B(k)u k + w k , (1)
where x k , u k , w k denote the state, control variable and additive disturbance at time k. These variables satisfy

x k ∈ X, u k ∈ U, w k ∈ W, (2)
where the constraint sets X ⊂ R dx , U ⊂ R du , W ⊂ R dx are assumed to be polytopes, containing the origin in their interior, with given d x , d u ∈ N >0 . This assumption ensures that the origin as the equilibrium point satisfies the above constraints. In case the origin is not the equilibrium point, the system can be translated into the frame of the equilibrium point and the construction can be easily adapted. This assumption will be used later in (4) to build a control Lyapunov function which is only equal to 0 at the equilibrium point. Also, the state-space matrices A(k), B(k) are assumed to belong to a given polytope, defined as below:

[A(k) B(k)] ∈ Ψ := conv {[A 1 B 1 ] , . . . , [A L B L ]} . (3)
This paper aims to construct a new family of control Lyapunov functions, also referred to as convex liftings in the present framework. In particular, the control Lapunov functions presented in this paper are more general than the piecewise linear family proposed in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF]; [START_REF] Rakovic | Parameterized robust control invariant sets for linear systems: Theoretical advances and computational remarks[END_REF], since besides their convex, piecewise affine properties, they are defined over the N -step controllable set, known to be non-necessarily contractive.

Construction of control Lyapunov functions

Before describing the main result, we need to recall some important ingredients which are instrumental for the proposed construction of control Lyapunov functions.

Positive invariance concept has been investigated in many studies [START_REF] Aubin | Differential inclusions: setvalued maps and viability theory[END_REF]; Bitsoris (1988b); [START_REF] Bitsoris | Constrained regulation of linear systems[END_REF]; Bitsoris (1988a); [START_REF] Blanchini | Set-theoretic methods in control[END_REF] and used in different control strategy designs. In case the system is affected by disturbances, the robust positive invariance concept is of use instead.

Definition 4.1 Given an admissible linear control law u k = Kx k ∈ U, a set Ω ⊆ X is called robust positively invariant with respect to system (1) subject to (2) and (3) iff

(A(k) + B(k)K)Ω ⊕ W ⊆ Ω, ∀ [A(k) B(k)] ∈ Ψ.
In order to compute such a set Ω, one should determine a local controller u = Kx ∈ U, which can cope with both the model uncertainty (3) and additive disturbances in W. Such a gain K can be obtained by different methods, see among the others [START_REF] Kothare | Robust constrained model predictive control using linear matrix inequalities[END_REF]; [START_REF] Daafouz | Parameter dependent Lyapunov functions for discrete time systems with time varying parametric uncertainties[END_REF].

According to such a local controller, one can use existing algorithms to compute a robust positively invariant set, e.g., in [START_REF] Gilbert | Linear systems with state and control constraints: The theory and application of maximal output admissible sets[END_REF]; [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF]; [START_REF] Nguyen | Implicit improved vertex control for uncertain, timevarying linear discrete-time systems with state and control constraints[END_REF].

The definition of another important ingredient N -step controllable set is recalled below.

Definition 4.2 Consider system (1) subject to model uncertainty (3) and constraints (2). Given a robust positively invariant set Ω and N ∈ N >0 , a set denoted by K N (Ω) ⊆ X is called the N -step controllable set if any point belonging to this set can reach Ω in N steps in the presence of suitable admissible controller, while staying inside X despite any disturbances in W and model uncertainties in Ψ. It is mathematically defined below for all i ∈ I N :

K 0 (Ω) = Ω, K i (Ω) = x k ∈ X : ∃u k ∈ U s.t. ∀ [A(k) B(k)] ∈ Ψ, (A(k)x k + B(k)u k ) ⊕ W ⊆ K i-1 (Ω) .
To determine K N (Ω), the computation of the 1-step controllable set should be performed, i.e. K i (Ω) should be computed according to K i-1 (Ω). Indeed, if one defines an intermediate set

S := x T u T T ∈ R dx+du : x ∈ X, u ∈ U, (A j x + B j u) ⊕ W ⊆ K i-1 (Ω), ∀j ∈ I L ,
then K i (Ω) can be determined as the orthogonal projection of the set S defined above onto the space of x. Similar computation is repeated until i = N to obtain K N (Ω).

The interested reader is referred to Section 2.6 in [START_REF] Kerrigan | Robust constraint satisfaction: Invari-ant sets and predictive control[END_REF] for further detail. One can easily see that if Ω is empty, then so are K i (Ω) for i ∈ N. Consequently, the proposed method cannot apply in this case. Also, if Ω = ∅ is not of full-dimension, then K i (Ω) might not be of fulldimension either. To illustrate this point, we consider the following simple system:

x k+1 = 1 0 0 α k x k + 1 α k u k ,
where uncertainty

α k ∈ [-1, 1] and control variable u k ∈ [-1, 1]. If Ω = {0}, then one can easily compute K i (Ω) = [y 1 y 2 ] T ∈ R 2 : y 1 = y 2 ∈ [-1, 1] for all i ∈ N >0 .
Although the proposed method can still apply in this case, we exclusively consider the case as presented in Assumption 1 to ensure that K i (Ω) for i ∈ N are of full-dimension.

Assumption 1 Ω is a full-dimensional polytope in R dx .
Note that 0 ∈ int(Ω) since the origin is assumed to be the equilibrium point and a full-dimensional set Ω is robust positively invariant. Furthermore, since Ω satisfies Assumption 1 and X, U, W are polytopes, then K i (Ω) for any finite i ∈ N is also a full-dimensional polytope. Therefore, the existence of a full-dimensional K N (Ω) depends on the existence of a full-dimensional Ω, since they fulfill the following property.

Lemma 4.1 Given a robust positively invariant set Ω satisfying Assumption 1, then K i-1 (Ω) ⊆ K i (Ω) for all i ∈ N.

Clearly, the sequence {K i (Ω)} ∞ i=0 is increasing and bounded above by X, accordingly the limit exists. Note that if the limit of this sequence is finitely determined, there exists

N * ∈ N >0 such that K N * -1 (Ω) ⊂ K N * (Ω) = K N * +1 (Ω).
In this case, any positive integer N < N * is suitable for the proposed construction to avoid

V(K N +1 (Ω))\K N (Ω) = ∅. Otherwise, if the limit of {K i (Ω)} ∞ i=0
is not finitely determined, this end may not be a polytope. In this case, one can always ensure for any N < +∞ that K N (Ω) is a polytope and V(K N +1 (Ω))\K N (Ω) = ∅. As a consequence, any positive integer N can be used in the proposed construction.

Before presenting the main results of the paper, a parametric linear programming (pLP) problem is recalled in the sequel:

max x c T x s.t. Hx ≤ Gλ + b,
where x denotes the decision variable, λ denotes the parameter and H, G, b, c denote matrices of suitable dimensions. Let Λ be the set of λ such that the above pLP problem has a finite, optimal solution for each λ ∈ Λ and no optimal solution for λ / ∈ Λ. Note that the set of constraints Hx ≤ Gλ+b can be easily transformed into equality constraints by introducing an auxiliary variable y of suitable dimension. Then the above pLP problem can be equivalently written as below:

z * (λ) = max x,y c T 0 T x y s.t. [H I] x y = Gλ+b, y ≥ 0,
where I denotes a suitable identity matrix. In the above form, [START_REF] Gal | Postoptimal analyses, parametric programming and related topics[END_REF] shows important properties of z * (λ) via Theorems IV-3 and IV-4 therein, they are recalled below for completeness, while their proof is referred to this reference for more detail.

Theorem 4.1 z * (λ) is a concave function over Λ.

Theorem 4.2 z * (λ) is continuous over Λ.

Given a robust positively invariant set Ω satisfying Assumption 1 and a constant h 0 > 0, we define the following:

V 0 := v T h 0 T : v ∈ V(Ω) ∪ {0} ⊂ R dx+1 , Π 0 := conv( V 0 ), 0 (x) := arg min z z s.t. x T z T ∈ Π 0 .
(4)

x z 0 Π 0 -1 1 1 0 (x) 0 Fig. 1. Illustration for the construction of 0(x) in (4).
Note that a Lyapunov function is equal to 0 only at the origin, therefore the origin is given a height equal to 0, while the vertices of Ω are given a height equal to h 0 > 0 in (4). Subsequently, the augmented set Π 0 is computed as the convex hull of these augmented points. Finally, the lower boundary 0 (x) of Π 0 is computed, since this function represents a control Lyapunov candidate over Ω. This observation will be formally proven in Proposition 5.2. An illustration of construction ( 4) is presented in Fig. 1, in which Ω is defined as interval [-1, 1] , h 0 = 1, and function 0 (x) is represented by the blue segments. Important properties of function 0 (x), defined in (4), are presented in the following lemma.

Lemma 4.2 Function 0 (x) : Ω → R defined in (4)

(1) is a convex, continuous, piecewise affine function;

(2) satisfies 0 ≤ 0 (x) ≤ h 0 for all x ∈ Ω;

(3) satisfies 0 (x) = 0 only if x = 0;

(4) satisfies 0 (x) = h 0 for all x ∈ V(Ω);

(5) satisfies 0 (x) = h 0 for all x ∈ ∂Ω;

(6) satisfies 0 (x) < h 0 for all x ∈ int(Ω).

For reading ease, the proof is referred to Subsection 8.1.

Remark 4.1 We remark that one can choose any value of constant h 0 > 0 to construct 0 (x) as in (4) without affecting its properties shown in Lemma 4.2.

We now present the construction of a control Lyapunov function defined over the N -step controllable set K N (Ω).

For ease of presentation, let N (x) denote such a control Lyapunov function. Function N (x) is expected to be convex, continuous, piecewise affine. Also, the value of N (x) at each vertex v ∈ V(K i (Ω))\K i-1 (Ω) should be strictly larger than the maximal value of N (x) over K i-1 (Ω). This requirement will be used later (more precisely ( 13)-( 15) in Proposition 5.1) to prove the strict decrease of N (x) along the trajectories outside Ω. For simplicity, one can choose the same height

h i for all v ∈ V(K i (Ω))\K i-1 (Ω) such that N (v) = h i for all v ∈ V(K i (Ω))\K i-1 (Ω).
Further, function N (x) should satisfy N (x) = 0 (x) for all x ∈ Ω. In order to find N (x), we construct step by step intermediate functions i (x) for i ∈ I N , defined over

K i (Ω), such that these functions are convex, i (v) = h i for all v ∈ V(K i (Ω))\K i-1 (Ω), and i (x) = i-1 (x) for all x ∈ K i-1 (Ω).
As proven in Lemma 4.2, 0 (x) defined in (4), is a convex, continuous and piecewise affine function. For ease of presentation, let

X (0) j j∈I M (0)
denote the polytopic partition of Ω associated with 0 (x) (the definition of a polytopic partition is referred to Nguyen et al. ( 2017b)) and we denote 0 (x) as follows:

0 (x) = f (0) j T x + g (0) j for x ∈ X (0) j , (5)
where

f (0) j ∈ R dx , g (0) j ∈ R, j ∈ I M (0)
for a suitable M (0) ∈ N >0 . One can further prove that M (0) is equal to the number of facets of Ω. Moreover, g (0) j can be proven to be 0 for all j ∈ I M (0) . Lemma 4.3 Let function 0 (x) be defined in (4) and denoted as in (5). Then,

0 (x) = f (0) j T x for all x ∈ X (0) j .
The proof is referred to Subsection 8.2. This lemma shows that all the regions X (0) j for j ∈ I M (0) share a common vertex as the origin. For ease of presentation, define function:

0 (x) := max j∈I M (0) f (0) j T x + g (0) j for all x ∈ R dx . (6)
In a general manner, the ingredients i (x), i (x), h i , V i , Π i for i ∈ I N are defined step-by-step (recursively) in the sequel. Let function i-1 (x), known as a piecewise affine function over K i-1 (Ω), be given in the following form:

i-1 (x) = f (i-1) j T x + g (i-1) j for x ∈ X (i-1) j
, where

X (i-1) j j∈I M (i-1)
denotes the polytopic partition of

K i-1 (Ω) associated with i-1 (x) and M (i-1) ∈ N >0 . As i-1 (x) is restrictively defined over K i-1 (Ω), we define the following function i-1 (x) := max j∈I M (i-1) f (i-1) j T x + g (i-1) j , ∀x ∈ R dx . (7)
With respect to a given scalar constant > 0, the function i (x) over K i (Ω) is constructed as follows:

i (x) := arg min z z s.t. x T z T ∈ Π i Π i := conv( V i ) (8) V i := v T h i T : v ∈ V(K i (Ω)) ∪ V i-1 h i := min h h s.t. i-1 (v) + ≤ h, ∀v ∈ V(K i (Ω))\K i-1 (Ω).
The constructions ( 7) and ( 8) are iteratively repeated until i = N. Recall that i (x) at each vertex v ∈ V(K i (Ω))\K i-1 (Ω) is required to be strictly larger than the maximal value of i-1 (x) over K i-1 (Ω), this requirement is enforced by the constraints

h i ≥ i-1 (v) + for all v ∈ V(K i (Ω))\K i-1 (Ω),
with a suitable > 0. The proof for

h i > max x∈Ki-1(Ω)
i-1 (x) will be shown later in Lemma 4.6. From the geometric viewpoint, the insertion of a constant > 0 in ( 8) is to ensure that all the vertices v ∈ V(K i (Ω))\K i-1 (Ω) are lifted onto a level above all the hyperplanes composing i-1 (x).

To illustrate the construction (8), consider a simple system T ,

x k+1 = x k + u k + w k , where x k ∈ [-10, 10] , u k ∈ [-0.5, 0.5] , w k ∈ [-0.1, 0.1] . It can easily be observed that a robust positively invariant set is Ω = [-0.2, 0.2] with the local controller u k = -0.5x k . Accordingly, one determines K 1 (Ω) = [-0.6, 0.6] , K 2 (Ω) = [-1, 1] . As shown in
[-0.2 0.1] T , [0.2 0.1] T , [-0.6 0.4] T , [0.6 0.4] T .
Similarly, (8) returns h 2 = 0.8, function 2 (x) is thus composed of the colored segments above the x-axis, the blue ones correspond to z = ±0.5x, while the brown segments are represented by z = ±0.75x -0.05 and the green segments are described by z = ±x -0.2. Some essential properties of functions i (x), i (x) are presented in the sequel. PROOF. As a convex, piecewise affine function over Ω, 0 (x) can also be represented by:

0 (x) = max j∈I M (0) f (0) j T x + g (0) j for all x ∈ Ω. (9)
Therefore, the definition of 0 (x) in (6) satisfies 0 (x) = 0 (x) for any x ∈ Ω. Following similar argument and using Lemma 4.4 lead to i (x) = i (x) for all x ∈ K i (Ω). 2

We now prove that the functions i (x), defined in (8), satisfy the aforementioned requirements. This is formally stated in the following result whose proof is presented in Subsection 8.3 for reading ease. Lemma 4.6 Functions i (x) for i ∈ I N satisfy:

(1) 0 ≤ i (x) ≤ h i for all x ∈ K i (Ω);

(2) h N > . . . > h i > . . . > h 0 > 0;

(3) i (x) = h i for all x ∈ V(K i (Ω))\K i-1 (Ω);

(4) i (x) = i-1 (x) for all x ∈ K i-1 (Ω);

(5) i (x) = h j for x ∈ V(K j (Ω))\K j-1 (Ω), ∀j ∈ I i ;

(6) i (x) > 0 for all x ∈ K i (Ω)\{0};

(7) i (βx) ≤ β i (x) for all x ∈ K i (Ω) and 0 ≤ β ≤ 1.

It can also be proven that such a function i (x) represents a convex lifting for the associated polytopic partition

X (i) j j∈I M (i) of K i (Ω).
The interested reader is referred to [START_REF] Nguyen | Explicit robust constrained control for linear systems: analysis, implementation and design based on optimization[END_REF]; Nguyen et al. (2015a,b) for a precise definition of convex lifting and its existence conditions. For simplicity, we skip this analogy in this paper and stress that those existence conditions are fulfilled in the present context.

Lemma 4.7 Function i (x) : K i (Ω) → R represents a con- vex lifting of the polytopic partition X (i) j j∈I M (i) of K i (Ω) for each i ∈ I N .
PROOF. The proof is referred to Lemma II.8 in Nguyen et al. (2017a) or Lemma 4.7 in Nguyen et al. (2017c). 2 Remark 4.2 From the mathematical point of view, the value of in (8) can be arbitrarily chosen as long as it is positive.

However, the value of h N can be relatively large in some cases, since a subset of the vertices of K N (Ω) and K N -1 (Ω) may be very close to each other, when N becomes larger. This end may present numerical sensitivity in determining N (x), as the precision of convex hull operation is known to be limited, see [START_REF] Avis | A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra[END_REF]. Therefore, in (8) should be chosen small enough to reduce the value of h N and large enough to ensure the precision of the convex hull computation. A value in interval [10 -5 , 1] is basically a good choice.

Control design procedure

The construction and prominent properties of a family of control Lyapunov candidates have been presented in the preceding section. In this section, we need to prove that they satisfy the conditions of control Lyapunov functions. Before proceeding the proof, we recall a definition of a local input-to-state stability (ISS) Lyapunov function from [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]; [START_REF] Khalil | Nonlinear systems[END_REF]:

Definition 5.1 Consider system (1) subject to constraints (2), (3), the feasible region X ⊆ X and a control law

u = κ(x) ∈ U. A continuous function V (x) : X → R + is called a local ISS-Lyapunov function, if the following conditions hold • χ 1 ( x ) ≤ V (x) ≤ χ 2 ( x ) for all x ∈ X , and • V (Ax + Bκ(x) + w) -V (x) ≤ -χ 3 ( x ) + σ( κ(x) )
for all x ∈ X and time-varying

[A B] ∈ Ψ, w ∈ W,
where χ 1 (•), χ 2 (•), χ 3 (•), σ(•) are class K-functions and • denotes any vector norm.

Note that as N (x), constructed in the preceding section, is a convex, continuous, piecewise affine function defined over a bounded set K N (Ω), then there always exist two suitable class K-functions satisfying the first condition in Definition 5.1. This observation is formally stated in the following result.

Lemma 5.1 Given a function N (x) constructed as in (8), then there exist two constant scalars

c 2 ≥ c 1 > 0 such that c 1 x ∞ ≤ N (x) ≤ c 2 x ∞ for all x ∈ K N (Ω). ( 10
)
The proof is presented in Subsection 8.4. Lemma 5.1 shows that N (x) satisfies the first condition in Definition 5.1. Accordingly, in order to prove N (x) to be a control Lyapunov function, one needs to show that there exists a controller u = κ(x) ∈ U for all x ∈ K N (Ω) such that N (x) satisfies the second condition in Definition 5.1. Such a controller can be determined by means of Algorithm 1. The main idea of ( 11) is to minimize the worst case of function N (x k+1 ) despite any [A(k) B(k)] ∈ Ψ and w k ∈ W, by means of minimizing an auxiliary variable γ ≥ 0. Also, constraints

Algorithm 1 Control design procedure Input: N (x), K N (Ω).

Output: optimal control action u * (x k ) at each instant k. 1: Compute N (x k ) 2: Solve the following problem:

γ * (u * k ) T T = arg min γ, u k γ s.t. N (A i x k + B i u k + w) ≤ γ N (x k ), A i x k + B i u k + w ∈ K N (Ω), u k ∈ U, γ ≥ 0, ∀i ∈ I L , ∀w ∈ V(W). ( 11
) 3: u * (x k ) = u * k . 4: k ← k + 1. Return to step 1. A i x k + B i u k + w ∈ K N (Ω)
for all i ∈ I L and w ∈ V(W) are aimed at ensuring that function N (x) is exclusively defined over K N (Ω), i.e. x k+1 stays inside K N (Ω) for all realizations of disturbances and uncertainties. These constraints can be removed if the set Ξ := x ∈ R dx : N (x) ≤ h N and K N (Ω) are identical.

Next, we concentrate on the proof that the convex lifting N (x) and the controller designed in Algorithm 1 satisfy the second condition in Definition 5.1. More precisely, we prove that N (x k ) is strictly decreasing along the trajectories outside Ω. This is formally stated in the following proposition.

Proposition 5.1 Consider function N (x) defined in (8) and the controller designed in Algorithm 1, then for all

[A(k) B(k)] ∈ Ψ and w k ∈ W, any x k ∈ K N (Ω)\Ω satisfies: N (A(k)x k + B(k)u * (x k ) + w k ) < N (x k ). PROOF. Consider any x ∈ K N (Ω), we define i(x) := arg min j∈{0}∪I N j s.t. x ∈ K j (Ω). (12)
It can be seen that the definition of the index i(x), for each x ∈ K N (Ω), is unique. Consider now any point v ∈ j∈I N V(K j (Ω))\Ω; according to definition (12) and claim 5) of Lemma 4.6, it yields

N (v) = h i(v) . (13)
As v ∈ K i(v) (Ω), there exists a controller, denoted by

u(v) ∈ U, such that A(k)v + B(k)u(v) + w k ∈ K i(v)-1 (Ω)
for all [A(k) B(k)] ∈ Ψ and all w k ∈ W according to its definition. Consequently, claim 1) of Lemma 4.6 leads to:

N (A(k)v + B(k)u(v) + w k ) ≤ h i(v)-1 . (14)
Inclusions ( 13), ( 14) and claim 2) of Lemma 4.6 yield:

N (A(k)v + B(k)u(v) + w k ) < N (v), (15)
for all [A(k) B(k)] ∈ Ψ and all w k ∈ W. Otherwise, if v ∈ V(Ω), there exists a control law u(v) ∈ U such that A(k)v + B(k)u(v) + w k ∈ Ω, since Ω is robust positively invariant. Accordingly, we obtain

N (A(k)v + B(k)u(v) + w k ) ≤ h 0 = N (v). (16)
Note that for any point x k ∈ K N (Ω)\Ω, there exists a region X

(N ) j in the polytopic partition

X (N ) j j∈I M (N )
, associated with N (x), such that x k ∈ X (N ) j

. Accordingly, x k can be written in the following form:

x k = v∈V(X (N ) j ) α(v)v, α(v) ≥ 0, v∈V(X (N ) j ) α(v) = 1. As x k ∈ K N (Ω)\Ω, then there exists at least one vertex v ∈ V(X (N ) j
)\Ω such that α(v) > 0. Accordingly, it yields:

N (x k ) = v∈V(X (N ) j ) α(v) N (v) (17a) > v∈V(X (N ) j ) α(v) N (A(k)v + B(k)u(v) + w k ) (17b) ≥ N (A(k)x k + B(k) v∈V(X (N ) j ) α(v)u(v) + w k ) (17c) ≥ N (A(k)x k + B(k)u * (x k ) + w k ). (17d)
Note that inclusion (17a) is obtained because function N (x) is affine over X (N ) j

. Meanwhile, inclusion (17b) is obtained due to inclusions (15) and ( 16). Also, inequality (17c) is induced from the convexity of N (x), as proven in Lemma 4.4. Finally, inclusion (17d) is guaranteed since U is convex and v∈V(X (N ) j ) α(v)u(v) ∈ U. The proof is complete. 2 Proposition 5.2 Consider function N (x) defined in (8) and the controller designed in Algorithm 1, then for all

[A(k) B(k)] ∈ Ψ and w k ∈ W, any point x k ∈ Ω satisfies N (A(k)x k + B(k)u * (x k ) + w k ) ≤ N (x k ) + N (w k ). PROOF. As x k ∈ Ω, there exists a region X (0) j in the poly- topic partition X (0) j j∈I M (0)
of Ω, associated with 0 (x), such that x k ∈ X (0) j . Such a region X (0) j has its vertices as vertices of Ω and the origin. Therefore, x k can be written in the following form:

x k = v∈V(X (0) j ) α(v)v, α(v) ≥ 0, v∈V(X (0) j ) α(v) = 1. Also, since Ω represents a robust positively invariant set, then for any v ∈ V(Ω), there exists a control law u(v) ∈ U such that A(k)v+B(k)u(v)+w k ∈ Ω for all [A(k) B(k)] ∈ Ψ and w k ∈ W.
In other words, in-clusion ( 16) holds at v ∈ V(Ω). Accordingly, we obtain

N (x k ) = v∈V(X (0) j ) α(v) N (v) ≥ v∈V(X (0) j )\{0} α(v) N (A(k)v + B(k)u(v) + w k ). (18)
If we choose u(0) = 0, then inclusion (18) yields

N (x k ) + N (w k ) ≥ N (x k ) + α(0) N (w k ) ≥ v∈V(X (0) j ) α(v) N (A(k)v + B(k)u(v) + w k ). (19)
The convexity of N (x) also leads to

v∈V(X (0) j ) α(v) N (A(k)v + B(k)u(v) + w k ) ≥ N (A(k)x k + B(k) v∈V(X (0) j ) α(v)u(v) + w k ) (20a) ≥ N (A(k)x k + B(k)u * (x k ) + w k ). ( 20b
)
Again, inclusion (20b) follows the convexity of U and v∈V(X

(0) j ) α(v)u(v) ∈ U.
Accordingly, incorporating inclusions ( 19) and (20b) yields

N (x k ) + N (w k ) ≥ N (A(k)x k + B(k)u * (x k ) + w k ), for all [A(k) B(k)] ∈ Ψ and w k ∈ W. 2
Remark 5.1 Note that the choice of Ω is arbitrary as long as its polytopic topology and robust positively invariant properties are guaranteed. The minimal and maximal robust positively invariant sets, e.g. in [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF] are preferable.

The main result related to control design of the paper is formally stated in the following theorem: Theorem 5.3 Given system (1) subject to uncertainty (3) and constraint (2), the controller designed in Algorithm 1 guarantees the recursive feasibility and input-to-state stability.

PROOF. The recursive feasibility is guaranteed by sufficiently large γ > 0 according to Propositions 5.1 and 5.2. Also, N (x), constructed in (8), satisfies the conditions of an ISS Lyapunov function according to Definition 5.1. Therefore, input-to-state stability is proven. More intuitively, Proposition 5.1 shows that {

N (x k )} ∞ k=0 outside Ω is strictly decreasing and bounded in [h 0 , h N ], therefore lim k→∞ N (x k ) = h 0 , consequently lim k→∞ ρ Ω (x k ) = 0. The proof is complete. 2
Remark 5.2 It is worth stressing that neither finite cost function usually employed in model predictive control [START_REF] Mayne | Constrained model predictive control: Stability and optimality[END_REF], nor the control Lyapunov function computed in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF] can guarantee their strict decrease outside a given robust positively invariant set Ω like function N (x) as proven in Proposition 5.1.

Remark 5.3 Although the construction of K N (Ω) is much cheaper than the one of the maximal λ-contractive set, since the latter one is computed mainly relying on iterative procedures, see [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF]; [START_REF] Kerrigan | Robust constraint satisfaction: Invari-ant sets and predictive control[END_REF], the construction of N (x) is however performed at the price of solving multiple parametric linear programming problems, see among the others [START_REF] Grancharova | Explicit nonlinear model predictive control: theory and applications[END_REF]. Also, since the proposed construction relies on the vertex representation, it is thus limited to small-dimensional systems.

Remark 5.4 The number of regions in the polytopic partition associated with N (x) may be relatively large in comparison to the methods in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF]; Nguyen et al. (2017c), because its construction relies on combinations of the vertices of the sets K i (Ω) for i ∈ {0} ∪ I N . This may thus slow down the online computation. However, since this problem is a linear program, it is not expensive.

Numerical example

This section aims to illustrate the above construction via a numerical example. To this end, we consider an angular antenna positioning system, presented in [START_REF] Kothare | Robust constrained model predictive control using linear matrix inequalities[END_REF],

x k+1 = 1 0.1 0 1 -0.1β k x k + 0 0.0787 u k , (21)
where the uncertain parameter β k ranges in [0.1 10] . Also, the state, control variables are subject to the following constraints:

x k ∞ ≤ 1, u k ∞ ≤ 2. (22)
A stabilizing local controller u = [-3.9922 -6.5135] x is chosen to compute Ω as the maximal output admissible set. Accordingly, Ω and the 10-step controllable set are presented in Fig. 3. In this example, we choose h 0 = 0.1 and = 10 -4 . The constructed control Lyapunov function is shown in Fig. 4, where the strict decrease of its value along the closed-loop dynamics outside the robust positively invariant set Ω is also proven with the controller designed by Algorithm 1. The closed-loop stability in the sense of Lyapunov is more clearly illustrated in Fig. 3. Also, the design problem ( 11) is composed of 887 constraints, as N (x) associates with a partition of 402 regions. Note however that the method presented in [START_REF] Blanchini | Ultimate boundedness control for uncertain discrete-time systems via set-induced Lyapunov functions. Automatic Control[END_REF] exploiting the maximal 0.999-contractive set for the same numerical example can formulate a control problem with 40 constraints, while the method in Nguyen et al. (2017c) with the same 0.999-contractive set requires solving a problem consisting of 61 constraints. We remark that Fig. 5 represents the 1-step controllable set K 1 (Ω), which is not contractive since two vertices of Ω, surrounded by two red circles, lie on the boundary of K 1 (Ω). This implies that for some points on the boundary of K 1 (Ω), there does not exist a controller satisfying the given constraints such that their next state stay in the interior of K 1 (Ω). Also, to illustrate more clearly control Lyapunov function 10 (x), we show in Fig. 6 function 1 (x) defined over set K 1 (Ω). Finally, the computation of N (x) is performed in the environment of MPT 3.0 [START_REF] Herceg | Multi-parametric toolbox 3.0[END_REF].

Conclusions

This paper presented a novel method to construct a family of piecewise affine control Lyapunov functions. The most remarkable property of this construction is that it does not require the contractivity of the domain over which these control Lyapunov functions are defined. Accordingly, this construction was shown to lead to a simple robust control design procedure as a linear programming problem. The construction was finally illustrated via a numerical example.

Ω K 1 (Ω) x 1 x 2 Fig. 5. K1(Ω) is not contractive. x 2 x 1 1 (x)
Acknowledgment: The authors would like to thank the associate editor and the reviewers for their helpful comments.

Appendix

Proof of Lemma 4.2

For claim 1), 0 (x) is a piecewise affine function, since it is resulted from a pLP problem. Also, the convexity and continuity of 0 (x) follow as direct consequences of Theorems 4.1 and 4.2. For claim 2), it is observed that any point x T z T ∈ Π 0 can be expressed as a convex combination of the points in V 0 , defined in (4), i.e.,

x T z T = α(0)0 + v∈V(Ω) α(v) v T h 0 T α(0), α(v) ≥ 0, α(0) + v∈V(Ω) α(v) = 1. (23)
According to inclusion (23), claim 2) is easily deduced by the argument: 0 ≤ z = v∈V(Ω) α(v)h 0 ≤ h 0 . Also, z = 0 holds only if α(v) = 0 for all v ∈ V(Ω) and α(0) = 1, leading to claim 3).

For claim 4), one can see that any x ∈ V(Ω) can be expressed as in ( 23) only if α(v) = 0 for v ∈ V(Ω)\{x} and α(v) = 1 for v = x, since x represents a vertex of Ω. In other words, 0 (x) = z = h 0 for x ∈ V(Ω). To prove claim 5), we note that any point x ∈ ∂Ω can be written in form ( 23) only if α(0) = 0. This end leads to 0 (x) = z = h 0 .

Finally, consider any point x ∈ int(Ω), then there exists a point y(x) ∈ ∂Ω and 0 ≤ β < 1 such that x = βy(x). Accordingly, the convexity of 0 (x) yields: 0 (x) = 0 (βy(x)) ≤ β 0 (y(x)) + (1 -β) 0 (0) = βh 0 < h 0 . This end completes the proof of claim 6). 2

8.2 Proof of Lemma 4.3

Consider any point x T z T ∈ Π 0 such that x ∈ Ω\{0},
this augmented point is represented by a convex combination of the points of V 0 as shown in ( 23). Accordingly, one obtains z = (1 -α(0))h 0 . It can be seen that the minimal value 0 (x) of z holds at the maximal value of α(0), denoted by α * (0). Since x ∈ Ω\{0}, then α * (0) < 1. Let α * (v) denote the optimal value of α(v) such that z = 0 (x), we will prove that x/(1 -α * (0)) ∈ ∂Ω. Indeed, one can see that 0 (x) = (1 -α * (0))h 0 . Note also that

v∈V(Ω) α * (v) = 1 -α * (0), leading to x/(1 -α * (0)) = (1 -α * (0)) -1 v∈V(Ω) α * (v)v ∈ Ω.
Therefore, the convexity of 0 (x) yields

h 0 = 0 (x)/(1 -α * (0)) ≤ 0 (x/(1 -α * (0))) ≤ h 0 .
According to Lemma 4.2, x/(1 -α * (0)) ∈ ∂Ω. Consequently, there exists a facet of Ω denoted by F such that x/(1 -α * (0)) ∈ F. Geometrically, the point x/(1 -α * (0)) can be obtained by the intersection between ∂Ω and the line going through 0, x such that x lies between 0 and x/(1 -α * (0)). Also, x T 0 (x) T can be expressed by

x 0 (x) = α * (0)0 + (1 -α * (0)) v∈V(F ) β(v) v 0 (v) β(v) ≥ 0, v∈V(F ) β(v) = 1. (24)
Similarly, for any x ∈ conv ({0} ∪ V(F )) , the point x T 0 (x) T is described as in (24). In other words, conv ({0} ∪ V(F )) represents a region in the polytopic partition X

(0) j j∈I M (0)
of Ω associated with 0 (x). Furthermore, the regions in this partition share a common point 0. Also, 0 (0) = 0 according to claim 3) of Lemma 4.2, it yields g (0) j = 0 for all j ∈ I M (0) . 2 8.3 Proof of Lemma 4.6

First, we will prove claims 1), 2) and 3) simultaneously. Obviously, 0 (x) = h 0 for x ∈ V(Ω), we need to prove that 1

(x) = h 1 for all x ∈ V(K 1 (Ω))\Ω, 0 ≤ 1 (x) ≤ h 1 and h 1 > h 0 . Indeed, any point x ∈ K 1 (Ω) is written as a convex combination of the vertices of K 1 (Ω), i.e. x = v∈V(K1(Ω)) α(v)v, α(v) ≥ 0, v∈V(K1(Ω)) α(v) = 1.
According to the definition of h 1 in (8), we obtain

v∈V(K1(Ω))\Ω α(v)h 1 + v∈V(K1(Ω))∩Ω α(v)h 0 ≥ v∈V(K1(Ω))\Ω α(v) + v∈V(K1(Ω)) α(v) 0 (v) ≥ v∈V(K1(Ω))\Ω α(v) + 0 (x). (25)
Note that the last inequality in inclusion ( 25) is obtained due to the convexity of function 0 (x) over R dx . Since this inclusion holds true for all x ∈ K 1 (Ω), thus if one chooses x ∈ V(Ω) such that there exists at least one v ∈ V(K 1 (Ω))\Ω, α(v) > 0, then it yields h 1 ≥ h 0 + > h 0 . This end leads to 0 ≤ 1 (x) ≤ h 1 according to the definition of Π 1 , 1 (x) in (8) and similar argument as in the proof for claim 2) of Lemma 4.2.

To prove that 1 (x) = h 1 for all x ∈ V(K 1 (Ω))\Ω, it can easily be seen that 1 (x) is a convex function over K 1 (Ω) as proven in Lemma 4.4, then it attains its maximal value h 1 at its vertices. Note however that for the points x ∈

V(K 1 (Ω))∩Ω, 1 (x) ≤ h 0 < h 1 . Therefore, 1 (x) = h 1 for x ∈ V(K 1 (Ω))\Ω.
For the other i ∈ I N , the proof follows similar arguments.

To prove claim 4), first we prove it holds with i = 1. Indeed, consider any point x T z T ∈ Π 1 ; this point can be expressed as a convex combination of the points in V 1 as:

x T z T = β(0)0 + 1 j=0 v∈V(Kj (Ω)) α j (v) v T h j T β(0), α j (v) ≥ 0, β(0) + 1 j=0 v∈V(Kj (Ω)) α j (v) = 1. (26)
The above inclusion yields

z = 1 j=0 v∈V(Kj (Ω)) α j (v)h j ≥ v∈V(Ω) α 0 (v) 0 (v) + v∈V(K1(Ω)) α 1 (v)( 0 (v) + ) ≥ 0 (x) + v∈V(K1(Ω)) α 1 (v) ≥ 0 (x). (27)
The equality in ( 27) holds only if α 1 (v) = 0 for all v ∈ V(K 1 (Ω)), leading to x ∈ Ω. Suppose claim 4) holds true until i -1, we will prove it holds for i. In fact, claim 4) holding true until i -1 leads to j (x) = j-1 (x) for x ∈ K j-1 (Ω), j ∈ I i-1 . Accordingly, consider any point x T z T ∈ Π i , then it can be written in the following form:

x T z T = β(0)0 + i j=0 v∈V(Kj (Ω)) α j (v) v T h j T β(0), α j (v) ≥ 0, β(0) + i j=0 v∈V(Kj (Ω)) α j (v) = 1. (28) Note that each v ∈ V(K j (Ω)) satisfies h j ≥ j (v) = i-1 (v) for j ∈ {0} ∪ I i-1 , therefore inclusion (28) leads to z = i j=0 v∈V(Kj (Ω)) α j (v)h j (29) ≥ i-1 j=0 v∈V(Kj (Ω)) α j (v) i-1 (v) + v∈V(Ki(Ω)) α i (v)h i .
According to Lemma 4.5 and definition (8), we have:

i-1 (v) = i-1 (v) for all v ∈ K i-1 (Ω), i-1 (v) + ≤ h i for all v ∈ V(K i (Ω)). (30)
Recall that proving h i ≥ h i-1 + follows similar arguments as showing h 1 ≥ h 0 + . According to claim 1), we obtain

h i ≥ i-1 (v)+ = i-1 (v)+ for v ∈ V(K i (Ω))∩K i-1 (Ω).
Therefore, inclusions (29) and (30) yield:

z ≥ i-1 j=0 v∈V(Kj (Ω)) α j (v) i-1 (v) + v∈V(Ki(Ω)) α i (v)( i-1 (v) + ). (31)
As i-1 (x) is a convex function, it satisfies:

β(0) i-1 (0) + i-1 j=0 v∈V(Kj (Ω)) α j (v) i-1 (v) + v∈V(Ki(Ω)) α i (v) i-1 (v) ≥ i-1 (x). (32)
Inclusions ( 31) and (32) yield:

z ≥ i-1 (x) + v∈V(Ki(Ω)) α i (v) ≥ i-1 (x). (33)
The equality in (33) holds only if α i (v) = 0 for all v ∈ V(K i (Ω)). In other words, z = i-1 (x) only if x = β(0)0 + i-1 j=0

v∈V(Kj (Ω)) α j (v)v ∈ K i-1 (Ω). More precisely, i (x) = i-1 (x) = i-1 (x) takes place only if x ∈ K i-1 (Ω). Following similar argument, we conclude that claim 4) holds true for all i ∈ I N .

Claim 5) follows as a direct consequence of claims 3) and 4). Indeed, i (x) = j (x) for all x ∈ K j (Ω) and j ∈ {0} ∪ I i . On the other hand, claim 3) shows that j (x) = h j for x ∈ V(K j (Ω))\K j-1 (Ω). Therefore, i (x) = h j for x ∈ V(K j (Ω))\K j-1 (Ω) and j ∈ I i .

Claim 6) can easily be deduced from equation ( 29) that z = β(0)0 + i j=0 v∈V(Kj (Ω)) α j (v)h j ≥ 0, (34) as h j > 0 for all j ∈ {0} ∪ I i , thus the inequality (34) becomes equality only if α j (v) = 0 for all v ∈ i j=0 V(K j (Ω)) and β(0) = 1, leading to x = 0.

Claim 7) is obtained by the convexity of i (x), i.e. i (βx) ≤ β i (x) + (1 -β) i (0) = β i (x). The proof is complete. 2 8.4 Proof of Lemma 5.1

Consider the polytopic partition X ), there exists a c

c (j) 1 v ∞ ≤ f (N ) j T v + g (N ) j
, ∀v ∈ V(X (N ) j

). (36)

Note that any x ∈ X (N ) j can be expressed as follows:

x = v∈V(X (N ) j ) α(v)v, α(v) ≥ 0, v∈V(X (N ) j ) α(v) = 1.
The convexity of function

• ∞ yields c (j) 1 x ∞ ≤ v∈V(X (N ) j ) c (j) 1 α(v) v ∞ ≤ v∈V(X (N ) j ) α(v) f (N ) j T v + g (N ) j = f (N ) j T x + g (N ) j . ( 37
)
Recall that the existence of c 1 > 0, then we can ensure that c 1 x ∞ ≤ N (x) for all x ∈ K N (Ω).

In order to prove the existence of c 2 , consider again region X (N ) j

. First, we observe that f (N ) j T x ≤ f (N ) j 1 x ∞ . Also, since this region is a polytope due to the boundedness of X, then there exists a constant b j > 0 such that g (N ) j ≤ b j x ∞ for all x ∈ X (N ) j .

(38)

Note that if X (N ) j contains the origin, then g 2 > 0, then we ensure that N (x) ≤ c 2 x ∞ for all x ∈ K N (Ω). 2

  Fig. 2, we choose h 0 = 0.1, = 0.1; polytope Π 0 is thus computed by the convex hull of [0 0] T , [-0.2 0.1] T , [0.2 0.1] T and 0 (x) is represented by two blue segments above the interval [-0.2, 0.2] . Subsequently, the construction (8) results in h 1 = 0.4 and the vertices of Π 1 consist of [0 0]

Fig. 2 .

 2 Fig. 2. Illustration for the construction of N (x).

Fig. 3 .

 3 Fig. 3. Maximal output admissible set Ω, its 10-step controllable set K10(Ω) and the closed-loop dynamics.

Fig. 4 .

 4 Fig. 4. The constructed control Lyapunov function with h0 = 0.1, = 10 -4 and its strict decrease along the closed-loop dynamics outside Ω.

Fig. 6 .

 6 Fig. 6. Function 1(x) defined over K1(Ω).

  associated with N (x). Denote also such a function N (x) as follows: N (x) = f partition. Due to the boundedness and positivity of N (x) except at 0, for each v ∈ V(X (N ) j

  Moreover, if v = 0, then inclusion (35) holds true for all c (j) 1 (0) ∈ R, since its two sides are equal to 0. As a consequence, if one chooses c

  therefore any b j > 0 satisfies (38). Otherwise, one can choose b j as follows: b j = min x∈X b j satisfies (38) and b j > 0. We now define c


Carla Maria Pontes

Giancarlo Queiroz Pellegrino
email: giancarloqpellegrino@cefetmg.br

Analysis of quantum phase transition in some different Curie-Weiss models: a unified approach *

Keywords: numbers: 05.70.Fh, 73.20.Mf, 03.65.Sq Keywords: semiclassical limit, quantum phase transition

A unified approach to the analysis of quantum phase transitions in some different Curie-Weiss models is proposed such that they are treated and analyzed under the same general scheme. This approach takes three steps: balancing the quantum Hamiltonian by an appropriate factor, rewriting the Hamiltonian in terms of SU (2) operators only, and obtention of a classical Hamiltonian. SU (2) operators are obtained from creation and annihilation operators as linear combinations in the case of fermions and as an inverse Holstein-Primakoff transformation in the case of bosons. This scheme is successfully applied to Lipkin, pairing, Jaynes-Cummings, bilayer, and Heisenberg models.

Introduction

Quantum phase transition (QPT) has been given much attention and scientific effort during the last years since it proved to be a collective phenomenon observed in very different areas [START_REF] Sachdev | Quantum Phase Transitions[END_REF]. More recently, the concept -originally defined as a property of the ground-state energy -was enlarged to encompass also the manifestations found in excited states (ESQPT) [START_REF] Caprio | Excited state quantum phase transitions in many-body systems[END_REF][START_REF] Cejnar | Impact of quantum phase transitions on excitedlevel dynamics[END_REF]. One of the major interests in this field comes from the possibility of a better understanding of the intricate behavior of many-body systems, observed in very different physical situations, although "direct comparison of treatments using different complicated notations is difficult", as described by Lipkin, Meshkov, and Glick in their famous 1965 paper [START_REF] Lipkin | Validity of many-body approximation methods for a solvable model[END_REF].

We propose here a unified approach to analyze the occurrence of QPT in a set of models coming from different areas: the Lipkin and pairing models from nuclear physics [START_REF] Lipkin | Validity of many-body approximation methods for a solvable model[END_REF][START_REF] Kerman | Accuracy of the superconductivity approximation for pairing forces in nuclei[END_REF], the N -atom Jaynes-Cummings model from quantum optics [START_REF] Dicke | Coherence in spontaneous radiation processes[END_REF][START_REF] Shore | The Jaynes-Cummings Model[END_REF][START_REF] Garraway | The Dicke model in quantum optics: Dicke model revisited[END_REF], the bilayer model from condensed matter [START_REF] Moreira | Entanglement and classical instabilities: Fingerprints of electron-hole-to-exciton phase transition in a simple model[END_REF][START_REF] Figueiredo | Classical analysis of quantum phase transitions in a bilayer model[END_REF], and the Heisenberg model from magnetism [START_REF] Blundell | Magnetism in Condensed Matter[END_REF]. Although they are very different in their physical natures, these models share the common feature that the expected mean values of their observables admit an expansion in powers of 1/N , where N is, generically speaking, the number of constituents of the system; in short, they are Curie-Weiss models [START_REF] Yaffe | Large N -limits as classical mechanics[END_REF]. This property is deeply connected with the possibility of a (semi)classical description based on a proper Hamiltonian function corresponding to the quantum Hamiltonian operator for the model in question. Several quantum-classical connections can then be made, which lead to semiclassical characterization and analysis of QPT in these models. Such connections range from quantum spectra and classical orbits to critical value of a parameter at QPT and the change in (in)stability in phase space.

With this perspective, we propose the following scheme to analyze the above-cited models in a unified approach.

Firstly, the original quantum Hamiltonian is rewritten in terms of SU (2) operators only. This is achieved either by defining pseudo-spin or pseudoangular momentum operators as suitable combinations of the original fermionic creation-annihilation ones, or by performing an (inverse) Holstein-Primakoff transformation [START_REF] Holstein | Field dependence of the intrinsic domain magnetization of a ferromagnet[END_REF] on the original bosonic creation-annihilation operators to produce again SU (2) operators.

Secondly, the quantum Hamiltonian, written as

H = H 0 + gH int , (1)
is balanced by a factor N s [14], where N counts the (constant) number of particles in the system and

s = [H int ] -[H 0 ] . (2)
In this expression [H a ] is the highest power of SU (2) operators in H a . In a final step, a classical Hamiltonian function is obtained from these SU (2)-based quantum Hamiltonians, by following a prescription set by Lieb [START_REF] Lieb | The classical limit of quantum spin systems[END_REF]. With the quantum and classical Hamiltonians, it is possible to obtain, and compare, results associated with QPT both on quantum and classical levels. These steps are summarized as

H a, a † ; b, b † → H (J z ; J ± ) → H balanced (J z ; J ± ) → H (p, q) , (3)
where a and a † are generic fermionic operators, b and b † are bosonic ones, J z and J ± are sets of SU (2) operators, and p and q are pairs of canonical conjugate classical momenta and coordinates. H and H stand for quantum and classical Hamiltonians respectively. Moreover, since all the models we will treat can be put in the form

H = H 0 + g N s H int , (4)
where H 0 sums up the energy of the free constituents and H int takes into account the contribution to the energy given by the interaction between those elements, it will be seen that the classical Hamiltonians can be put in the form

H = p n + λH int (p, q) , (5)
where n is 0, 1 or 2, depending on the model, and λ = λ(g) relates the quantum interaction parameter g to the corresponding classical one λ. It will then be seen that, in these models, QPT presents universal characteristics concerning quantum spectra and the behavior of mean values as functions of the interaction parameter g, and also concerning the structure and behavior of classical orbits with respect to the parameter λ.

That scheme is illustrated in some detail with the Lipkin model and then is applied comparatively to the other models. In each case, a brief description of the model is given, a quantum treatment characterizes QPT in terms of the energy spectra and of a given mean value, both as functions of the interaction parameter g; finally, a classical treatment analyzes the classical phase space in order to characterize QPT in terms of the orbits and (in)stability of the critical points.

In the next section, the models are briefly presented, the above-described scheme is applied to each model, and the results are shown. In section 3, the results are compared and discussed from the perspective of common features eventually present. A short section 4 concludes this article.

Models, treatments, and results

In this section, we take the Lipkin model as a prototype to illustrate our procedure. The other models are treated more briefly within the same scheme. In all cases, sets of SU (2) operators are designed to accomplish for the interactions between the constituents of the model and to satisfy the usual commutation relations

[J z , J ± ] = ±J ± [J + , J -] = 2J z . (6)
Also, we illustrate the quantum results in the same mathematical situation: the operators H are diagonalized within subspaces of the same dimension 2J + 1 with J = 100 and with the interaction parameter g varying in [0, 3]. Quantum mean values are calculated for the ground state. Classical phase spaces are shown for values of the interaction parameter λ before and after the transition. We note here that the results presented below agree with those eventually found in the literature and obtained by other methods, when this is the case.

The Lipkin model

The Lipkin-Meshkov-Glick1 model appears for the first time in Nuclear Physics as a test for the validity of some many-body approximation methods [START_REF] Lipkin | Validity of many-body approximation methods for a solvable model[END_REF]. Since then, it has become a treatable and very useful model in many branches [START_REF] Leyvraz | Large-n Scaling Behavior of the Lipkin-Meshkov-Glick Model[END_REF][START_REF] Rosensteel | Equations of motion for a spectrumgenerating algebra: Lipkin-Meshkov-Glick model[END_REF][START_REF] Morrison | Collective spin systems in dispersive optical cavity QED: Quantum phase transitions and entanglement[END_REF][START_REF] Ribeiro | Exact spectrum of the Lipkin-Meshlov-Glick model in the thermodynamic limit and finite size corrections[END_REF][START_REF] Ping | Decay of Loschmidt echo at a critical point in the Lipkin-Meshkov-Glick model[END_REF]. Particularly interesting is the fact that, being exactly solvable, the Lipkin model can serve as a reliable laboratory to test ideas and methods related to QPT [START_REF] Heiss | The large N behavior of the Lipkin model and exceptional points[END_REF][START_REF] Heiss | On the thermodynamic limit of the Lipkin model[END_REF][START_REF] Chen | Unconventional quantum phase transition in the finite-size Lipkin-Meshkov-Glick model[END_REF].

We consider a system of N fermions, each occupying one of two levels separated by an energy ε. This situation configures a N -fold degenerate two-level system with single-particle states given by quantum numbers k = 1, 2, . . . , N for the particular degenerate state, and σ = ±1 for the upper or lower level in the state k. Taking a two-body interaction between the fermions, such that pairs of particles are scattered up or down from one level to the other without changing the values of k, this system can be described by the quantum Hamiltonian [START_REF] Lipkin | Validity of many-body approximation methods for a solvable model[END_REF]

H = 1 2 ε kσ σa † kσ a kσ + 1 2 V kk σ a † kσ a † k σ a k (-σ) a k(-σ) . (7)
The first term just counts the difference between the number of particles in the upper and lower levels, while the second term describes the change in energy occurring when a pair of particles goes from the same level -σ (with different k and k ) to the other level +σ.

Defining SU (2) pseudo-spin operators as

J + = k a † k(+1) a k(-1) = J † - (8)
J z = 1 2 kσ σa † kσ a kσ , (9)
this quantum Hamiltonian becomes

H = εJ z + 1 2 V J 2 + + J 2 -. (10)
In order to put in evidence the critical value of the interaction at QPT, we take the form [START_REF] Heiss | The large N behavior of the Lipkin model and exceptional points[END_REF]14]

H = J z + g 2N J 2 + + J 2 -, (11)
where the whole Hamiltonian ( 10) is scaled by the constant factor ε. Among the different values of J allowed by addition of angular momenta, the ground state of Hamiltonian ( 7) is realized with the eigenvalues J = N/2 and m z = -J = -N/2 for the operators J 2 and J z . The next step is the obtention of a classical Hamiltonian with the usual definitions of variables [START_REF] Lieb | The classical limit of quantum spin systems[END_REF]

j k = lim J→∞ J k J (k = +, -, z) . (12)
The relations

j x = 1 2 (j + + j -) = 1 -j z 2 cos φ j y = 1 2i (j + -j -) = 1 -j z 2 sin φ j z = cos θ (13)
make the variables p = j z and q = φ canonical conjugate variables in the sense of the classical Hamilton equations of motion. In this limit, from the quantum Hamiltonian [START_REF] Blundell | Magnetism in Condensed Matter[END_REF], we get the classical Hamiltonian

H = p + λ 1 -p 2 cos (2q) (14)
with

-1 ≤ p ≤ 1.
Figures 1(a) and 1(b) show respectively quantum results for energy spectra and the mean value J z for different values of the parameter g in eq. ( 11). Figures 1(c) and 1(d) show the classical phase space for two different values of the parameter λ in eq. ( 14).

We observe that the quantum spectra are symmetric with respect to their centers, a property that appears correspondingly in the classical phase spaces. Inflection points, characteristic of QPT in these Curie-Weiss models, are present in the spectra of the Lipkin model for g ≥ 1.0 and at energies E = ±100. These are precisely the energies of the separatrices that appear between rotations and librations, for λ ≥ 1.0. We refer to these orbits as open and closed orbits, respectively. The appearance of the closed orbits around equilibrium points (maximum at q = 0 and minimum at q = ±π/2) signals that a new phase is accessible to the system. The quantum mean value J z indeed shows a sudden rise in the collective excitation. We note, furthermore, that the inflection points which are seen in the upward half of the spectra predict a corresponding collective de-excitation for the highestenergy state.

The pairing model

The pairing model was conceived after a suggestive analogy between the spectra of nuclei and those of superconducting metallic states, where interacting pairs of particles with equal and opposite momenta are the major interest [START_REF] Bohr | Possible analogy between the excitation spectra of nuclei and those of the superconducting metallic state[END_REF]. We take here the form studied by Krieger and Goeke [START_REF] Krieger | Application of the adiabatic, time-dependent Hartree-Bogolyubov approximation to a solvable model[END_REF], in which two N -fold degenerate levels have an energy separation ε, just like in the Lipkin model. These could be two nuclear j-shells j 1 and j 2 , with states labeled as kσ, where σ = -1 for j 1 and σ = +1 for j 2 identify the shells and where k specifies the z-component of an angular momentum. In the pairing model, one is interested in the scattering of pairs of particles which are coupled with opposite values k and -k.

In order to take into account both the two levels and these pairing interactions, we write the quantum Hamiltonian Here a † and a are fermion creation and annihilation operators,

H = 1 2 ε kσ σa † kσ a kσ - 1 4 g kσ a † kσ a † kσ k σ ak σ a k σ . (15)
a † k(+1) = (-1) j 2 -k a † (-k)(+1) and a † k(-1) = (-1) j 1 -k a † (-k)(-1)
. Differently from the Lipkin model, in this pairing model the interaction between the particles may change the value of k, if a coupled pair in states (k, -k) is scattered to a new state (k , -k ). Due to this coupling, we define two sets of pseudo-spin operators, one for the upper level and another for the lower one, as

J 1+ = k>0 a † k(-1) a † k(-1) = J † 1- (16)
J 1z = 1 2 k a † k(-1) a k(-1) - N 4 ( 17
)
and

J 2+ = k>0 a † k(+1) a † k(+1) = J † 2- (18)
J 2z = 1 2 k a † k(+1) a k(+1) - N 4 . (19)
These are SU (2) sets of operators and [J 1α , J 2β ] = 0 for any α, β = -, +, z.

Applying these definitions to Hamiltonian [START_REF] Lieb | The classical limit of quantum spin systems[END_REF], we get

H = ε (J 2z -J 1z ) - g N (J 2+ + J 1+ ) (J 2-+ J 1-) . (20)
Moreover, the operators J 1 2 , J 2 2 and (J 1z + J 2z ) commute with the Hamiltonian and their mean values are constants of motion. As in the Lipkin model, the ground state here is realized for eigenvalues J 1 = J 2 = N 4 , and m z = -N 4 for both J 1z and J 2z . Since J 1z refers to the lower particle levels, with energies -ε 2 , it is a maximum for this ground state, while J 2z is a minimum. In this way, one needs only to diagonalize H in the subspace where J 1z + J 2z = 0. Also, using the same definitions ( 12) and ( 13), we get from the quantum Hamiltonian [START_REF] Ping | Decay of Loschmidt echo at a critical point in the Lipkin-Meshkov-Glick model[END_REF] the classical Hamiltonian

H = (p 2 -p 1 ) -Λ 2 -p 2 2 + p 1 2 + 2 1 -p 2 2 1 -p 1 2 cos (q 2 -q 1 ) . (21)
The following canonical transformations

p a = p 2 + p 1 2 q a = q 2 + q 1 ( 22
)
p b = p 2 -p 1 2 q b = q 2 -q 1
rewrites ( 21) as

H = 2p b -Λ 2 -(p a + p b ) 2 + (p a -p b ) 2 +2 1 -(p a + p b ) 2 1 -(p a -p b ) 2 cos (q b ) . (23)
Now, the restriction J 1z + J 2z = 0 takes the classical form p 2 + p 1 = 0 = p a , so H becomes

H = p + λ 1 -p 2 cos 2 q 2 , (24)
in appropriate energy unity, and where the index b is suppressed as it is no longer necessary. Note that, here, λ assumes negative values in order to follow the interaction assumed in [START_REF] Lieb | The classical limit of quantum spin systems[END_REF]. Again, Figures 2(a) and 2(b) show respectively quantum results for energy spectra and the mean value J 2z for different values of the parameter g in eq. ( 20). Figures 2(c) and 2(d) show the classical phase space for two different values of the parameter λ in eq. [START_REF] Bohr | Possible analogy between the excitation spectra of nuclei and those of the superconducting metallic state[END_REF].

It is seen that the same general structure found in the Lipkin model reappears with the pairing model. In this case, the spectra are not symmetric but still repeat the lower half of Lipkin spectra, with the inflections points associated with the quantum phase transition for g ≥ 1.0. Also, there are the closed orbits (librations) around an equilibrium point, for λ ≤ -1.0. The mean value J 2z measures the number of elements in the shell j 2 and shows a collective population of this shell for g above the critical value.

The N -atom Jaynes-Cummings model

The Jaynes-Cummings model [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF] with N atoms -also named in the literature as Dicke model [START_REF] Dicke | Coherence in spontaneous radiation processes[END_REF] or Tavis-Cummings model [START_REF] Tavis | Exact Solution for an n-Molecule-Radiation-Field[END_REF] -considers the interaction of N two-level atoms with a single mode radiation field of frequency ω 0 (see also [START_REF] Shore | The Jaynes-Cummings Model[END_REF] and [START_REF] Garraway | The Dicke model in quantum optics: Dicke model revisited[END_REF]). For each atom, the two levels are separated by an energy ε. Since its beginning, the model takes the N atoms (or molecules) as a single spin system described by SU (2) operators J 1+ , J 1-, J 1z with N = 2J 1 , and the radiation field described by bosonic creation and annihilation operators b, b † , and b † b. The quantum Hamiltonian for this spin-boson system, with hω 0 = ε, is written as where the first two terms count the free energies of atoms and field, and the last two terms describe the interactions between atoms and field in the rotating wave and non-rotating wave approximations.

H = ε b † b + εJ 1z + g bJ 1+ + b † J 1-+ g b † J 1+ + bJ 1-, (25)
In order to get a quantum Hamiltonian fully written in terms of sets of SU (2) operators, we apply an inverse Holstein-Primakoff transformation [START_REF] Holstein | Field dependence of the intrinsic domain magnetization of a ferromagnet[END_REF][START_REF] Carneiro | Application of an inverse Holstein-Primakoff transformation to the Jaynes-Cummings model[END_REF] on the bosonic operators b, b † , and b † b, so as to spinorize them. We define

b † b = J 2 1 + J 2z b † = J 2+ 1 √ J 2 1 -J 2z b = 1 √ J 2 1 -J 2z J 2-, (26)
and take |m = |ψ n , where |ψ n are the eigenstates of b † b with eigenvalues n = J 2 +m. Note that |ψ 0 = |-J 2 and b † b counts the number of excitations, starting at the lowest eigenvalue of J 2z .

Applying the transformation [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF] to Hamiltonian (25), we get

H = J 2 1 + (J 1z + J 2z ) + g √ N 1 √ J 2 1 -J 2z J 2-J 1+ + J 2+ 1 √ J 2 1 -J 2z J 1- (27)
+ g √ N J 2+ 1 √ J 2 1 -J 2z J 1+ + 1 √ J 2 1 -J 2z J 2-J 1-,
where, again, we use ε as the energy unity.

In the rotating wave approximation, with g = 0, the observable b † b + J 1z commutes with H and is a constant of motion; so (J 2z + J 1z ) also does. In the counter-rotating wave approximation, with g = 0, the observable b † b -J 1z commutes with H, as well as (J 2z -J 1z ). In diagonalizing Hamiltonian [START_REF] Tavis | Exact Solution for an n-Molecule-Radiation-Field[END_REF], we take advantage of these constants to reduce the size of the matrix.

Again, taking definitions ( 12) and ( 13), followed by the canonical transformations [START_REF] Heiss | On the thermodynamic limit of the Lipkin model[END_REF], the classical Hamiltonian reads

H = A + p a + λ 1 + (p a + p b ) 1 -(p a -p b ) 2 cos (q b ) (28) +λ 1 + (p a + p b ) 1 -(p a -p b ) 2 cos (q a ) .
In this expression, A is a constant with no contribution to dynamics, and λ, λ are appropriate parameters derived from g and g respectively. For the Jaynes-Cummings model, we separate two cases of interest. In the rotating wave approximation, λ = 0 and the chosen constant J 2z + J 1z = 0 leads to p 2 + p 1 = 0 = p a . In this case, we have

H λ =0 = λ (1 + p) 1 -p cos (q) . (29)
It is seen that the parameter λ is a global scale factor, so the phase space remains unaltered for any value λ > 0.

In the counter-rotating wave approximation, with λ = 0, we take J 2z -J 1z = 0 to get p 2 -p 1 = 0 = p b and

H λ=0 = p + λ (1 + p) 1 -p cos (q) . (30)
In the sequel, we present two sets of figures for these two cases (Figs. 3 and4). Both sets bring spectra and mean value of an observable as quantum results, and phase space pictures as classical results.

In the rotating-wave case, due to the commutation of J 1z + J 2z with the Hamiltonian, the interaction parameter g is a global scale factor, and the classical phase space is invariant with respect to λ. The mean value J 1z does not show any change for g > 0. One can observe that there is, in fact, an inflection point at the center of the spectra, but there is no accompanying change in stability in the classical phase space. Indeed, the mean value J 1z , calculated for the energy eigenstate at the center of the spectrum has the same constant behavior as that shown in Fig. 3(b).

In the counter-rotating-wave case, however, the whole phase-transition structure observed in Lipkin and pairing models is again present. One observes the inflection points in the lower part of the spectra, the collective excitation of the atoms for g > 1.0 and the appearance of closed orbits around the lower critical point in phase space. Note that the other upper critical point never disappears.

The bilayer model

The bilayer model, in which two separate layers are each occupied by a twodimensional fermionic gas, has been used to test Einsteins' prediction that at low temperatures all the bosons in the system should condense into the same quantum state. These bilayer models provide the necessary bosons in the form of excitons given by bound states of electron-hole pairs [START_REF] Eisenstein | Evidence for spontaneous interlayer phase coherence in a bilayer quantum Hall exciton condensate[END_REF]. In 2004, an experiment with a bilayer system under strong magnetic field described exciton condensation in electron-electron parallel layers [START_REF] Eisenstein | Bose-Einstein condensation of excitons in bilayer electron systems[END_REF]. A model Hamiltonian was proposed to mimic the dynamics of creation and annihilation of excitons formed from the N/2 pairs of fermions, as suggested by that experiment [START_REF] Moreira | Entanglement and classical instabilities: Fingerprints of electron-hole-to-exciton phase transition in a simple model[END_REF][START_REF] Figueiredo | Classical analysis of quantum phase transitions in a bilayer model[END_REF]. We rewrite it here as where b, b † stand for boson and a, a † stand for fermion creation and annihilation operators. This Hamiltonian can be written in terms of SU (2) operators only with the aid of definitions similar to [START_REF] Leyvraz | Large-n Scaling Behavior of the Lipkin-Meshkov-Glick Model[END_REF][START_REF] Rosensteel | Equations of motion for a spectrumgenerating algebra: Lipkin-Meshkov-Glick model[END_REF] for the fermionic operators

H = δb † b + g N/2 k=1 a † k1 a † k2 b + b † a k2 a k1 , (31)
J 1+ = N/2 k=1 a † k1 a † k2 = J † 1- (32)
J 1z = 1 2 N/2 k=1 a † k1 a k1 + a † k2 a k2 - N 4 , (33)
and of the inverse Holstein-Primakoff introduced for the Jaynes-Cummings model [START_REF] Jaynes | Comparison of quantum and semiclassical radiation theories with application to the beam maser[END_REF] for the bosonic ones

b † b = J 2 1 + J 2z b † = J 2+ 1 √ J 2 1 -J 2z b = 1 √ J 2 1 -J 2z J 2-. (34)
With these new operators, and taking δ = 1 as energy unity, the bilayer Hamiltonian reads

H = (J 2 1 + J 2z )+ g √ N J 1+ 1 √ J 2 1 -J 2z J 2-+ J 2+ 1 √ J 2 1 -J 2z J 1-. (35)
In this model, the total number N of fermions is constant and is given by the operator

N op = N/2 k=1 a † k1 a k1 + a † k2 a k2 + 2b † b, (36) therefore,
N op = 4J 2 1 + 4(J 1z + J 2z ) (37)
and N = 4J 2 with (J 1z + J 2z ) also commuting with Hamiltonian [START_REF] Reis | Relation between quantum phase transitions and classical instability points in the pairing model[END_REF]. We work in the subspace J 1z + J 2z = 0 in this model. Once again, definitions ( 12) and ( 13), followed by the canonical transformations [START_REF] Heiss | On the thermodynamic limit of the Lipkin model[END_REF], put the classical Hamiltonian in the form

H = (1 + p a ) + p b + λ 1 + (p a + p b ) 1 -(p a -p b ) 2 cos (q b ) . (38)
Since we choose J 1z + J 2z = 0, we get p a = 0 and, neglecting the constant term, we arrive at It is seen that the bilayer model presents results with the same structure as does the N -atom Jaynes-Cummings model in the counter-rotating case. For the bilayer model, at the critical value g = 1/ √ 2, there happens a collective formation of excitons. The difference between the quantum parameters is due to the fact that in the bilayer model one has N = 4J 2 , whereas in the Jaynes-Cummings model N = 2J 1 . Some care, however, must be exercised in comparing these models. The interaction part of the quantum bilayer Hamiltonian is, in fact, equal to the interaction part of the rotating-wave Hamiltonian for the Jaynes-Cummings model.

H = p + λ (1 + p) 1 -p cos (q) (39)

The Heisenberg model

The Heisenberg model is a spin model, usually studied in statistical mechanics to show magnetic phase transitions [START_REF] Blundell | Magnetism in Condensed Matter[END_REF]. However, it is an attractive model also in the few-body domain, where possible integrable and chaotic versions were considered [START_REF] Nakamura | Quantum irregular spectrum in a corresponding classically chaotic lattice spin system[END_REF][START_REF] Nakamura | Classical dynamics and quantum spectra for a nonintegrable three-spin system[END_REF][START_REF] Pellegrino | Quantal and Classical Manifestations of Symmetry Breaking in the Integrable Two-spin System[END_REF][START_REF] Pellegrino | On the interplay between symmetry breaking, integrability, and chaos in the semiclassical limit of the Heisenberg system[END_REF].

Here we look at the anisotropic two-spin Heisenberg model, often written as

H = F (S 1 • S 2 + σS 1z S 2z ) (40)
or, in our scheme,

H = J 1z J 2z + g 2 (J 1+ J 2-+ J 1-J 2+ ) . (41)
In this model, (J 1z + J 2z ) commutes with H and provides a constant of motion, which we take as J 1z + J 2z = 0. Relations ( 12) and ( 13), followed by the same canonical transformations [START_REF] Heiss | On the thermodynamic limit of the Lipkin model[END_REF], produce the classical Hamiltonian

H = (p a -p b ) (p a + p b ) + λ 1 -(p a -p b ) 2 1 -(p a + p b ) 2 cos (q b ) . ( 42
)
From J 1z + J 2z = 0, we have p a = 0, so we end with

H = p 2 + λ 1 -p 2 cos (q) , (43)
where λ and g have opposite signs. Figure 6, showing quantum and classical results for the Heisenberg model, appears below. Observe that the quantum mean-value curve is symmetric, due to the quadratic term J 1z J 2z in Hamiltonian (41). The Heisenberg model shows a richer structure, which, in a sense, combine characteristics seen in the rotating-wave case of the Jaynes-Cummings model with those seen in the other models. The spectra present the usual inflection point associated with quantum phase transition -observed in the mean value J 2z (Fig. 6(b)) -and another one at energy zero for g > -1.0. Accordingly, there are two types of critical points in the classical phase space, but only one of them reflects QPT, with a corresponding change in stability. One type, at q = ±π, is always present and its surrounding librations are confined by a separatrix orbit -as seen in Fig. 6(c) -while the other one, at q = 0, appears as λ goes above 1.0, with closed orbits around it and the -now rectangular -separatrix between the critical points. Therefore, Fig. 6(c) is associated with the mean value J 2z for g > -1.0, and Fig. 6(d) with that part with g < -1.0. Once they exist, these critical points are static with respect to variations in the interaction parameter. This is a different behavior in comparison with the other models, where the critical point moves up or down as a function of λ. As a final remark, for values g > -1.0, the inflection at energy zero has precisely the energy of the separatrix seen in Fig. 6(c). However, like the inflection point found in the rotating-wave case of the Jaynes-Cummings model, it is not accompanied by a change in stability in phase space.

Discussion

A number of universal characteristics have been pointed out in the literature, concerning quantum phase transitions (QPT) in these Curie-Weiss models [START_REF] Heiss | On the thermodynamic limit of the Lipkin model[END_REF][START_REF] Reis | Relation between quantum phase transitions and classical instability points in the pairing model[END_REF][START_REF] Nemes | Quantum entanglement and fixed point Hopf bifurcation[END_REF][START_REF] Emary | Quantum Chaos Triggered by Precursors of a Quantum Phase Transition: The Dicke Model[END_REF][START_REF] Emary | Chaos and the quantum phase transition in the Dicke model[END_REF][START_REF] Moreira | Entanglement and classical instabilities: Fingerprints of electron-hole-to-exciton phase transition in a simple model[END_REF][START_REF] Figueiredo | Classical analysis of quantum phase transitions in a bilayer model[END_REF][START_REF] Pellegrino | Quantal and Classical Manifestations of Symmetry Breaking in the Integrable Two-spin System[END_REF]. In the following, we comment on some of them.

QPT is signaled by the presence of an inflection point in the spectrum. As long as the interaction parameter g is varied, this inflection point moves upwards in the spectrum and a critical parameter for QPT is usually taken as the lowest energy level becomes the inflection point, which then disappears. This situation is seen in all models here, except for the Jaynes-Cummings model in the rotating wave approximation. Also, as the lower excited energy levels become the inflection point, the corresponding quantum states manifest QPT. This is a clear illustration of a more general definition of QPT, known as excited state quantum phase transition (ESQPT) [START_REF] Caprio | Excited state quantum phase transitions in many-body systems[END_REF][START_REF] Cejnar | Impact of quantum phase transitions on excitedlevel dynamics[END_REF].

At QPT, there happens a conspicuous increase in the mean value of an observable associated with the number of one type of the elements in the systems. Again, this is seen in all but one of our cases.

Finally, the energy of the inflection point in the quantum spectrum is the energy of a separatrix orbit in the classical phase space, and the appearance of this separatrix, as the classical interaction parameter λ crosses the critical value, marks the transition. In addition to these characteristics, the unified approach presented here allows a number of new observations and comparisons.

Starting with the obvious resemblance in the curves for the mean values (part (b) in the figures), one realizes that the lower half of the spectra have the same structure. This resemblance shows up also in the structure of the classical phase space. In fact, since the open orbits are clearly associated with the free term H 0 of the quantum Hamiltonian, and the closed ones with the interaction term H int , the competition between these terms responds for the mechanism under which QPT takes place. In this way, the appearance of closed orbits around an equilibrium point -a sudden change in the global (in)stability of phase space [START_REF] Hines | Quantum entanglement and fixed-point bifurcations[END_REF][START_REF] Nemes | Quantum entanglement and fixed point Hopf bifurcation[END_REF] -signals that a new quantum regime is accessible to the system as a whole, the collective formation of excitons in the bilayer model or the collective excitation of fermions in the Lipkin model, as examples.

Moreover, a unified approach may bring up interesting questions, concerning physical effects in different systems with similar behavior in the general scheme. Counter-rotating Jaynes-Cummings model (eq. ( 30)) and bilayer model (eq. ( 39)) have both the same classical expression. Do the physical effects seen in one system have analogs in the other? This same question arises for the pair formed by Lipkin model (14) and pairing model [START_REF] Bohr | Possible analogy between the excitation spectra of nuclei and those of the superconducting metallic state[END_REF]. In the Heisenberg model, QPT assumes a somewhat different classical manifestation. The (dis)appearance of closed orbits around a critical point does not follow the movement of the equilibrium point with respect to the border of the phase space. This is an otherwise static critical point, with respect to the interaction parameter. Static critical points are present also in the rotating-wave Jaynes-Cummings model but, in this case, there is no change in classical phase space neither in the quantum mean value. Are there physical effects differentiating the quantum phase transitions in this Heisenberg model and the other ones?

Conclusion

In this work, we propose a unified approach to the analysis of quantum phase transitions in some different Curie-Weiss models. Two points are distinctive in the proposal: the balancing of free and interaction Hamiltonians by a proper factor, and the transformation of fermionic and bosonic creation and annihilation operators into SU (2) operators. The first point avoids an eventual QPT from occurring for an interaction parameter value too close to zero, which would make difficult its observation. The second point puts all these models under the same general scheme, which permits treating different physical situations, making comparisons, and raising analogous questions for different areas.

We do not present here a detailed analysis of the physical results for these models, mainly of the possible translations of physical questions usually present in one model but not in a similar other, as pointed out by their general form in this scheme. Work along this line is in progress and will appear in due future.

Figure 1 :

 1 Figure 1: Lipkin model: (a) quantum spectra for interaction parameter values g = 0, 0.5, . . . , 3; (b) quantum mean value J z ; (c) classical phase space for λ = 0.5; (d) classical phase space for λ = 1.5. The phase transition occurs at g = 1.0 and λ = 1.0.

Figure 2 :

 2 Figure 2: Pairing model: (a) quantum spectra for interaction parameter values g = 0, 0.5, . . . , 3; (b) quantum mean value J 2z ; (c) classical phase space for λ = -0.5; (d) classical phase space for λ = -1.5. The phase transition occurs at g = 1.0 and λ = -1.0.

Figure 3 :

 3 Figure 3: Rotating-wave Jaynes-Cummings model: (a) quantum spectra for interaction parameter values g = 0, 0.5, . . . , 3; (b) quantum mean value J 1z ; (c) classical phase space for any λ > 0.

Figure 4 :

 4 Figure 4: Counter-rotating-wave Jaynes-Cummings model: (a) quantum spectra for interaction parameter values g = 0, 0.5, . . . , 3; (b) quantum mean value J 1z ; (c) classical phase space for λ = 0.7; (d) classical phase space for λ = 0.8. The phase transition occurs at g = 1.0 and λ = 1/ √ 2.

Figure 5 :

 5 Figure 5: Bilayer model: (a) quantum spectra for interaction parameter values g = 0, 0.5, . . . , 3; (b) quantum mean value J 2z ; (c) classical phase space for λ = 0.7; (d) classical phase space for λ = 0.8. The phase transition occurs at g = 1/ √ 2 and λ = 1/ √ 2.

Figure 6 :

 6 Figure 6: Heisenberg model: (a) quantum spectra for interaction parameter values g = -3, -2.5, . . . , -0.5; (b) quantum mean value J 2z ; (c) classical phase space for λ = 0.5; (d) classical phase space for λ = 1.5. The phase transition occurs at g = -1.0 and λ = 1.0.

We follow the established literature and refer to this model as the Lipkin model.

Acknowledgements

The authors wish to express their warm gratitude to José Geraldo Peixoto de Faria and Helen Barreto Lara for many illuminating discussions. C.M.P.C. thanks to CAPES-Brazil.


Y S Choi

S T Noau

FORCED PERIODIC VIBRATION OF UNSYMMETRIC PIECEWISE-LINEAR. SYSTEMS

The forced steady response of a single degree of freedom system involving a large non-linearity, represented by unsymmetric piccewise•linear stiffness, is determined by a harmonic balance Newton-Raphson method with the application of the fast Fourier transformation (FFT) algorithm. All possible subharmonic, harmonic, and superharmonic responses arc sought. The responses of two systems involving subharmonic and superharmonic dominant motions are calculated by the newly developed method and compared with previously published findings. The results obtained by using this method reveal the details of the response of the systems more efficiently than previous methods.

INTRODUCTION

Vibration problems involving small non-linearities lend themselves to closed form solutions obtained by using perturbation techniques. \Vhen the nonlinearities cannot be assumed small, general solution methods are usually not available. The dynamic charac• teristics of general non-linear systems are too complex to predict as compared to those of linear systems. For example, the response to periodic excitation may result in multiplevalued response, chaotic response or jump phenomena.

A common solution procedure for strong non-linear vibration problems is the direct numerical integration. This procedure can yield both transient and steady state response for given initial conditions and specified frequencies. The disadvantages of this approach include the following: (i) it would be too expensive and tedious to completely characterize the dynamic behavior of the system in relatively large ranges of their various parameters; (ii) in a multiple-valued response region, a11 possible solutions could be found only by changing the initial conditions in a trial and error process; (iii) to guarantee the convergence of the solution, the time-step integration method would have to be chosen carefully; this is one of the drawbacks to direct numerical integration.

The determination of closed form solutions for piecewise-linear systems was attempted by Maezawa [ 1] using Fourier series expansions to represent stable, unstable and multivalued responses. Unfortunately, a solution method based on assumed response patterns often suffers restrictions. Unprcdicted types of vibration, not included in the assumed response shapes may be excluded in these methods. Furthermore, it is almost impossible to find the exact duration of contact, 0 0 , in some ranges of 8 0 (which was treated as a parameter in Maezawa•s solution procedure) because of the s~vere non-linear nature of the relations involving 0 0 • Yamauchi [START_REF] Yamauchi | [END_REF] and Saito [3] used fast Fourier transformation (FFf) to determine successfully the non-linear response of a modified Jeficott model for a rotating disk. Their solution did not include possible subharmonic response which might occur. Saito finally solved a set of iterative linear simultaneous equations in complex variable form which were produced from the Taylor series expansion. The complex form introduced unduly complicated the solution procedure.

In this paper, the dynamic characteristics of an unsymmetric piecewise-linear springmass system under harmonic excitation is examined. All possible steady periodic, subharmonic, superharmonic, and harmonic responses are included. To simplify the calculation of the simultaneous complex equations only the real variables are used instead of the complex variables used in reference [3].

The calculated subharmonic response by the newly developed method is compared with the results of direct numerical integration. Also, the superharmonic response of an unsymmetric stiffness system with an offset is obtained and compared with the results given in reference [ 1 ]. The results for both cases closely agree with the previously published results, and in addition, provide information about the unstable periodic solutions which cannot be obtained by a numerical integration procedure.

EQUATION OF MOTION

The system considered consists of a mass, a viscous damper, and an unsymmetric piecewise-linear spring system, and is excited by a harmonic force. The equation of motion can be written as

(I)
where m, c, t, and x denote the mass, viscous damping coefficient, time and displacement, respectively. The parameters A and w denote the amplitude and frequency of the forcing function. The unsymmetric piecewise-linear force, g(x ), is depicted in Figure 1 and is expressed as follows:

x<h} X~ lz ' (2)
where K 1 and K 2 are the spring constants, h JS an offset value, and 5 is defined as

5 = (K 1 -K 2 )h.
There is no restriction on the magnitude of the K 1 and K 2 values, so that g(x) exists for all real x, although the derivative of g(x) with respect to x is of a finite discontinuity at x =h. Because of this arbitrary piecewise-linear nature of g(x) the equation of motion (I) describes a system involving a large non-linearity.

For the free undamped vibration of the system and no offset (h = 0), the period of this unsymmetric system will be composed of two half sine waves, so that the assumed period, T, can be expressed as T= rr(Jm/ K,+Jm/ K 2 ). The bilinear radian frequency can therefore be written as w" =2rr/T=JK/m, where K is an equivalent stiffness, defined as (4) In order to determine the subharmonic vibration of order 1/ v, v being a positive integer, a dimensionless time, 0, is adopted as the new independent variable,

vO = wt, ( 5
)
where v is taken as unity when the harmonic or superharmonic vibrations arc sought.

The period of the subharmonic vibration will be v-times the period of the exciting force. \Vith the above equivalent stiffness and the consideration of subharmonic vibrations, equation ( 1) can be non-dimensionalizcd by introducing the following parameters, as in reference [ 4]: a stiffness ratio, a:::: K 2 / K., a damping ratio,~= c/2mwn, and a frequency ratio fJ = w/ Wn. Equation (1) can then be written as

2TJ

),

[ (I +.Ja]2 t:* r- y+ u , 2TJv a y<{* y=xK/A, g*=hK/A, (7)
and 8* is defined so as to make g* continuous at s*. The prime denotes differentiation with respect to the scaled time, 0.

STEADY STATE SOLUTION

The steady periodic solution of equation ( 6), including subharmonic, superharmonic and harmonic vibrations can be written as the Fourier series

N y(O)= a 0 +2 ~ (an cos nO-b" sin 110), ( 8
) n-1
where N is the number of harmonics to be taken into account under the assumption of small bandwidth in the final solution. The multiple of 2 and the minus sign in bn are adopted to facilitate accommodating complex Fourier coefficients when FFf is used. In the same way, the non-linear restoring force, g( 0), can also be expressed as a Fourier series of the following form:

N g*(O)= e 0 +2 L (en cos nfJ-dn sin nO). ( 9
)
n=l Since this non-linear restoring force was derived from the piecewise-linear form of the stiffness in equation {6), each of e 0 , e 11 , and dn will be functions of a 0 , a,, and bn.

Substituting equation ( 8), and equation ( 9) into equation ( 6) and applying a harmonic balance procedure yields 2N + 1 non-linear simultaneous equations, involving 4N + 2 unknowns. The harmonic balance procedure gives, for the constant term,

v 2 c 0 = 0,
for the cosine terms, and for the sine terms, n = I, ... , N,

jf ll ~ V} ifn=v n=I, ... ,N. (lOa) (lOb) ( IOc)
Since Co, cM 'and d" arc functions of G 0 , an, and bn, then once a 0 , an, and b" arc determined, c 0 , c", and d" can be readily calculated based on the relations given by equation ( 6). This means that if a 0 , a,., and b" are given, y( 8) can be found, and with the relation of equation (6) g*(O) can then be generated. The application of the FFT procedure could generate c 0 , c"' and d" through its efficient computational algorithm.

Therefore, the relation given by equation ( 6) with the application of the FFf algorithm could regenerate the 2N +I equations which gives the relations implicitly between the Fouricr coefficients of y( 0), a 0 , an, and bn and those of g*( 0 ), c 0 , c,, and d,. Consequently, equation ( 10) would have 2N + l unknowns, a 0 , a", and b".

The resulting non-linear simultaneous equations can be handled by using a Newton-Raphson iteration method. If initial values are given properly, the solution converges very rapidly. A guess of initial values for the unknowns G 0 , an, and bn can be made with a linear approximation. From equation (6) the linear steady state solutions can be found with equivalent stiffness as shown in equation ( 4) instead of the piecewise-linear stiffness. In the case of multiple steady state solutions for a given frequency ratio, all solutions found from the equivalent stHTness calculation arc pursued as possible starting points for the harmonic balance method. Small increments .1a 0 , Ja", and Jbn arc assumed to be close to the solution in question, so that one can write Gn =a~+ aa,, (11) where the superscript zero indicates initial guessed values for the first iteration and old values for the next iteration.

Similarly, the nonlinear restoring force terms eo, en, and d" can also be expressed as follows:

(12) Since c 0 , c", and d" depend on ao, an. and b"' the increments can be expanded as

(13)
where the partial derivatives are to be calculated at the current state value. In order to account for the large non-linearity of the problem, the increments ~a0 , .Jb"' and ~b" must be chosen comparatively smaller than a 0 , ant and b" whenever numerical differentiation is performed.

The numerical differentiation procedure can be done by using forward differentiation as follows: (I) a 0 is changed by adding the preset small Lla 0 , but other Fourier coefficients, a"' b" arc not changed; (2) the newly assumed discrete displacement, y( 8) is generated by using a discrete inverse Fourier transform; (3) the discrete non-linear force, g*(O), is calculated by using equation ( 7); (4) the Fourier coefficients of the discrete non-linear force, c 0 , c"' and cl, are calculated by using the fast Fourier transformation (FFf); (5) the differential values of oc 0 / iia 0 , Ac;/ c1a 0 iid;/ ria 0 are caJculated as ac,/iJao = (c,-c?>l .1.ao, where i =I, 2, ... , N; (6) the calculations in 1 to 5 are repeated for each a, and b,.

Substituting equation (11), equation ( 12), and equation (13) into equation (10) results in the following incremental equations:

2 o 2(v o 2 o 11:! = -n b + IJ -a -,, d + 0 • 5 - n t1 n 2'

T] T]

( 15)

Here 0•51' 2 /7] 2 is only included when n equals to v. Equation (15) can be put in the matrix form where rlco ila 0

1 ilc 1 "-ria 0 l r1d 1 "- [K) = ila 0 2 ilc,.... 1 1 - ila 0 2 ild,'\1 " - aa 0 [ K ]{~y} = {R}, ilc 0 ilc 0 ilc 0 i'la• ilbl ila,, 2 ,1c 1 2{v , ac, 2 ac.
-I+" ---+ v • -

1 1 - ila 1 '1 ilbl a aN 2lv 2 ild 1 2 fld, 1 fld, --+ v - l+" - " - '1 ila 1 ilbl a aN . . 2 ilcN 2 .lcN z 2 ilc, .. , " -
,.

-

-N + v - a a, ab,
,la,...

{R}=

Here [ K] corresponds to the Jacobian matrix whose elements are calculated at every step. This iteration will continue until all the components of the correction vector { R} become sufficiently close to zero. The calculation procedure presented here is basically an iterative process. Initial guessed values for a new case can be obtained by choosing those previously computed values for the conditions close to the new case. To find multiple solutions, the initial values are calculated differently for the assumed Jinear equivalent stiffness case since the nondimensionalized equation ( 6) will be different depending on the subharmonic ratio v. Then, for a certain case the response will trace a different path. The method will yield both stable and unstable possible periodic solutions.

NUMERICAL RESULTS AND DISCUSSION

A system with stiffness of a> 1 was considered as the simple unsymmetric bilinear system of Figure I for zero gap length. Figure 2 shows an overview of the response of this system. The numerical results obtained by the FFf method are compared with those obtained by using direct numerical integration. The average steady state ampJitude, y, defined as half of the peak-to-peak of the response y, is plotted as a function of the frequency ratio, 17, with the damping ratio, C. and stiffness ratio, a, fixed.

The direct numerical integration results were obtained by the fourth-order Runge-Kutta method with initial conditions, y 0 and y~, set to zero. The direct numerical integration continued until the motion reached steady state periodic motion, so that y(21Tj) = y(2 ;r(j + ")),

y'(2rrj) = y'(21T(j+ v)), ( 17
)
where v is the order of the subharmonic vibration and a prime denotes differentiation with respect to the dimensionless time, 0. In this numerical marching scheme, the time increment should be smal1 in order to find exact contact time at every contact point. Any large deviation of the contact time will result in a false steady state solution. At some range of TJ the time increment is given by one five-hundredth of the exciting period.

Figure 2 shows good agreement between the results obtained by using the FFf method and those obtained by direct numerical integration. The only difference is that the harmonic response is shown to exist as given by the FFT approach throughout the whole range of frequency ratios where the harmonic response is stable or unstable. This is because the direct numerical integration can give only stable solutions. Also, the results of direct numerical integration depend on given initial conditions, making it difficult to find all possible solutions. The FFT method yields, on the other hand, all possible equilibrium solutions without consideration of the initial conditions. ------------------------------------, To find an adequate number of harmonic components to be retained for the analysis, numerical experiments were carried out and their results arc shown in Figure 3. The steady state second order subharmonic response at a frequency ratio of 2 was found by direct numerical integration, and by the FFT method with 2, 4, and 8 harmonic components. In order to avoid the aliasing effects, the number of discrete data were taken as four times that of the harmonic components. As shown in Figure 3, when 4 or 8 harmonic components were considered, good agreement between the results with those of direct numerical integration has occurred. The result shows that four harmonic components are sufficient for this case.

1•4 1:>.. 1•2 c ~ 1•0 Q) 0 0 0.. -~ "C Cl) 08 g 0•6 ~ Q.l > et 0-4 0•2 Frequency roti0,'7]
The variation of subharmonic response with respect to variation in the stiffness ratio is shown in Figure 4. The figure shows that only a harmonic response exists in the case where the stiffness ratio is unity, i.e., linear vibration, and no fourth order subharmonic response exists for a stifiness ratio of five. This behavior was also shown by Thompson and Elvcy [ 4]. ------------------- Frequency rotio, 71

Figure 4. Effect of stiffness ratio as determined by FFf method.'; 0 • 1; •• -,a= 10; --, o = 5; +,a= I.

1 75 ---------------------------------- The dependence of the response on the damping ratio is shown in Figure 5, where the stiffness ratio is fixed as 10 and the damping ratio takes the values of 0•1, 0• 3 and 0• 5. The figure shows clearly that higher damping suppresses the higher order subharmonic responses.

1•5 •~ l 1•25 ~ Q. 1•0 ~ '-a 0•75 . • ~ 0 05 ... • > <1. 025 0•0 1•0 }\ I \ I I 1'._1 I t•0•1 r v \ I I ,, \ \ \ \ t•0•3 /'\ J \ I ' ' ' ....
Cases involving supcrharmonic responses were also examined and compared with the results of Maeza\va [ 1]. The system considered has an unsymmctric piecewise-linear stiffness, an offset, and no damping. As shown in Figure 6, the results obtained by using the FFf technique are very close to those obtained by Maczawa. The results show the superharmonic resonance to occur in the neighborhood of0•35 and 0•55 frequency ratios. The shift of the main resonance near the frequency ratio of unity is due to the adoption of K 1 as a reference stiffness for the calculation of the natural frequency. As compared to the complicated derivations and lengthy calculations by Maezawa's approach, the FFf method employed in the present work proves to be simpler and more versatile. In Figure 7, a comparison between the results obtained by the two methods is made for different non-dimensionalized values of A/ K 1 d. The trend shows close agreement with that of Maezawa.

CONCLUSION

A numerical algorithm for determining the steady forced vibration of a highly non-linear system in which an FFf technique is used has been described. The algorithm is comparatively easy to formulate and results in saving of computer time as compared to that for a direct numerical integration method.

The newly developed method was applied to unsymmetric piecewise-linear systems. The approach yielded results which are in good agreement with previously published findings, and can be used eiTectively to determine all possible steady state solutions. For subharmonic resonanccs, the results show that reducing the stiffness ratio or increasing the damping ratio would lead to suppressing the higher order subharmonic responses. The results also show that the superharmonic responses of an unsymmetric piecewiselinear system occur only in some specific frequency range.

The approach developed here is being extended to the analysis of non-conservative multi-degree-of-freedom systems.

Figure

  Figure I. Unsymmelric piecewise-linear restoring force.

  ./a)]2 1

  {Jy} = (Jau. Ja 1 , Jb 1 , . . . , Ja...,, Jb.d r

Figure 2 .

 2 Figure 2. Comparison bel\\•een numerical and FFT resuhs in frequency domain. 0, Numerical solution; --, FFf solution.

Figure 3 .

 3 Figure 3. Ellcct of retained number of Fourier coefficients.-, Numerical solution; ---. N ""2; 0, N :::: 4; +, N=8.

Figure 5 .

 5 Figure5. Effect of damping ratio as determined by FFT method. a= 10.

, 0 0~Figure 7 .

 07 Figure 7. Effect of non.dimcnsionalized displacement. A/ K 1 d, in supcrharmonic resonance.+, AI K.f d == 2; 0, A/ K. 1 d::: 1.

ACKNOWLEDGI\1ENT

This work is being carried out as part of a research project supported by NASA-Marshall space flight center under contract No. NASS-36182. This support is gratefully acknowledged. The authors arc also grateful to Thomas Fox of NASA for his input as to similar methods used in non-linear control analysis and for his interest.