Stability and Boundedness of the Solutions of Multi-Parameter Dynamical Systems with Circulatory Forces

: We consider a linear dynamical system under the action of potential and circulatory forces. The matrix of potential forces is positive deﬁnite, and the main question is when the circulatory forces induce instability to the system. Different approaches to studying the problem are discussed and illustrated by examples. The case of multiple eigenvalues also is considered, and sufﬁcient conditions of instability are obtained. Some issues of the dynamics of a nonlinear system with an unstable linear approximation are discussed. The behavior of trajectories in the case of unstable equilibrium is investigated, and an example of the chaotic behavior versus the case of bounded solutions is presented and discussed.


Introduction
In the theory of dynamical systems, interest in the properties of symmetry is steadily growing. In recent decades, a systematic study of systems possessing the so-called reversible symmetry has become widespread [1][2][3][4][5]. Consider an autonomous dynamical system with continuous time, governed by the ordinary differential equation: dx It is said [1] that the map S : R n → R n is a reversing symmetry of (1) if: In this case, System (1) is invariant with respect to transformation (t, x) → (−t, Sx). The stability problem for non-conservative dynamical systems is one of the fundamental problems in modern science and technology. It is related to many real-life applications and is mainly matched with the response of mechanical systems to parametric excitation and flutter phenomenon exhibited by aerospace systems undergoing interaction between a structure and fluid flows, as well as playing a crucial role in the control of walking robots and many others [6][7][8][9][10][11].
behavior. Although the motion is unstable in the Lyapunov sense, it can however be bounded, which can prove to be acceptable enough in engineering applications.

Formulation of the Problem
The non-conservative undamped linear systems are mostly presented in the following form: where M and K are positive definite symmetric real square matrices, N is real skew-symmetric, and x ∈ R n represents a state vector. In order to underline the absence of damping, such a system is called an undamped non-conservative or circulatory system. Because 3 is a linear autonomous ODE system, its characteristic polynomial f (λ) has order 2n and contains only even degrees of λ. Consequently, the system (3) is stable if all roots of the polynomial: f (Λ) = a n Λ n + a n−1 Λ n−1 + · · · + a 1 Λ + a 0 (4) are real and negative, Λ = λ 2 . The polynomial (4) is referred to as the reduced characteristic polynomial related to (3). Three main cases are possible: (1) The system is marginally or weakly stable (Lyapunov stability) if all solutions x(t) of (3) are bounded for t ∈ [0, +∞).This means that every eigenvalue λ j is on the imaginary axis and semi-simple (there are s linearly independent eigenvectors associated with the eigenvalue of multiplicity s); (2) The system (3) is unstable by flutter (or oscillatory unstable) if at least one of the eigenvalues λ j is complex with the nonzero real part (there exists an oscillating motion with exponentially growing amplitude); (3) The system (3) is unstable by divergence (or statically unstable) if all of the Λ are real and at least one of them is real positive (there is an aperiodic, exponentially growing motion).
The negativeness of the roots Λ j (j = 1, · · · , n) may be validated with the help of Routh-Hurwitz or Lienard-Chipard criteria, and the absence of complex roots may be detected by at least three approaches. Below, we give their brief description.
Theorem 1. [29,37]. Given a polynomial f (Λ) with real coefficients, a necessary and sufficient condition for the polynomial to have only real roots is that the elements of the discriminant sequence are all non-negative D j ≥ 0, j = 1, · · · , n.
The second approach (B-approach [31,32]) is connected with consideration of the quadratic form x τ Px, where: Note that the elements p j are connected with the coefficients of f (Λ) by the following formulas: p j = −(p j−1 a n−1 + · · · + p 1 a n−j+1 + ja n−j ), j = 1, · · · , n.
(c) The system (3) is unstable by flutter iff x τ Px, can take negative values.
For example, for n = 3, these formulas take the form (with a 3 = 1): For application purposes, the following corollary is helpful: The number of real roots of the polynomial (4) (with real coefficients) different from each other is equal to µ S − ν S , where µ S is the "sign permanence" index and ν S is the "sign variability" index (i.e., quantities of the values with the same sign and changes of the sign, respectively) in the sequence (in fact, the first element of the sequence is one, and it is followed by sequentially recorded principal minors ∆ (j) S of matrix S till order 2r − 2): Here, r is the rank of the Hankel form: for polynomial f (λ).

