# Lyapunov definition and stability of regular or chaotic vibration modes in systems with several equilibrium positions

Yuri Mikhlin, T.V. Shmatko, Gayana Manucharyan

► To cite this version:

Yuri Mikhlin, T.V. Shmatko, Gayana Manucharyan. Lyapunov definition and stability of regular or chaotic vibration modes in systems with several equilibrium positions. Computers & Structures, 2004, 82 (31–32), pp.2733-2742. 10.1016/j.compstruc.2004.03.082. hal-01348406

HAL Id: hal-01348406

https://hal.science/hal-01348406

Submitted on 23 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LYAPUNOV DEFINITION AND STABILITY OF REGULAR OR CHAOTIC
VIBRATION MODES IN SYSTEMS WITH SEVERAL EQUILIBRIUM
POSITIONS

YU.V.MIKHLIN

*Department of Applied Mathematics, National Technical University, 21 Frunze str.,*
*Kharkov, 61002, Ukraine*

E-mail: muv@kpi.kharkov.ua

T.V.SHMATKO

*Department of Higher Mathematics, National Technical University, 21 Frunze str.,*
*Kharkov, 61002, Ukraine*

E-mail: ktv_ua@yahoo.com

G.V.MANUCHARYAN

*Department of Applied Mathematics, National Technical University, 21 Frunze str.,*
*Kharkov, 61002, Ukraine*

E-mail: gayane@kpi.kharkov.ua

This paper deals with forced vibrations of two-DOF systems with more
than one equilibrium positions. Such systems may be obtained by dig-
itization of elastic post-buckling systems. A vibration mode, which is
periodic at small force amplitudes and becomes chaotic as the force am-
plitudes are slowly increased, is selected. It is possible to formulate and
solve the problem of stability of a periodic or chaotic vibration mode
in a space with greater dimension using the classical Lyapunov sta-
bility definition and some calculating procedures. Instability of phase
trajectories is used as a criterion of the chaotic behavior in dynamical
systems. Trajectories with very close initial values are compared. Use
of the Lyapunov stability definition shows mutual stability/instability
of the trajectories. Calculations permit to observe an appearance and
enlargement of the chaotic behavior regions. Specific results are ob-
tained for the nonautonomous Duffing equation and pendulum system.

1 Introduction

Models under consideration may be obtained by digitization of some elastic systems that
have lost stability under a constant compressive force. If displacements of the nonlinear
elastic system are approximated by a single harmonic of the Fourier series expansion for
space coordinates, a system having a single degree of freedom is obtained. Behavior of the
model described by the nonautonomous Duffing equation was examined in papers [1,2 et al.].
At small force amplitudes, periodic motions close to one of two stable equilibrium positions
are observed but irregular, apparently chaotic motions begin as the force amplitudes are
slowly increased.

If the first two harmonics of the Fourier series for space coordinates are used, one obtains
a set of two-second order ODEs, coupled in nonlinear terms only. In this case an energy
”pumping” from one vibration mode to another mode is possible. Thus one can formulate a problem of the stability of periodic or chaotic vibration mode in a space with a greater dimension. There are many publications on the general and special problems of stability of periodical solutions, e.g., fundamental works by Lyapunov [3]; Ince [4]; Andronov, Vitt and Khaikin [5]; Minorsky [6] etc. But a problem of stability of chaotic modes in a space with a greater dimension has no analytical solutions in a general case.

We make use of the classical Lyapunov stability criterion for a case when initial variations are not arbitrary small that is they are limited below (Section 3). As a result we obtain some test suitable for a computer realization. It is discussed a connection of this test with an approach used by Schiehlen [7], and with so-called Lyapunov exponents (Section 3). Calculations, corresponding to the proposed test, are realized as long as boundaries of stability and instability regions in the space of system parameters are variable. Stability of regular or chaotic modes of nonlinear beams and shells (Section 4) are considered. Note that first results obtained for this problem were published in papers [8,9].

The chaotic behavior of nonlinear dynamical systems is an object of recent and current investigations. The first ideas on the chaotic behavior of dynamical systems can be found in works by Poincare' [10] and Birkhoff [11 et al.].

Criteria for the onset of chaos beginning in dynamical systems are very diverse. A detailed survey of various such criteria of the chaos beginning was published by Moon [2]. The most important analytical or numerical-analytical criterions of the chaos beginning: criteria are based on the configuration of the homoclinic trajectories, period doubling, overlapping of resonances, analysis of spectra, calculation of Lyapunov exponents, “Smale horseshoe” et al [12-15]. All trajectories contained in the “strange attractor” region are unstable, that is trajectories which are close to each other at some point of time then diverge exponentially [2,16,17].

We suggest here an approach to determine the onset of chaos based on the Lyapunov definition and the proposed calculating test (Section 5). Mutual instability of phase trajectories is used as criterion of chaotic behavior in dynamical systems. One compares trajectories that are initially very close. The proposed test shows the mutual instability of the trajectories. Calculations permit to observe a process of appearance and enlargement of the chaotic behavior regions if some selected parameters of the dynamical system under consideration are changing. Specific results were obtained for the nonautonomous Duffing equation and pendulum system when an amplitude of the external periodic action increases.

