Existence of analytical solution, stability assessment and periodic response of vibrating systems with time varying parameters

The paper is focused on the solution of a vibrating system with one-degree-of-freedom with the objective to deal with the method for periodical response calculation (if exists) reminding Harmonic Balance Method of linear systems having time dependent parameters of mass, damping and stiffness under arbitrary periodical excitation. As a starting point of the investigation, a periodic Green’s function (PGF) construction of the stationary part of the original differential equation is used. The PGF then enables a transformation of the differential equation to the integro-differential one whose analytical solution is given in this paper. Such solution exists only in the case that the investigated system is stable and can be expressed in exact form. The second goal of the paper is to assess the stability and solution existence. For this purpose, a methodology of (in)stable parametric domain border determination is developed. c © 2020 University of West Bohemia. All rights reserved.


Introduction
Periodic systems have been subjected to the continuing interest in many branches of engineering. Mainly the stability of these systems has been intensively investigated and, therefore, a lot of works covering this topic can be found in the literature. Hill [6] seems to be one of the pioneers who laid the first mathematical foundations for the stability assessment of parametric systems. Since then various techniques have been developed. Methods used for the stability assessment of periodic systems can be divided into several categories.
First of all, there are perturbation methods [7,15,16]. According to them, the solution is given by the first several terms of an asymptotic expansion. The standard is usually not to use more than two terms. These approaches are known as first-order and second-order perturbation techniques. It should be mentioned in the connection with a stability assessment that the perturbation methods are primarily suitable for systems containing only small fluctuation parameters. Further, the various approaches of the perturbation techniques can be shown in the relation to the assessment of stability. They include the methods of Wentzel-Kramers-Brillouin (WKB) [14] or Poincare-Lindstedt [5].
The second category of assessment methods includes those based on the Lyapunov theory [7,18]. Lyapunov was the first to use linearisation of nonlinear systems for stability assessment near a point of equilibrium (work point). A linearisation using the Jacobian matrix corresponds to the expansion of the function on the right hand side to a power series and all no-linear terms are omitted. Whether the system is stable or not then depends on the eigenvalue properties of the Jacobian matrix. This technique is known as the First method of Lyapunov. The Second method of Lyapunov uses a Lyapunov function having an analogy to the potential function of physical systems.
The last and quite broad category of methods includes those based on the Floquet theory [3]. This theory allows to study the stability of linear time periodic systems. For a nonlinear periodic system, the stability of the solution can be analysed by means of a linearised mathematical model (e.g., see [7,11]). However, sometimes it gives less accurate prediction of stability boundaries when individual steps of the solution are approximated by linear systems. In the context of the Floquet theory, further methods were developed such as the State Transition Matrix (STM) methods [4,16] and the Infinite Determinant (ID) methods [7,21]. The STM methods are essentially numerical in nature and are applicable without restrictions to the order of the system and the magnitude of the parameters. This is one of the main reasons why these methods are widely used for the calculation of characteristic exponents. In contrast, the ID methods are more commonly applied to second-order systems with small periodic parameters. The use of the ID method for systems of any order and with arbitrary coefficients is demonstrated in [1]. However, in these cases, the calculation is rather cumbersome because of the need to compute a quite large determinant. As alternative to ID methods, the method of multiple scales and the method of strained parameters can be found in [16].
Besides stability assessment, the second important task in the investigation of ordinary differential equations (ODEs) with time periodic coefficients is to find their periodic solutions (if exist). From a practical point of view, one of the most important is the second-order ODE known as Hill's equation, which was first studied by Hill in 1886. Whittaker and Watson in [21] provided the periodic solutions involving Hill's ID. Further, more general results on the Hill's homogeneous equation can be found in [10]. Shadman and Mehri [19] investigated the existence of periodic solutions of forced Hill's equation and extended the problem to the case of a non-homogeneous matrix valued Hill's equation. Additionally, there are several works such as [16] that use the Harmonic Balance Method (HBM) to find an approximate solution of the Hill's equation. In special cases where Hill's equation has a specific form, this one becomes the Mathieu's equation [12], the Whittaker-Hill's equation [20] or the Meissner's equation [13]. Among them, the best known is probably the Mathieu's equation which in the past was studied by many researchers including MacLachlan [9], Newland [17] or Jordan and Smith [7]. Kourdis and Vakakis [8] solved the equation of linear oscillator under parametric excitation of general type and used an analytical approach based on amplitude-phase decomposition of the response.
The current work builds on and significantly expands our previous study [2], where we presented an analytical solution and stability assessment of a one-degree-of-freedom (1 DOF) linear vibrating system with periodic stiffness that was excited by a periodic force. In this paper, the equation of motion with other periodically varying parameters under periodic excitation is studied. To find the steady-state solution, presented approach uses the periodic Green's function (PGF) as a response to the Dirac chain. The application of the dynamical compliance method leads to an integro-differential equation with degenerated kernel. This approach may resemble the use of HBM to calculate the periodic solution of linear differential equations with varying parameters. The main difference between the classical HBM and the presented method is the introduction of a system matrix depending on the fluctuation matrices of stiffness, mass and damping. The sign of a 2T -periodic characteristic matrix determinant decides on (in)stability and existence of the periodical steady-state solution. The ability of the proposed methodology is illustrated through examples.