Main Results: Circulatory System with Three Degrees of Freedom
Let us consider the system (3) with the following matrices: Firstly, we try the B-approach. The coefficients a j are: 12 + n 2 13 + n 2 23 , a 0 = k 1 k 2 k 3 + k 1 n 2 23 + k 2 n 2 13 + k 3 n 2 12 .
The expressions for elements p j , according to their definition, are: where * is the Frobenius norm of a matrix * , i.e., K = ∑ n j,s=1 k 2 js . These formulas give some idea of how the matrices K and N influence the elements p j , and allow describing stability/instability conditions in terms related to the structure of these matrices (for a 2-DOF system, such conditions were given in [19,31]). However, when n ≥ 3, for computational purposes, it is more convenient to express p 4 by the coefficients of f (Λ) and after that by p 1 , p 2 , p 3 according to Formula (6). Because all coefficients of f (Λ) are positive, Case (b) of Theorem 2 is unachievable, and the semi-positiveness of matrix P defines which, Case (a) or (c), takes place. With this observation, we have two expressions to analyze: ∆ 2 = 3p 2 − p 2 1 and: Obviously, if p 2 ≤ 0, then System (3) is unstable. Noting, in addition, that ∆ 3 may be presented in the form: and if ∆ 2 ≤ 0, then ∆ 3 ≤ 0 as well. This means that the positiveness of ∆ 3 implies the positiveness of ∆ 2 , and Case (a) takes place if p 2 > 0, ∆ 3 ≥ 0. Considering the expression (13) as a polynomial of second order with respect to p 3 ; its positiveness leads to the following inequalities: and requirement p 2 > p 2 1 /3 or: stands for the necessary condition of stability.
To compare the effectiveness of the B-approach, let us apply the technique of the G-approach. Taking into account the formulas (9), we get: S is a quadratic polynomial on a 0 with a negative coefficient. Then, independently of other coefficients, the necessary condition of the positiveness of ∆ S is the non-negativeness of its discriminant (with respect to a 0 ): Because of this fact and taking into account that a 1 > a 2 2 /3, both ∆ S and ∆ S are negative. This means that there is no possibility of sign sequence is not equal to zero, and there are only two possible sign combinations: {+ + ++} and {+ + −−}. In the first case, µ = 3, ν = 0, so f has three distinct real roots. The second case follows as µ = 2, ν = 1, and f has one real root and a pair of complex conjugate roots.
Eventually, we come to the single condition of stability ∆ S > 0. Substituting the coefficients (11), we have: The current expression is large enough for direct analysis, though it has a symmetric form (is invariant toward cyclic permutation of the indexes 1, 2, 3). Suppose that: Introducing the parameter transformation: and regarding: we come to the following expression for discriminant (with respect to λ): Being considered as the quadratic polynomial with respect to p 2 , discriminant D may be positive only if expression: is positive. Then, the condition D > 0 is equivalent to: where parameters b 1 , b 2 , q 2 are connected to the coefficients of System (3) by Formula (15). The bounding surfaces from (17) are presented in Figure 1. Finally, when d = 0, the polynomial f has three real roots, i.e., one of multiplicity two, if ∆ To demonstrate the approach suggested, let us consider an example of a mechanical system. There are not many systems in the form of (3) with uncertain parameters for n > 2 in the literature (with the analytical study of stability conditions) [30,38], so in order to underline the key-points, we take an example from [30]. Example 1. Let us consider a three link pendulum subjected to the follower force. The last has two components: tangential P F and constant non-tangential P C (Figure 2, non-dissipative case). In the small vicinity of equilibrium θ j = 0, j = 1, 2, 3, the linearized motion equations lead to System (3) with matrices: Figure 2. Dynamic model of the system subjected to a follower load (adopted from [30]).
The dimensionless characteristic equation contains two parameters: α and β = Pl/k. By the help of the MD-approach, the stability domain was obtained, and six expressions a 2 (α, β), were involved in the process. Using the G-approach, only three expressions a 2 , a 2 a 1 − a 3 a 0 , and D = D 3 are needed. In order to bring some contrast to the results of [30], we assume that P C = 0 (hence, α = 1), but at the same time, the stiffnesses of the joints k 1 , k 2 , k 3 differ from one another. Introducing the dimensionless parameters and time: we have the following reduced characteristic polynomial: The domain of stability is formed as the intersection of three sets: a 2 > 0, and D f > 0 (the corresponding expression is given in Appendix A). The discriminant is the polynomial of sixth order and does not allow for any analytical reduction. Nevertheless, the corresponding 3D surfaces are smooth enough, and plotting does not present difficulties. The results are presented in Figure 3.
With the standard assumption of the equal stiffnesses at the joints k 1 = k 2 = k 3 , the value of β crit ≈ 1.305 (Figure 3). With different values of k j , this value can be markedly increased (up to 1.5 times). The optimal configuration is when k 3 is bigger and k 2 is smaller than k 1 (see Figure 4).