## 2 Principal Model

Consider the following system that can be obtained by discretization of equations of nonlinear dynamics of some elastic systems:

$$ 
\begin{aligned}
\ddot{y}_1 + \delta \dot{y}_1 - \alpha y_1 + \beta y_1^3 + c y_1 y_2^2 &= f \cos \omega t, \\
\ddot{y}_2 + \delta \dot{y}_2 + ay_2 + by_2^3 + c y_2 y_1^2 &= 0,
\end{aligned}
\quad (1) $$

where $y_1(t)$ and $y_2(t)$ are unknown functions; $\delta$ is the coefficient determining friction; all coefficients are positive (excepting may be the coefficient $\alpha$).

The system (1) can be obtained in the following problem. Consider nonlinear bending vibrations of a beam with a length equal to $l$. Within framework of the Kirchhoff hypothesis
the vibrations are described by the following equation [18]:

$$ \mu \frac{\partial^2 w}{\partial t^2} + EI \frac{\partial^4 w}{\partial x^4} - \frac{ES}{2l} \int_0^l \left( \left( \frac{\partial w}{\partial x} \right)^2 dx \right) \frac{\partial^2 w}{\partial x^2} + \Gamma \frac{\partial^2 w}{\partial x^2} = F(x, t), \quad (2) $$

where $x$ is the space coordinate; $w(x,t)$ is the bending deflection; $\mu$ is the density of the beam material; $E$ is the Young's modulus; $S$ is the area of the transversal cross-section; $I$ is the corresponding moment of inertia; $\Gamma$ is the axial compressive force; friction will be introduced later; $F(x,t)$ is the distributed external effect.

Let $w = w_{xx} = 0$ if $x = 0$ and $x = l$. Let the function $F(x,t)$ be presented as $F(x,t) = f_1 \cos \omega t \sin(\frac{\pi x}{l})$. Representing the deflection as truncated Fourier series (by sines) and making use of the Bubnov-Galerkin procedure for (2), one comes to a system of nonlinear ODEs connecting coefficients of the series. Taking into account two harmonics of the Fourier series expansion,

$$ w = y_1(t) \sin \frac{\pi x}{l} + y_2(t) \sin \frac{n\pi x}{l}, \quad (3) $$

one obtains the system (1) for determination of $y_1(t)$ and $y_2(t)$, where $\delta$ is the coefficient determining friction;

$$ \alpha = -\frac{1}{\mu} \left[ EI \left(\frac{\pi}{l}\right)^4 - \left(\frac{\pi}{l}\right)^2 \Gamma \right]; \quad \beta = \frac{ES}{4\mu l} \left(\frac{\pi}{l}\right)^4; \quad c = \frac{ES}{4\mu l} \left(\frac{\pi}{l}\right)^4 n^2; $$

$$ a = \frac{1}{\mu} \left[ EI \left( \frac{n\pi}{l} \right)^4 - \left( \frac{n\pi}{l} \right)^2 \Gamma \right]; \quad b = \frac{ES}{4\mu l} \left( \frac{n\pi}{l} \right)^4; \quad f = \frac{1}{\mu} f_1. $$

In the post-buckling dynamics we have $\Gamma > \frac{\pi^2 EI}{l^2}$, and the coefficient $\alpha > 0$ too.