Periodic solution to equation of motion
Let us assume that the behaviour of a system is described by the equation of motion where some parameters and excitation are time periodic. It means that the following conditions of periodicity are satisfied: All these functions can be expressed in the form of the Fourier series while using the compact notation e n (t) = e inωt with ω = 2π/T and where ω is the basic angular frequency and T is the period corresponding to the frequency ω. The parameters μ m , μ b , μ k in (1) represent the measures of mass, damping and stiffness modulation, respectively. The remaining symbols have the following meaning: m s , b s , k s are stationary mass, damping and stiffness parameters, respectively. Further, let us assume that the number N is sufficiently large for a good approximation of the periodic functions. PGF should be taken as a starting point for finding a solution to the equation of motion (1). This function can be obtained as the periodic response of the stationary part of the system to a Dirac chain excitation and has the following form [2] The solution of (1) can be written in the form of convolution integrals (excitation consists of the parametric and external parts) Substituting (2)  into (5) and considering e −n (s)e j (s) = e j−n (s), the function q(t) can be briefly written Let us introduce matrices Equation (7) represents an integro-differential equation with degenerated kernel for the function q(t). Using (8), this integro-differential equation can be rewritten into the form where we used the condition of orthogonality The matrix L represents the dynamic compliance of the stationary part of the system that corresponds to individual harmonics. Let us introduce the following notations so that (9) can be rewritten into the form The first and the second time derivatives of (15) arė where The derivatives of the function q(t) can be expressed aṡ To obtain α, β, γ, we pre-multiply (15), (18), (19) by , respectively, and integrate from 0 to T . Taking into consideration (14), the equations can be further rewritten as The matrix products A mÎ , A bÎ , A kÎ can be replaced as follows: which can be easily proven by direct substitution. All matrices . .
are Hermitean, i.e., which follows from the fact that x −j = x j for j = 0, x 0 ∈ R. The bar notation means complex conjugated expression. Matrices H m , H b and H k express fluctuation of mass, damping and stiffness, respectively. Taking (22) and (20), the system of equations (20) can be converted into the matrix form where Let us denote A as a T -periodic system matrix and vector b as a T -periodic system excitation.
Vector y contains coefficients of linear combination of approximation functions. Parameter μ can be freely chosen from the set μ ∈ {μ k , μ b , μ m }, while the values of c α , c β , c γ have to be changed accordingly, for example, The solution of (25) has the form where I is an identity matrix of the same type as A. The matrix in round brackets can be called as a T -periodic characteristic matrix. Taking (28), the solution (15) can be rewritten as and subsequently as where The use of (29) to the final solution of (1) then yields The accuracy of this analytical periodic solution depends on the number of N harmonics taken into account. Let us denote the solution approximated by 2N + 1 terms of the Fourier series as q N (t). Thereafter, we assume the solution q N (t) → q(t) if the condition is satisfied, where ε is a small positive number.