Special Case: Multiple Eigenvalues of Matrix K
This section was inspired by the recent publication of F. Udwadia [34], where he generalized the well-known theorem of Merkin [24]. The last one states that when matrix K has eigenvalue Λ of multiplicity n, an arbitrary small circulatory force brings flutter instability to the system. In order to generalize this theorem, as it was pointed out by Udwadia, the commutativity of matrices K and N is required. Note that the condition of instability in case KN = NK was first obtained by Bulatovic [29]. Definition 3. [34]. Consider an n-by-n block diagonal matrix: in which the s-th (square) diagonal block A s has dimensions i s by i s with i 1 ≥ i 2 ≥ · · · ≥ i k , and another n-by-n block diagonal matrix: for which the p-th (square) diagonal block B p has dimensions j p by j p with j 1 ≥ j 2 ≥ · · · ≥ j r . We shall say that matrices A and B have the same block diagonal structure if: Theorem 4. [34]. Let K = λ 1 I 1 + λ 2 I 2 + · · · + λ k I k , and T be the orthogonal transformation such that the diagonal matrix Ω = T τ KT has repeated eigenvalues. The matrix N = diag(N 1 , · · · , N k ) that has the same block diagonal structure as Ω and where N j , j = 1, · · · , k are arbitrary skew-symmetric matrices will always commute with Ω, and the matrix N T = TNT τ will always commute with K. Arbitrarily small elements in the matrices N j will lead to a flutter instability of the systemÿ + (K + N T )y = 0.
Though this result is interesting from a mathematical viewpoint, at first blush, it seems not to be very helpful in direct applications. The reason is that for the system (3) (as M = I), the requirement of the same block diagonal structure implies the decomposition of the whole system into unstable and stable (probably) subsystems. This observation yields the question: Why not consider each subsystem separately?
However, on second thought, the sense of this result opens the possibility to detect the instability for non-linear dynamical systems, especially systems of high dimension and/or multi-parameter ones. Here, the instability of a small (say, k = 2, 3) subsystem as a rule leads to the instability of the whole system, which may lead to serious reduction of the computational effort. Bearing this idea in mind, we shall try to introduce some modifications of the above-mentioned result based on the following way.
Let B = (b 1 , b 2 , · · · , b n ) τ be a rectangular matrix and b j (j = 1, · · · , n) the rows. We "bisect" matrix B into two sub-matrices with q and n − q rows, respectively, and shall use the notation B q , B \ B q for them below.
Suppose that system MKN allows the representation: where matrices M j , K j (j = 1, 2) are square symmetric and positive definite and N, N j for skew-symmetric matrices. The dimensions of matrices M 1 , K 1 , N 1 are equal to l ≤ n, and symbol " " denotes the direct sum of matrices. Let further J M 1 be the Jordan normal form of the matrix M 1 , i.e., is symmetric and positive definite. Consequently, there exists an orthogonal matrix T K 1 such that The following sufficient condition of instability for system MKN is given by the following theorem. has the eigenvalue (semi-simple) of multiplicity q). Denote by T M 1 and T K 1 orthogonal matrices, the columns of which are formed by the eigenvectors of the corresponding matrix, and the first q columns are associated with multiple eigenvalues. Let the matrix equality: hold. Then, if: the system MKN is unstable (flutter instability).
Proof. We remind that as block matrices M 1 , K 1 are symmetric and positive definite, there eigenvalues determine the transformation matrices T M 1 , T K 1 , and: where J M 1 J K 1 are Jordan normal forms of those matrices.
These forms are diagonal, and without loss of generality, we may count left upper cells of dimension q as containing proportional elements corresponding to eigenvalue Λ 1 . As matrices T M 1 , T K 1 are orthogonal, then their inverse matrices are equal to transpose ones; therefore, systems: M 1ẏ + (K 1 + N 1 )y = 0 and: are equivalent, i.e., have the identical fundamental matrices. According to the assumption, we have J K 1 = Λ 1 I q diag(Λ q+1 , · · · , Λ l ). Moreover, the linear nonsingular transformation ξ = (J M 1 ) −1/2 T K 1 e −Λ 1 t η leads to the following system: and due to condition (21) (N 1 ) q = N 1 0 q , the equality holds. This means that subsystem: is separated. Using the Lyapunov function V = η τ q η q , one may see that: but because of rankN 1 = 0, the system (25) has at least one pair of purely imaginary eigenvalues. This implies the existence of complex eigenvalues for System (24) and resulting in the flutter instability of System (3). Now, the calculation procedure will be explained by the following case study.