If $y_2 = 0$, one obtains a nonautonomous Duffing equation with two potential wells (for a case $\alpha > 0$). It was shown [1,2] that as the force amplitude $f$ increases, the behavior of the system changes. A sequential doubling of the vibration period results, and beginning with some values of $f$ it is possible to observe irregular, apparently chaotic motions (for example, if $\alpha = 10, \beta = 100, \delta = 1, \omega = 3.76$, the chaotic vibrations appear if the force amplitude $f$ is a little greater than 1). Additional analyses (spectral, and also based on Poincare' mapping) show that a region of the phase space exhibits the properties of a "strange attractor".

Equations which are similar to the system (1) can be obtained also in problem of pre-buckling or post-buckling nonlinear dynamics of plates or shells.

## 3 Lyapunov stability definition and its computer realization

One consider a problem of stability of the periodic or chaotic vibration mode $y_2 = 0$ of the system (1). Instability of the mode means energy "pumping" from one harmonic of the Fourier series to another.

The variables $y_2, \dot{y}_2$ emerge in the form of variations (or perturbations). Assume first that values of the variations are significantly smaller than values of the variable $y_1$ (in the
region of stability of the mode $y_2 = 0$). This corresponds to a conventional formulation of the problem of local stability in configuration space.

Consider the well-known Lyapunov stability criterion [3,6,19] stating that: the solution $y = 0$ is stable if for all $\varepsilon > 0$ there exists a $\delta > 0$ such that for all $y_0 \in N_\delta(0)$ and $t \ge 0$ we have $y(t) \in N_\varepsilon(0)$. Here $N_\delta^{(0)}$ and $N_\varepsilon^{(0)}$ mean $\delta$- and $\varepsilon$- neighborhoods of the mode $y = 0$. For example, $N_\gamma = \{y \in \mathbb{R}^n \mid |y| \le \gamma\}$ (the neighborhoods could be calculated differently).

Introduce a relation between the quantity $\varepsilon$ and the initial value of the variable $y$. Let

$$ \varepsilon = \rho |y_0| \le \rho \delta \quad (\rho = \text{const}) \qquad (4) $$

This means that a value of $\delta$ is not arbitrarily small because $\delta \ge \frac{\varepsilon}{\rho}$. Note that this assumption does not contradict to the Lyapunov definition meaning because in this definition the initial values of $\delta$ can not tend to zero. One has from (4) that $\rho \ge \frac{\varepsilon}{\delta}$, that is $\rho$ is a high limit of the fraction $\frac{\varepsilon}{\delta}$. Besides, one obtains from the Lyapunov stability criterion taking into account the inequality (4) that $|y(t)| \le \rho |y_0|$.

One has from the preceding the following stability test for the system (1):

Instability of the mode $y_2 = 0$ is established if

$$ \max_{t \in [0, T]} |y_2(t)| \ge \rho |y_2(0)| \qquad (5) $$

The inequality (5) gives us a possibility to use a computer in the next stability analysis.

It is necessary to choose values of $\rho$ and $T$. Here a value $\rho^{-1}$ is a measure of smallness of initial variations with respect to maximum admissible variations. An increase of the value $\rho$ means that feasible initial values of variation decrease. There is some arbitrariness in a choice of the value $\rho$. In engineering one says that some quantity $\alpha$ is greater than another quantity $\beta$ “in order of magnitude” if $\frac{\alpha}{\beta} \approx 10$. We can choose, for example, $\rho \approx 10$ which constitutes as an “engineering test”. The discussion on the choice of the constant $\rho$ is continued in the Subsection 3.2.

Calculations are conducted as long as boundaries of stability and instability regions (in a chosen scale) on the system parameter space are variable. This is the principal criterion for the choice of the calculation time $T$. Note that all calculations are realized by using the standard Runge-Kutta procedure. The calculations are made at points on some chosen mesh of the system parameter space. It is clear that if the mesh widths decrease and the number of mesh points increases infinitely, the interval of time $T$ tends to infinity.

We now discuss the dependence (or independence) of the stability analysis on the variations initial conditions. The linear stability results are not dependent on initial conditions, this is a property of linear variational equations. But it is known [20,21] that additional nonlinear instability regions (obtained if we take into account nonlinear terms) have a smaller dimension in parameter space than instability regions obtained by the linearized stability analysis. Examples verify that the stability analysis is independent of initial variations if the initial variations are small.

## 3.1 Schiehlen's analysis of stability

In 1993 W. Schiehlen [7] proposed an approach for the stability assessment of nonlinear oscillations “from an engineering point of view”. His approach was based on the Lyapunov
stability definition. Schiehlen introduced the “inverse stability measure”

$$IS1(x_0, t_0) = \begin{cases} \frac{\|x(t)\|_{\infty}}{\|x_0\|} & \text{for } x_0 \neq 0 \\ 1 & \text{for } x_0 = 0 \end{cases} \qquad (6)$$

where

$$\|x(t)\| = \max_{1 \le i \le n} |x_i(t)|, \\ \|x(t)\|_{\infty} = \lim_{T \to \infty} \|x(t)\|_{T} = \lim_{T \to \infty} \max_{t \in [t_0, T]} \|x(t)\|.$$

The measure $IS1(x_0, t_0)$ “characterizes the ratio between a given initial state and the corresponding maximal displacement of the trajectory”. Here $x(t)$ is the time-dependent n-state vector.

The inverse stability measure

$$IS2(r, t_0) = \max_{x_0 \in \{x: \|x\|=r\}} IS1(x_0, t_0)$$

“characterizes the maximal displacement of all trajectories starting out of the initial conditions’ subspace which is by definition a hyper-cube with respect to the equilibrium point $x=0$”.

These measures were used in a numerical calculations where the integration interval is limited, $0 \le t \le T$. So, the measure (4) was replaced by

$$IS1_T(x_0, t_0) = \frac{\|\mathbf{x}(t)\|_T}{\|\mathbf{x}_0\|}.$$

There is a relation between the stability definition and the inverse stability measure. In the unstable case it is possible to choose the following series of initial conditions: $\|\mathbf{x}_{0_1}\| > \|\mathbf{x}_{0_2}\| > \dots > \|\mathbf{x}_{0_n}\|$. One has $\lim_{n \to \infty} \|\mathbf{x}_{0_n}\| = 0$ and, as an result, the $IS1_{\infty}(x_{0_n}, t_0) \to \infty$. But if the inverse stability measure is limited for all $x_0$, then $\|\mathbf{x}(t)\|_{\infty} \le s \|\mathbf{x}_0\|$. This means that the system is stable by Lyapunov.

The approach was used by Schiehlen to analyze a stability of the single and double pendulums. It is possible to select both stability and instability regions in place of the systems initial values. But it should be noticed that the choice of the calculation time $T$ is not discussed in paper [7].

## 3.2 The stability test (4) and Lyapunov exponents

“Lyapunov characteristic exponents” are introduced in some problems of stability or analysis of the “strange attractors” [2,15,17,22,23]. Positiveness of the exponents means the global mutual instability of phase trajectories; this is accepted as a criterion of the chaotic behavior onset in the nonlinear system under consideration. The largest Lyapunov exponent can be defined, for instance, as the following:

$$\lambda = \lim_{t \to \infty} \frac{1}{t} \ln \frac{|y(t)|}{|y(0)|} \qquad (7)$$
One has from here

$$ \frac{|y(t)|}{|y(0)|} = e^{(\lambda+\Delta)t}, \quad \text{where } \Delta \to 0 \text{ for } t \to \infty \qquad (8) $$

Comparing (8) and (5), one obtains the next inequalities

$$ \rho \le e^{(\lambda+\Delta)T}, \quad T \ge \frac{\ln \rho}{\lambda+\Delta} (\lambda > 0), \qquad (9) $$

where $T$ is the maximal calculation time. Thus the constants $\rho$ and $T$ presented in the criterion (5) are linked by the inequality (9). It is clear that a decrease of the value $\rho$ gives the corresponding decrease of the calculation time $T$ (it is confirmed by concrete computations too). Near the boundaries of the stability/instability regions values of the Lyapunov exponent $\lambda$ decrease, and the calculation time increase. But the value $\rho$ can not be chosen very small because near the boundaries values of $\lambda$ are small, and the calculation time is depended on the uncertain value $\Delta$. The significant arbitrariness in a choice of the value $\rho$ will be showed in the next section, too. We note that in the work by Rugonyi and Bathe [23] the concept of Lyapunov stability is connected with the Lyapunov characteristic exponents calculation. The effect of velocities (see the next Remark here) in the exponents calculation is analyzed too. Besides, the procedure proposed in [23] is used in the stability investigation of some continuous nonlinear systems, namely, the buckled beam and fluid-structure interaction systems. Additional discussion on the proposed approach and the Lyapunov exponents calculation approach is contained in the section Discussion.

REMARK. We can use in place of the test (5) the following “energetic” test. Namely, instability is fixed if

$$ y_2^2(t) + \dot{y}_2^2(t) \le \rho (y_2^2(0) + \dot{y}_2^2(0)) \quad (0 \le t \le T). \qquad (10) $$

It is clear, and checking calculations confirm it, that regions of instability in the system parameter space, obtained if we use (10) are greater that those obtained if we use test (5). This follows from the fact that there are regions where the variations $y_2(t)$ are small but the corresponding velocities $\dot{y}_2(t)$ are large.

### 3.3 Example

We show now the significant arbitrariness in a choice of the value $\rho$. For definiteness we consider the following system which can be obtained from (1) if we use some values of the physical and geometrical characteristics of the corresponding elastic system and use some scale transformations:

$$ \ddot{y}_1 + \delta \dot{y}_1 + [1 - q + y_1^2 + 4\epsilon y_2^2] y_1 = f \cos \omega t, \\ \ddot{y}_2 + \delta \dot{y}_2 + 4[4 - q + \epsilon y_1^2 + 4y_2^2] y_2 = 0. \qquad (11) $$

Consider the stability of the mode $y_2 = 0$, which can be both regular and chaotic. One selects the following domain of plane of parameters $(f, \epsilon)$: $0 \le f \le 10, 0 \le \epsilon \le 20$. Here $\delta = 0.1$. Mesh width in the domain was chosen as: $\Delta f = 1/3, \Delta \epsilon = 1$.

In the chosen scale of mesh, boundaries of stability and instability regions for the chosen region of the parameters plane stabilized in the interval $0 \le t \le T$ ($T \approx 250$).
Increasing a value of the constant $\rho$ presented in the test (5), we consider results of the calculations for $\rho = 10$; 5; 2.5. It is important that the stability regions are the same for the different values of $\rho$ and $q$. Corresponding results are presented in Fig. 1, where $q = 0.9$. A smoothing function is used in tracing of the graphic. The results show that for the system under consideration the localization of the vibration mode is the same as the stability of the mode.

<figure><img src="image_2.png" /><figcaption>Figure 1: Stability (S) and instability (U) regions of the mode $y = 0$ of equations (11), obtained using the stability test (6). Here $q=0.9$. Boundaries of stability and instability regions for the chosen region of parameters plane are stabilized on the interval of calculations, $0 \le t \le T$ ($T \approx 250$). Results are the same for $\rho = 10$; 5; 2.5.</figcaption></figure>

# 4 Stability of Forced Vibration Modes in Nonlinear Shells

Consider equations of the nonlinear shallow-shell theory connecting the normal deflection $w(x, y, t)$ and the stress function $\Phi(x, y, t)$ [24-27]:

$$ \frac{D}{h} \nabla^4 w = L(w, \Phi) + \frac{1}{R} \frac{\partial^2 \Phi}{\partial x^2} - \frac{\rho}{g} \frac{\partial^2 w}{\partial t^2} + q, \\ \frac{1}{E} \nabla^4 \Phi = -(1/2) L(w, w) - \frac{1}{R} \frac{\partial^2 w}{\partial x^2}, \quad (12) $$

where $x, y$ are respectively the longitudinal and circumferential coordinates; $E$ is the Young's modulus; the cylindrical stiffness $D = \frac{Eh^3}{12(1-\nu^2)}$; $\nu$ is the Poisson ratio; $h$ is the shell thickness; $R$ is the radius of the shell; $\rho$ is the density of the shell; $q(x, y, t)$ is the transversal loading;
the operator L is the following:

$$L(u, v) = \frac{\partial^2 u}{\partial x^2} \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 u}{\partial y^2} \frac{\partial^2 v}{\partial x^2} - 2 \frac{\partial^2 u}{\partial x \partial y} \frac{\partial^2 v}{\partial x \partial y}.$$