Existence of periodic solution and system stability assessment
If the aforementioned system is unstable, then the periodic solution (33) does not exist and, therefore, we focus mainly on the stability assessment. It is necessary to take into account resonant stage with parametric excitation having half frequency ω * = ω/2, see [2]. It corresponds to the 2T -periodic solution, which has to be identical with the T -periodic solution. It will be proven bellow.
The 2T -periodic solution can be expressed in the analogical form to (30). Then where Remark 1. The relations given above can be derived analogously to those in Section 2.1. We start by finding the 2T -periodic PGF defined as the response of the stationary part of the investigated system to a Dirac chain excitation with a 2T period. It follows that We make sure that the following relations hold where matrix J is defined in (A.2) and matrix H * x is given as Similarly, H x , see (24), as well as H * x is Hermitean for x ∈ {m, b, k}. The subscripts correspond to the original angular frequency ω. Applying the same procedure as in the T -periodic system, we can come to the equation for the calculation of the Fourier-coefficients vector. The form of this equation is analogous to (25), i.e., where where The total 2T -periodic solution has the form (35) and can be written in a similar form as (31) or (33), i.e., where and I is an identity matrix of the same type as the 2T -periodic system matrix A * . The matrix I − μA * can be called a 2T -periodic characteristic matrix.

Lemma 2. Let q(t) in (33) be the T -periodic solution to the equation of motion (1) in steadystate. Then, this solution can be also written in the form
using a T -periodic reduced system matrix Proofs of Lemmas 1 and 2 are provided in Appendixes B and C, respectively. Now, let us turn our attention to the investigation of the (in)stability border. It is evident from (48) that some solution q * (t) exists only if the matrix I − μA * is not singular. Otherwise, 1/μ is equal to one of the eigenvalues of the matrix A * . As shown below, the determinant of A * is real and the eigenvalues of A * are either real or complex conjugate pairs. Consequently, the real eigenvalues correspond to the changes of determinant sign and determine the (in)stability border in the parametric plane ω and μ.
First, it can be easily shown that the relations are valid because of the orthogonality of matrices P and Q ( P T P = I and Q T Q = I), see Appendix A. For this reason, A * , A and A are similar matrices with identical spectra. The symbol Σ denotes the set of eigenvalues of the relevant matrix. The proof of (53) follows from the relations between the matrix A * and the matrices A and A A = P T A * P and A = Q T P T A * PQ .
Lemma 3. The set of nonzero eigenvalues of the matrices A * , A and A is identical to the set of eigenvalues of a 2T -periodic reduced system matrix From this lemma follows the fact that in order to obtain the borders of system (in)stability, it is possible to solve the set of eigenvalues of matrix A * red instead of that of matrix A * .
Lemma 4. Eigenvalues of the 2T -periodic reduced system matrix A * red are either real or complex conjugate pairs. For this reason, the determinant of this matrix is a real value because the matrix determinant can be expressed as a product of all eigenvalues.
Proofs of Lemmas 3 and 4 are given in Appendixes D and E, respectively.

Results and discussion
The verification of the relations derived in previous sections will now be demonstrated on a simple problem. For this purpose, let us consider the system At this point, it is appropriate to transform this equation using the dimensionless time τ = Ωt into the form where and The derivations of the variable q with respect to dimensionless time can be denoted by () = d/dτ and the coefficients c α , c β , c γ are used in accordance with (28).

Stability problem
The stability assessment is performed for homogeneous equation (57). As mentioned earlier, the determination of the (in)stability boundaries corresponds to the eigenvalue problem of matrix A * red , see Lemma 3. This matrix is defined for a damped system (D = 0.01) by (55) and has a reduced form c α L * H * k − ω 2 c γ L * H * m N * 2 corresponding to the conservative system (D = 0). According to Lemma 4, the eigenvalues of matrix A * red are either real or come in complex conjugate pairs. However, only the real eigenvalues define the (in)stability borders because the investigated system has real periodic parameters. Two different value sets of parameters c α , c β , c γ are presented: c α = c β = c γ = 1 in the first case and c α = c β = 1.1, c γ = 1 in the second one. When analysing the undamped system, the parameter c β is omitted in both aforementioned cases. Numerical experiments showed that for N > 10, the number N has a negligible effect on the accuracy of values |μ| < 1. Therefore, all computations in this paper are performed for (in)stability borders with parameter settings N = 10.
The results of stability assessment are shown in Fig. 1. It can be noted that a relatively small change in parameters c α , c β leads to stability regions of considerably different shapes. This behaviour can be seen in the undamped ( Fig. 1(a) and (c)) as well as in the damped ( Fig. 1(b) and (d)) system. The most apparent differences are particularly obvious for the largest value of 1/η 2 (the smallest value of the parameter ω). The Floquet theory was applied to verify the correctness of the derived (in)stability borders. The characteristic exponents were computed for the combinations of 1/η 2 and μ by using the Runge-Kutta integration method. The dots in Fig. 1 represent points of stability. As the reader can see, the obtained results are in a good agreement. Considering only periodic stiffness in (57) (known as the Mathieu equation), the stability regions become symmetric as shown, e.g., in [2]. All numerical tests have shown that the determinant of matrix I − μA * red is positive in the regions of stability and negative in the regions of instability. In the later case, a periodic solution for steady state does not exist. It could be used as a criterion for determining the stability of linear systems.