Remark 1.
As follows from the formulation of the theorem, it gives the sufficient condition of instability. If the condition (22) is not satisfied, the stability of the system depends on the other components of matrix N. For instance, in Example 2, the roots of f (Λ) are: and if n 1 = 0, the system (26) is stable when n 2 2 < 1/2 (is spectrally stable when n 2 = ± √ 2/2).
Example 3. We consider one more example, which differs from the previous one. System (3) has a high dimension, and to avoid unnecessary conversions, we take M = I. Matrices K, N are: Even in this case, when matrices M, K are known, the dimension n = 5, and more importantly, the presence of five uncertain parameters may cause a serious obstacle for finding the stability conditions. The potential system has two pairs of multiple eigenvalues, and from a common viewpoint, it raises a query about the existence of a such matrix N (even with N 1), which would not bring the flutter instability. Again, the direct analysis of the characteristic polynomial is not so optimistic, i.e., the algorithm described in Section 3 still works, but the expressions calculated are very bulky. In order to omit the problem occurring, we try to employ Theorem 5.
The eigenvalues and eigenvectors of matrix K are as follows: As matrix K = K has two different multiple eigenvalues, we have to choose one of them to determine block-matrix K 1 . Let this matrix correspond to Λ 1 = 12, (actually, there is no difference; the choice of Λ = 24 leads to the same conditions), then with transformation matrix T and taking into account that q = 2, the equality (21) of the theorem imposes the following restrictions n 34 = −n 14 − n 24 , n 45 = 2n 14 + n 24 , and hence, the condition (22) takes the form: Thereby, if n 13 = 0, the system (MKN) is unstable (even with an infinitely small magnitude of circulatory force). However, if n 13 = 0, components n 14 , n 24 are responsible for stability or instability. As the system is decomposed, the characteristic polynomial is reduced to third order form: λ 3 + 7λ 2 + [16 + 3(n 14 + n 24 ) 2 ]λ + 12 + 51 4 n 2 14 + 12n 14 n 24 + 6n 2 24 .
Its coefficients are positive, and its discriminant: is positive in the domains shown in Figure 5. This result can be interpreted in the following way: when the potential system has multiple eigenvalues, the circulatory forces may contain "harsh" components, which break the stability no matter how small they are, and "soft" components are "dangerous" only if they are high enough.