Circular cylindrical shell complete around the circumference is considered here. Let the shell be compressed along an element of the cylinder by uniformly distributed static axial forces after the shell has lost stability. One presents here the transverse loading as

$$q = F \cos \omega t \sin \frac{m\pi x}{l} \sin \frac{ny}{R}.$$

It means that the external periodical force excites the primary vibration mode with respect to the space coordinates $(x, y)$. Vibrations in the mode may be chaotic if we consider the post-buckling case. Two variants are possible. If the mode is orbitally stable then we obtain a localization of energy. But If the mode is unstable, the transfer of energy to other space mode is possible.

## 4.1 Companion modes in shells: analysis of stability

Let here the deflection $w$ be approximated by the following truncated Fourier series [26-28]:

$$w = f_1 \cos sy \sin rx + f_2 \sin sy \sin rx + \frac{n^2}{4R}(f_1^2 + f_2^2) \sin^2 rx, \quad (13)$$

where $s = \frac{n}{R}$, $r = \frac{m\pi}{l}$; $n, m$ are integer positive numbers; $f_1, f_2$ are unknown functions of time; the transverse loading is chosen of the form $q = f \cos \omega\tau \cos sy \sin rx$.

To analyze here a principal problem of the vibration modes energy pumping, we suppose that an initial deflection $w_0$ of the cylindrical shell is equal to zero. Obtaining the stress function $\Phi$ in the form of truncated Fourier series with respect to the space coordinates and using the standard Bubnov-Galerkin procedure, we obtain the equations to determine $f_1(t)$ and $f_2(t)$ which were exposed, for example, in [26,28].