System response
The validity of (51) for the steady-state response is demonstrated for two types of excitation. Because the stability of the system given by (57) has already been proven, the response can be also solved. The two periodic excitations were chosen in the form of a linear cosine function which can be defined as and in the form of a saw function The functions f 1 (t) and f 2 (t) are depicted in Fig. 2(a) and (b), respectively. Due to the Tperiodicity of both functions, they can be expressed in the form of Fourier series. The original functions f 1 (t) and f 2 (t) can be transformed into the form and The parameter f S is assumed to be f S = k S q S in both cases of excitations, where q S is a static displacement corresponding to the f S value. An appropriate number of Fourier series terms N 1 and N 2 to describe the functions f 1 (τ ) and f 2 (τ ), respectively, was taken with respect to numerical experiments and to the negligible effect of higher terms on the total response. The values N 1 = 25 and N 2 = 50 are used for the subsequent calculations. Because only the periodic steady state response is required in this work, the sought solution corresponds only to the particular solution of the equation of motion (1) and is given by a convolute integral, see equation (5). Then the solution consists of the solution with no excitation and with null initial conditions. For this reason, the verification of numerical results by the Runge-Kutta integration method respects homogeneous initial conditions in the following examples. All analyses are then performed with respect to (57).
Recall that the presented solution q(t) given by (51) makes sense only in the stable region. If the parameters are chosen from the unstable domain, e.g., D = 0.01, η −2 = 2, μ = 0.7, c α = c β = 1.1 and c γ = 1.0 (see Fig. 1), then a periodic steady state response does not exist. This is shown in Fig. 3, where results are obtained by the Runge-Kutta integration method. It is also well known that any exitation of linear systems do not affect the stability assessment.
Results obtained from the analytical solution and the Runge-Kutta continuation are compared in