Some General Remarks
Theorem 5 leads to the following remarkable consequence: if the linear part of the nonlinear dynamical system: is decomposed into unstable and stable (even asymptotically stable) subsystems and these subsystems are connected by nonlinear terms, then the equilibrium of the system is unstable. Moreover, the linearly stable component may bring more "discords" into nonlinear dynamics. Let us explain this by employing the following example: 1, 2, 3), N = n 1 0 1 −1 0 The linear part of System (27) is taken from Example 2 (with n 1 = 0), so it is "semi-stable", i.e., it consists of unstable and stable subsystems, and their eigenvalues are presented in Remark 1. To be able to perform numerical integration of the equations considered, we take n 1 = 0.1, n 2 = 0.3, then λ 1,2,3,4 = ±0.0499 ± 1.0012i, λ 5,6 = ±1.7029i, λ 7,8 = ±1.4491i. If the initial perturbations are sufficiently small, then for a period of time, the motion is regular. For System (27) with initial values x(0) = (0.001, 0.03, −0.001, 0.001) τ ,ẋ(0) = (0.01, 0.012, 0.001, 0.02) τ , this interval is approximately estimated as [0, 110]. The maximal Lyapunov exponent is not very big (≈0.05), and the behavior with respect to x 1 , x 2 is similar to linear instability (flutter). However, at some point, t crit ≈ 111, the system dynamics very rapidly becomes chaotic, and in a few seconds, the trajectory whirls away (leaves the visible limits; and the computation is interrupted). The clear sign of chaotic behavior is the fact that the very small change of even one initial value (x 1 (0) from 0.001 to −0.001) leads to the divergence of trajectories (see Figure 6). However, in applied research, Lyapunov's instability is not always a "desperate" case. Often, for reliable operation of the device, the growth of solutions is allowed in some limits. In this case, even the positive Lyapunov exponent does not necessarily lead to the failure of the system. In order to illustrate the latter statement, let us consider the following system: Its equilibrium x = 0,ẋ = 0 is unstable; however, numerical simulations show that the motions are bounded, if the initial perturbations are small enough (see Figures 7 and 8). Furthermore, the behavior of component x 3 is remarkable. Contrary to Example 2, the motion is "almost" asymptotically stable with respect to this coordinate and its velocity. Speaking more strictly, it is unstable, but tends to zero as time increases. Note that time to time, there appear short intervals where the amplitude increases during two-three oscillations more then 50 times; nevertheless, it very quickly returns to small values (Figure 9).