One analyzes now a stability of the vibration mode $f_1 = f_1(t)$, $f_2 = 0$. Let us the variable $f_2$ emerges as a small variation. After a linearization by $f_2$ we can obtain the following equations:

$$\begin{gathered} \frac{d^2 f_1}{dt^2} \cdot (1 + 2\chi f_1^2) + \omega_0^2 f_1 + 2\chi f_1 \left(\frac{df_1}{dt}\right)^2 + \gamma_1 f_1^3 + g f_1^5 = f \cos \omega t, \\ \frac{d^2 f_2}{dt^2} + f_2 \left(\omega_0^2 + 2\chi \left(\left(\frac{df_1}{dt}\right)^2 + f_1 \frac{d^2 f_1}{dt^2}\right) + \gamma_1 f_1^2 + g f_1^4\right) = 0, \end{gathered} \quad (14)$$

where

$$\begin{gathered} \omega_0^2 = \frac{1}{\rho} \left[ \frac{D}{h} (s^2 + r^2)^2 + \frac{Er^4}{R^2 (s^2 + r^2)^2} \right], \quad 2\chi = \frac{2}{3} \left( \frac{n^2}{2R} \right)^2, \\ \gamma_1 = \frac{1}{\rho} \left[ \frac{E}{16} r^4 + \frac{Dn^4 r^4}{hR^2} - \frac{Er^4 s^4}{(s^2 + r^2)^2} \right], \quad g = \frac{3E}{16\rho} n^2 r^4 s^6 \left[ \frac{1}{(s^2 + r^2)^2} + \frac{1}{(s^2 + 9r^2)^2} \right]. \end{gathered}$$
One introduces the following change of variables in equations (14):

$$x = \frac{f_1}{h}, \quad y = \frac{f_2}{h}, \quad \tau = \omega_0 t, \quad \Omega = \frac{\omega}{\omega_0}, \quad F = \frac{f}{\omega^2 h}, \quad K = 2\chi h^2, \quad G = \frac{\gamma_1 h^2}{\omega_0^2}. \qquad (15)$$

We neglect later the small order terms $gf_1^5$ and $gf_1^4 f_2$. Besides, we introduce a coefficient $\alpha$ which is defined by the constant static axial force. In a case of $\alpha=0$ the force is absent; in a case of $\alpha=1$ (this is a critical loading) a trivial equilibrium position becomes unstable; in a case of $\alpha > 1$ (this is a post-buckling case) new equilibrium positions appear. Under external periodic excitation a transfer from periodic motions to irregular, apparently chaotic motions is possible. Introducing a dissipation, we can obtain, as a result, the following equations instead of the system (14):