Conclusions
In this paper, a methodology for analytical solution and stability assessment of differential equations with periodically varying parameters was presented. The differential equations can describe behaviour of mechanical vibrating systems containing time periodically varying mass, damping and stiffness parameters. The main idea of the presented parametric system investigation was the transformation of the differential equation with time dependent coefficients to the integro-differential equation, whose response was expressed by means of convolution of periodic Green's function and a sum of the external and parametric excitation. The paper made several statements with corresponding proofs. The first statement was related to the real valued determinant of the 2T -periodic characteristic matrix I − μA * , whose eigenvalues are either real values or complex conjugate pairs. In this work, it was proven that non-zero eigenvalues of the 2T -periodic system matrix A * can be calculated from the so-called 2T -periodic reduced system matrix A * red , whose order is three times smaller than that of the original matrix A * . The proposed approach for the determination of stability domain and its border was demonstrated on examples and the obtained results were checked by the Floquet method showing very good agreement. Furthermore, the performed numerical experiments have shown the same results as in [2], namely, that a positive sign of the determinant of the 2T -periodic characteristic matrix corresponds to the solution in the stable region. The periodical response obtained by presented analytical approach in stable parameter domain (do not confuse A red and A * red ) showed also good agreement with the results of the Runge-Kutta continuation. The main benefits and advantages of the presented approach (PA) can be summarised in the following points: • Unlike the classical harmonic balance method, the PA enables to perform stability assessment based on the knowledge of fluctuation matrices of mass, damping, and stiffness and the real valued determinant of the 2T -periodic characteristic matrix.
• The PA enables to exactly determine the border of (in)stability with arbitrary accuracy. This approach is not limited by the magnitude of parameters μ m , μ b and μ k .
• In case of stability, the presented methodology enables to find an analytical solution of the system response in the steady state.
• The PA can be used in sensitivity analysis required, e.g., for the solution of stochastic vibration or in an optimisation process.
• The PA can be extended for systems with n < ∞ DOF.
The Fourier-coefficients vectors can be expressed in changed order in such a way that upper subvector is formed by coefficients corresponding to the T -periodic solution and the remaining coefficients are assembled to lower subvector according to scheme (δ * ∈ {α * , β * , γ * }) Each vector i n (n = 1, . . . , 4N + 1) in (A.2) corresponds to the n-th column of the identity matrix I ∈ R 4N +1,4N +1 . Equations (A.1) can be also given in the forms: Due to the orthogonality of the permutation matrix P, the matrices J and U are orthogonal, and generate vector subspaces P 1 and P 2 , respectively. Both subspaces represent orthogonal supplements to each other, which can be written as Let us introduce an orthogonal matrix and multiply the right hand side of the equation y * = μA * y * + b * by the matrix P T . The use of (A.1) then yields y = P T y * = μ P T A * P y + P T b * , Let us put the following notations Then, equation (A.9) can be rewritten as In order to re-express the matrix A in the other forms, the following submatrices are analysed: The first of them can be so arranged as to be block diagonal because of J T L * H * k U = 0, as the result of orthogonality (A.7), and J T L * J = L and J T H * k J = H k . Similarly, it can be shown that By comparing matrices in (A.11) and (26), and by taking into account (A.14)-(A.16), it is possible to write and where notation of the submatrices A ij (i, j = 1, 2, 3) follows from (26). The vector b in (A.12) can be rewritten as and substitute (A.19) to (A.13). Types of the identity matrices I in the matrix Q correspond to dimensions of the subvectors α U , α L etc. Pre-multiplying (A.19) by the matrix Q T , we obtain due to its orthogonality (i) The terms μE * (t)y * and μE(t)y are compared. According to (A.1), we can write Further, with the help of (A.9) and (49), one obtains where while the vector y in (A.10) takes the form with respect to (A.27) and (A.29). Multiplying E(t) and y, and taking into consideration that e U (t) ≡ e(t), the following equalities hold: It means that the first terms on both sides of (B.1) are equal to each other. The second term on the left hand side of (B.1) can be rewritten into the form The validity of U T L * J = 0 follows from the orthogonality of matrices J and U, see (A.7), and from the fact that the matrix L * is the diagonal one.
These results indicate that the proof of identity (B.1) is completed and the identity (50) is also proven.
where the matrix A * is defined in (45). If the transformation is used, and the substitution (D.2) into (D.1) and the pre-multiplication (D.1) by the matrix (W * ) −1 are performed, the original eigenvalue problem can be reformulated as The set of eigenvalues of block triangular matrix consists of eigenvalue subsets corresponding to the individual diagonal block submatrices. However, this means that the spectral matrix of the matrix A * has the form The proof is hereby completed.

Appendix E Proof of Lemma 4
The independence of numbers c α , c β and c γ is assumed. This condition is sufficient for eigenvalues in the individual terms of (55) to be either real or complex conjugate pairs because where Λ R , Λ S and Λ T are the spectral matrices (diagonal matrices having eigenvalues on the diagonal) corresponding to the matrices R, S and T, respectively, which are analysed separately. where the matrixÎ is defined in (21), we can gradually write From (E.5) follows the fact that the spectrum of the matrix L * H * k is identical with the spectrum of its complex conjugated matrix L * H * k . For this reason, their eigenvalues have to be either real or complex conjugate pairs. Assuming that c α , c β , c γ are real constants, we can extend this statement to the matrix R = c α L * H * k .
ii) Properties of the matrix S, see (E.2) 2 : Taking the relationŝ (E. 6) and the relations in (E.4), we want to prove the equality of the first two terms in the equation because N * is real valued. By simple arrangements, we can come to Comparing the last terms in (E.7) and (E.8), we can say that the equality of the first two terms in (E.7) is proven. This equality means also that the spectrum of the matrix iωL * H * b N * and the matrix S consists of eigenvalues which are either real or complex conjugate pairs.
iii) Properties of the matrix T, see (E.2) 3 : ConsideringÎ we can come to When the spectrum of the complex conjugated matrix is identical with that one of the original matrix, we can say that it consists of either real or complex conjugate pairs of eigenvalues.
The proof of Lemma 4 is finished.