A Case Study: The Bounded Solutions of a 2-DoF System with Unstable Origin
Let D be some domain of R m containing the origin. Consider an ODE system: We assume that the function F is continuous, and the Cauchy problem for (29) has a unique solution x(t, t 0 , x 0 ) that continuously depends on the initial values.
For any t ∈ I = [t 0 , +∞), we consider the domains D 1 (t), D 2 (t) containing the origin x = 0, such that their boundaries ∂D 1 , ∂D 2 are continuous functions of time and ∂D l (t) ⊂ D 2 (t). The following statement holds. Lemma 1. [39]. Assume that there exists a continuously differentiable function U(t, x) : I × D 2 → R and a continuous function U 0 (t, x) ≤ U(t, x) such that: (A) For any t ≥ t 0 , the function U(t, x) is bounded for x ∈ ∂D 1 (t); (B) for any (t, x) ∈ I × D 2 (t)\D 1 (t), the time derivative of the function U(t, x) with respect to system (1) is non-positive; (C) for any pair t l , t 2 (t 2 > t 1 ) from I, the following inequality holds: Then, for any t 0 ∈ I, the premise x 0 ∈ D 1 (t 0 ) implies that x(t, t 0 , x) ∈ D 2 (t).
This statement may be applied for different purposes: to prove the stability (marginal or asymptotic), the boundedness of the solutions, the estimation of the rate of the asymptotic behavior, etc.
To illustrate it, let us consider the following system: Here, h, k 1 , k 2 , n are some positive dimensionless parameters. System (30) may be considered as the equations of motion of two identical non-linear oscillators connected due to the influence of circulatory force. This system has a single equilibrium x = 0,ẋ = 0, which is asymptotically stable if n = 0 and is unstable for any arbitrary small non-zero value of n. The eigenvalues of the linearized system are located at the tops of the square in the complex plane; their values are: Thus, the maximal Lyapunov exponent is equal to n/ √ 2R, and the origin is exponentially unstable. However, if the ratio n/k 1 is small enough, then all trajectories incipient in the vicinity of the origin are bounded. To prove this fact, we shall use the lemma.
First of all, we realign the variables, parameters, and time according to formulas: This allows formally putting the values of k 1 and k 2 in Equation (30) equal to one (for convenience, the symbol "∼" is omitted below). As the system considered is autonomous, we shall seek the polynomial positive definite auxiliary function as the sum of the quadratic part and the form of fourth order: U(x,ẋ) = U (2) + U (4) . Obviously, it is impossible to achieve the negative definiteness ofU (2) (the presence of a positive Lyapunov exponent); thus, we aim for the negativeness ofU (4) to satisfy Condition (B) of the lemma. The function U is taken in the following form: and its time derivative with respect to System (30) is: The upper bound for derivativeU(x,ẋ) may be derived in the form: Here, r j = x 2 j +ẋ 2 j (j = 1, 2), p 1 ≈ 1.4, and the value of p 2 (h) depends on parameter h and is defined as: The surface p = p 2 (h, ξ) is presented in Figure 10. From the other side, the function U satisfies the following inequalities: In the vicinity of the origin of the first quadrant of the plane Or 1 r 2 , the boundary of the set U(r 1 , r 2 ) = 0 is formed by two branches (solid lines in Figure 11a). The smaller one bounds the domain U > 0. Thus, we should choose the domain D 1 in such a way that ∂D 1 belongs to domain U ≤ 0. It can be taken as U 1 (r 1 , r 2 ) ≤ α(h, n), where α is some constant number. Then, the boundary of D 2 may be given by condition U 0 = α + ε, where ε is an arbitrary small positive number. Therefore, all conditions of the lemma are fulfilled, and all trajectories starting in D 1 are limited by the threshold of ∂D 2 . For values h = 4, n = 0.1, these domains are presented in Figure 11. This conclusion is confirmed by numerical integration of System (30). For different starting points of the phase space, the trajectories go to the limit cycle ( Figure 12). In fact, the estimation shown in Figure 11 is accurate enough.

Concluding Remarks
In this paper, we discussed some stability problems for circulatory dynamical systems. Although existing analytical methods allow us to solve the problem (obtaining stability conditions in the closed form) in the case of high order systems with uncertain parameters, these conditions may turn out to be too cumbersome for analysis. Therefore, it is advisable to use optimized algorithms from the point of view of avoiding "unnecessary" computational procedures, which allows reducing the computation time and, possibly, obtaining stability conditions in a simpler form. Some of these procedures were discussed in Section 2 and tested in Section 3. In Section 4, a case of multiple eigenvalues was studied. It is interesting from a theoretical viewpoint as well, as this case often occurs in different applications. A sufficient condition of instability was formulated. In Section 5, some interesting examples of nonlinear dynamics were considered. Then, it was demonstrated that the flutter instability of the linearized system may lead to various results, namely the motion may become chaotic with overwhelming growth of the amplitude and may be bounded. Actually, it may be even asymptotically stable with respect to part of the variables. However, the rigorous study of this question is based on the use of Lyapunov functions for a nonlinear multi-parameter system, which leads to extremely complicated transformations and requires a separate study.
The future research will consider the problem of free flexural vibrations and dynamic loads in an elastic rod with stepwise-changing stiffness, as well as the application of the stated results to the hemivariational principles such as damage [40,41]. +4a 1 a 3 a 2