$$x'' = (F \cos \Omega \tau - (1 - \alpha) x - Kxx'^2 - Gx^3 - \delta x') / (1 + Kx^2), \\ y'' = -y [(1 - \alpha) + K(x'^2 + xx'') + Gx^2 + \delta y'] \qquad (16)$$

The obtained equations (16) permit to investigate a stability both subcritical, and supercritical forced vibrations of cylindrical shells under the constant static axial force.

<figure><img src="image_3.png" /><figcaption>Figure 2: Instability regions of the vibration mode $x = x(t)$, $y = 0$ of the equations (16) in the plane $F, \Omega$: $0 \le F \le 1$, $0 \le \Omega \le 2$, obtained by using the test (6) ($x(0) = 0.5, \alpha = 0.5, \delta = 0.1, T = 400$).</figcaption></figure>

Next numerical results were obtained for the following values of geometrical and physical parameters of shells, and for the following circumferential and axial wavenumbers: $m = 1$, $n = 6$, $R = 1$, $h = 0.001$, $E = 2 \cdot 10^{11}$, $\rho = 7.8 \cdot 10^3$, $\nu = 0.3$, $l= 10$. A coefficient of dissipation was chosen as $\delta = 0, 0.05, 1$. A coefficient of the static loading was chosen as $\alpha = 0, 0.5, 0.9, 1.1$. It was considered the following values of the initial amplitude:
<figure><img src="image_4.png" /><figcaption>Figure 3: Instability regions of the vibration mode $x = x(t)$, $y = 0$ of the equations (16) in the plane $F, \Omega$: $0 \le F \le 1$, $0 \le \Omega \le 2$, obtained by using the test (6) ($x(0) = 0.5$, $\alpha = 0$, $\delta = 0.05$, $T = 1000$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Instability regions of the vibration mode $x = x(t)$, $y = 0$ of the equations (16) in the plane $F, \Omega$: $0 \le F \le 1$, $0 \le \Omega \le 2$, obtained by using the test (6) ($x(0) = 2$, $\alpha = 0$, $\delta = 0.05$, $T = 1000$).</figcaption></figure>

$x(0) = 0.5, 1, 2$; initial perturbations $y(0)$ are the following: $y(0) = 0.01 \cdot x(0)$. It is accepted that $x'(0) = y'(0) = 0$.

Calculations which were made by using the test (5) permit to determine stability/
<figure><img src="image_6.png" /><figcaption>Figure 5: Instability regions of the vibration mode $x = x(t)$, $y = 0$ of the equations (16) in the plane $F, \Omega$: $0 \le F \le 1$, $0 \le \Omega \le 2$, obtained by using the test (6) ($x(0) = 2$, $\alpha = 0.5$, $\delta = 0.05$, $T = 300$).</figcaption></figure>

instability regions of the vibration mode $x = x(t)$, $y = 0$ in the plane $F, \Omega$. One considers the following domain of the plane $F, \Omega$: $0 \le F \le 1$, $0 \le \Omega \le 2$. Mesh width in the domain was chosen as the following: $\Delta F = 0.05$, $\Delta \Omega = 0.09$.

In Figs. 2-5 a small part of the numerical calculation results is presented. Here points of the chosen mesh which correspond to the unstable regions, are marked by symbol. Time of the stability/ instability boundaries stabilization $T$ is pointed out too. In the instability regions only the coupled ($f_2 \neq 0$) companion modes (13), which were investigated by Kubenko et al. [28], are realised.

An increase of dissipation implies a corresponding decrease of the instability regions. Besides, an increase of the axial static stress coefficient $\alpha$ implies a corresponding increase of the instability regions. Results are the same for $x(0) = 0.5$ and $x(0) = 1$.

In a case of $\alpha > 1$ a trivial equilibrium position of the shell is unstable, and two new stable equilibrium positions appear, this is a post-buckling phase of the shell behavior. Under external periodical excitation chaotic vibrations can be appeared. The stability analysis by using the test (4) is not depends on the type (regular or chaotic) of the vibration mode under consideration. It is interesting that for $\alpha > 1$ both regular or chaotic vibration mode $f_1 = f_1(t)$, $f_2 = 0$ is unstable, that is only the coupled companion modes (13) are realised in this case.

Note that a similar analysis can be realized in a problem of stability of the forced vibration mode of the sloping arch which dynamics is described within a framework of the Kirchhoff hypothesis. A boundary of the stability/ instability regions is obtained, and the corresponding results were presented in paper [29].
# 5 Mutual Instability of Phase Trajectories

Consider the well-known nonautonomous Duffing equation,

$$ \dot{y} + \delta \dot{y} - \alpha y + \beta y^3 = f \cos \omega t \quad (17) $$

where the system parameters are: $\delta = 1$, $\alpha = 10$, $\beta = 100$, $\omega = 3.76$, and the value of the amplitude $f$ is variable. Examine the following region of phase plane of the equation (17): $0 \le y \le 0.66$; $0 \le \dot{y} \le 0.8$. Introduce some mesh in the defined region using the increments: $\Delta y = 0.02$, $\Delta \dot{y} = 0.016$.

<figure><img src="image_7.png" /><figcaption>Figure 6: Mutual stability/instability of phase trajectories of the nonautonomous Duffing equation (17). Corresponding initial values of the trajectories are presented on the following region of phase place: $0 \le y \le 0.66$; $0 \le \dot{y} \le 0.8$. The test (18) are used to determine the mutual stability /instability of trajectories outgoing from the near initial points (the time of stabilization here, $T < 100$). Results of the stability analysis are obtained for $f = 0.8$ (figure a), $f = 0.856$ (figure b) and $f = 0.9$ (figure c). The initial points of the chosen mesh which correspond to unstable trajectories are marked by dark squares.</figcaption></figure>

Points of the mesh $P_{ij}(y_{i0}, \dot{y}_{j0})$ will be chosen as initial points for the selected phase trajectories $y_{ij}^{(1)}(t)$ of the equation (17). Let us take other initial points near the chosen initial points $P_{ij}$, namely, $Q_{ij}(y_{i0} + \Delta y_{i0}, \dot{y}_{j0})$ where the value $\Delta y_{i0}$ is sufficiently small, and consider the other phase trajectory, $y_{ij}^{(2)}(t)$. One compares trajectories outgoing from the near initial points, using the limited Lyapunov test (5).

Instability of the outgoing trajectory is established if

$$ |y_{ij}^{(1)}(t) - y_{ij}^{(2)}(t)| \ge \rho |\Delta y_{i0}| \quad (0 \le t \le T). \quad (18) $$
We now select $\rho = 10$. Let $\Delta y_0 = 0.002$. Results of the stability analysis (the time of stabilization here, $T < 100$) are presented in Fig. 6 (a, b, c) for different values of the external amplitude: $f = 0.8$; 0.856; 0.9. In Fig. 6 the initial points of the chosen mesh, which correspond to unstable trajectories, are marked by dark squares. The calculations (in the chosen mesh of the equation phase place) show that for $f < 0.85$ the mutual instability of phase trajectories can be observed near the separatrix branches. Instability regions begin to extend if values of the external amplitude $f$ are slightly larger than 0.85, and the enlargement is very fast. The results generically correspond to known behavior of the nonautonomous Duffing system [1,2].

<figure><img src="image_8.png" /><figcaption>Figure 7: Mutual stability/instability of phase trajectories of the pendulum equation. There $-4 \le y \le 4$; $-2.4 \le \dot{y} \le 2.4$. The test (18) are used to determine the mutual stability/instability of trajectories outgoing from the near initial points (the time of stabilization here, $T = 50$). Results of the stability analysis are obtained for a) $f = 0.05$; b) $f = 0.2$; c) $f = 0.4$; d) $f = 0.6$; e) $f = 0.8$; f) $f = 1$. The initial points of the chosen mesh which correspond to unstable trajectories are marked by dark dots.</figcaption></figure>

In Fig. 7 the results of stability analysis for the equation of the pendulum with excited point of suspension are presented. The equation of motion is the following:

$$ \ddot{y} + \delta\dot{y} + (1 + f \cos \omega t) \sin y = 0. $$

The calculations were made with the friction $\delta = 0.1$, the force frequency $w = 2$, the time of calculations $T = 50$, the widths of the mesh $\Delta y = 0.1$, $\Delta \dot{y} = 0.1$ and $\Delta y_0 = 0.01$. $\rho$ equals 10. Obtained results correspond to the behavior of this system described in book of Neymark, Landa [30].
# 6 Discussion

It is appropriate to compare the presented approach with the well known methods of the Lyapunov characteristic exponents calculations. We can note that it is not necessary to analyze here variational equations and calculate corresponding eigenvalues because only some direct consequence from the classical Lyapunov stability definition is used in the proposed approach. The Lyapunov exponents are introduced, for example, in [22, p.67 and further] as a limit of the eigenvalues when the independent variable tends to infinity. Of course, we can obtain a smaller information about the system behavior than from the Lyapunov exponents calculation approach, but we avoid some known computing difficulties [22, p.73 and further]. Note that contemporary procedures to calculate the Lyapunov exponents both in discrete and continuous nonlinear systems are described in [23].

Next, we think that in numerous works on numerical approaches a problem of choice of the calculation time is not explained until now, although some discussion about a selection of the Lyapunov exponents calculation time is presented, for example, in [22, p.81]. Unfortunately, the principle "to calculate as long as possible" is dominated in many works. It seems to us that we proposed a simple test for determination of the calculation time starting from the real calculations at points on some chosen mesh of the system under consideration parameter space.

# 7 Conclusions

We analyzed in this paper the stability of regular or chaotic vibration modes in two-DOF multistable nonlinear systems. The stability problem for chaotic modes has no analytical solution because the chaotic motions cannot be described without a computer simulation. A new stability test is proposed here. The test permits to implement on computer the classical Lyapunov stability definition for a case when initial variations are not arbitrary small, that is they are limited below. Note that a realization of this approach is not difficult because only the standard Runge-Kutta procedure with an additional stability condition is used. Calculations by computer which correspond to the test are conducted for points of some chosen mesh of the system parameter space and continued until boundaries of stability and instability regions (in the mesh scale) are stabilized. This is the principal criterion for a choice of the calculation time. Stability of regular or chaotic vibration modes of some nonlinear elastic system is considered.

The methodology presented in this work is sufficiently general to be applicable to stability problems for vibration modes in different types of nonlinear systems.

The approach of determination of the chaos onset, which is based on the Lyapunov stability definition, is also proposed. Mutual instability of phase trajectories is used as a criterion of the chaotic behavior of dynamical system. An effectiveness of this approach is confirmed by calculations of the chaos onset of the nonautonomous Duffing equation and pendulum system. The methodology can be considered as an alternative to the well-known Poincare' map approach, and can be applicable to a determining of the chaos onset in other dynamical systems.
# 8 Acknowledgments

The authors are very thankful to Dr. Konstantin Avramov for a very useful discussion on the paper.

# 9 References

[1] P.J.Holmes, A nonlinear oscillator with a strange attractor, *Philos. Trans. Royal Soc.*, London **A292**, 419-448 (1979).

[2] F.C.Moon, *Chaotic Vibrations*, Wiley, New York 1987.

[3] A.M.Liapunov, *Stability of Motion*, Academic Press, New York 1966.

[4] E.L.Ince, *Ordinary Differential Equations*, Longmans Green, London 1926.

[5] A.A.Andronov, E.A.Vitt and S.E.Khaikin, *Theory of Oscillations*, Pergamon Press, Oxford 1966.

[6] N.Minorsky, *Nonlinear Oscillations*, Van Nostrand, Princeton 1962.

[7] W.Schiehlen, Nonlinear oscillations in multibody systems. In *Proc. of 1st European Nonlinear Oscillations Conference*. Hamburg, 1993, Akademie Verlag Gmbh, Berlin, (Mathematical Research, vol.72), 85-106 (1993).

[8] Yu.V.Mikhlin, Stability of periodical and chaotic vibrations in systems with more than one equilibrium positions, In *Proc.of MATHMOD*, Vienna, pp. 903-905 (1994).

[9] Yu.V.Mikhlin, A.V.Nemchin, Stability of periodic and chaotic modes in systems with several equilibrium positions, In *Proc. of the 4th Conference on Dynamical Systems – Theory and Applications*, Lodz (Poland), pp.235-240 (1997).

[10] H.Poincare, Les Methods Nouvelles de la Mecanique Celeste, 3 Vols. Gauthier-Villars: Paris 1899.

[11] G.D.Birhoff, Sur l'existence de regions d'instabilite en Dynamique, Ann. Inst. H. Poincare, 2, 369-386 (1932).

[12] E.N. Lorenz, Deterministic non-periodic flow, J. Atmos. Sci., 20, 130-141 (1963).

[13] D.Ruelle and F.Takens, On the nature of turbulence, Common. Math. Phys., 167-192 (1971).

[14] V.K.Melnikov, On the stability of the centre for time periodic perturbations, Trans. Moscow Math. Soc., 12, 1-57 (1963).

[15] J.Guckenheimer, P.Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag, New York 1993.

[16] Y.Ueda, Randomly transitional phenomena in the system governed by Duffing's equation, J.Stat.Phys 20,181-196 (1979).

[17] A.Wolf, Quantifying chaos with Lyapunov exponents, Nonlinear Sci. Theory Appl. Ed. A.V.Holden, Manchester Univ. Press (1984).

[18] H.Kauderer, *Nichtlinear Mechanik*. Springer-Verlag, Berlin 1958.

[19] L.Perko, *Differential Equations and Dynamical Systems*, Springer-Verlag, New York 1991.

[20] C.L.Siegel and J.K.Moser, *Lectures on Celestial Mechanics*, Springer-Verlag, New York. Heidelberg, Berlin 1971.

[21] G.Pecelli and E.S.Thomas, Normal modes, incoupling, and stability for a class of non-linear oscillators, *Quart.of Appl.Math.* **37**, 281-301 (1979).
[22] T.S.Parker and L.O.Chua, *Practical numerical algorithms for chaotic systems*, Springer-Verlag, New York 1989.

[23] S.Rugonyi and K.-J.Bathe, An evaluation of the Lyapunov characteristic exponent of chaotic continuous system, *Int.J.for Numerical Methods in Engineering. 56*, 145-163 (2003).

[24] Th.Karman and H.S.Tsien, The buckling of thin cylindrical shells under axial compression, *J.Aeron. Sci.* **8**, 303-312 (1941).

[25] E.H.Dowell and C.S.Ventres, Modal equations for the nonlinear flexural vibrations of a cylindrical shell, *Int. J. Solids and Structures* **4**, 975-991 (1968).

[26] *Thin-Shell Structures. Theory, Experiment and Design*, Ed. by Y.C.Fung and E.E.Sechler, Prentice-Hall, Englewood cliffs, New Jersey 1974.

[27] A.S.Volmir, *Nonlinear Dynamics of Plates and Shells*, Nauka, Moscow 1972 [In Russian].

[28] V.D.Kubenko, P.S.Kovalchuk, T.S.Krasnopolskaja, *Nonlinear interaction of the bending vibrations modes of cylindrical shells*, Naukova Dumka, Kiev 1984 [In Russian].

[29] Yu.V.Mikhlin, Stability of regular or chaotic post-buckling vibration modes of elastic shells, *Nonlinear Dynamics of Shells and Plates* (Ed. by M.P.Paidoussis, M.Amabili, P.B.Gonsalves), 2000 ASME International Mechanical Engineering Congress, Orlando, 95-103. 2000.

[30] Yu.I.Neymark, P.S.Landa, Stochastic and Chaotic vibrations, Nauka, Moscow, 1987 [In Russian].