Stability boundary approximation of periodic dynamics

TL;DR

This paper introduces an averaging-based method to approximate stability boundaries in parameter space for systems with parametric excitation, demonstrated on a 2-DOF pendulum with damping.

Contribution

It provides a novel approach for estimating stability regions using fourth order approximations of the monodromy matrix for systems with parametric excitation.

Findings

Small damping shifts stability boundaries upwards.

The method accurately predicts stability domains.

Application to a pendulum demonstrates effectiveness.

Abstract

We develop here the method for obtaining approximate stability boundaries in the space of parameters for systems with parametric excitation. The monodromy (Floquet) matrix of linearized system is found by averaging method. For system with 2 degrees of freedom (DOF) we derive general approximate stability conditions. We study domains of stability with the use of fourth order approximations of monodromy matrix on example of inverted position of a pendulum with vertically oscillating pivot. Addition of small damping shifts the stability boundaries upwards, thus resulting to both stabilization and destabilization effects.

Equations

$$\dot{\mathbf{X}}(t)=\mathbf{J}(t)\cdot \mathbf{X}(t),\mathbf{X}(0)=\mathbf{I},$$

$$\mathbf{J}(t)={\mathbf{J}}_0(t)+{\mathbf{J}}_1(t)+{\mathbf{J}}_2(t)+{\mathbf{J}}_3(t)+\dots ,$$

$$\dot{\mathbf{Y}}(t)=\mathbf{H}(t)\cdot \mathbf{Y}(t),\mathbf{Y}(0)=\mathbf{I},$$

$$\mathbf{H}(t)={\mathbf{H}}_1(t)+{\mathbf{H}}_2(t)+{\mathbf{H}}_3(t)+\dots ,$$

$$\mathbf{Y}(t)=\left(\mathbf{I}+{\mathbf{U}}_1(t)+{\mathbf{U}}_2(t)+\dots \right)\cdot \mathbf{Z}(t),$$

$$\dot{\mathbf{Z}}(t)=\mathbf{A}\cdot \mathbf{Z}(t),\mathbf{Z}(0)=\mathbf{I},$$

$$\mathbf{Z}(t)=\exp \left(\left[{\mathbf{A}}_1+{\mathbf{A}}_2+{\mathbf{A}}_3+\dots \right]t\right).$$

$$\dot{\mathbf{Y}}(t)=$$

$$\left({\mathbf{H}}_1(t)+{\mathbf{H}}_2(t)+\dots \right)\cdot \left(\mathbf{I}+{\mathbf{U}}_1(t)+{\mathbf{U}}_2(t)+\dots \right)\cdot \mathbf{Z}$$

$$=\left({\dot{\mathbf{U}}}_1(t)+{\dot{\mathbf{U}}}_2(t)+\dots \right)\cdot \mathbf{Z}+\left(\mathbf{I}+{\mathbf{U}}_1(t)+{\mathbf{U}}_2(t)+\dots \right)\cdot \left({\mathbf{A}}_1+{\mathbf{A}}_2+\dots \right)\cdot \mathbf{Z},$$

$${\mathbf{A}}_1=\frac{1}{T}{\int }_0^T{\mathbf{H}}_1(t)\mathrm{d}t,$$

$${\mathbf{U}}_1(t)={\int }_0^t($$

$${\mathbf{A}}_2=\frac{1}{T}{\int }_0^T($$

$${\mathbf{U}}_2(t)={\int }_0^t({\mathbf{H}}_2(\tau )-{\mathbf{A}}_2+$$

$${\mathbf{A}}_3=\frac{1}{T}{\int }_0^T({\mathbf{H}}_3(t)+$$

$$+$$

$${\mathbf{U}}_n(t)={\int }_0^t({\mathbf{H}}_n(\tau )-{\mathbf{A}}_n+$$

$$\dots$$

$$+$$

$${\mathbf{A}}_{n+1}=\frac{1}{T}{\int }_0^T({\mathbf{H}}_{n+1}(t)+$$

$$\dots$$

$$+$$

$$\mathbf{F}$$

$${\mathbf{F}}_1={\mathbf{F}}_0\cdot {\mathbf{Z}}_1(T)={\mathbf{F}}_0\cdot {\mathbf{A}}_{1}T.$$

$${\mathbf{F}}_2={\mathbf{F}}_0\cdot \left({\mathbf{A}}_{2}T+\frac{1}{2}{\mathbf{A}}_1^2{T}^2\right).$$

$$\displaystyle\mathbf{F}_{3}=\mathbf{F}_{0}\cdot\Big{(}\mathbf{A}_{3}T+\frac{1}{2}\left(\mathbf{A}_{1}\cdot\mathbf{A}_{2}+\mathbf{A}_{2}\cdot\mathbf{A}_{1}\right)T^{2}+\frac{1}{6}\mathbf{A}_{1}^{3}T^{3}$$

$$\displaystyle\mathbf{F}_{4}=\mathbf{F}_{0}\cdot\Big{(}\mathbf{A}_{4}T$$

$$\displaystyle+\frac{1}{6}(\mathbf{A}_{1}^{2}\mathbf{A}_{2}+\mathbf{A}_{1}\mathbf{A}_{2}\mathbf{A}_{1}+\mathbf{A}_{2}\mathbf{A}_{1}^{2})T^{3}+\frac{1}{24}\mathbf{A}_{1}^{4}T^{4}\Big{)},$$

$${\rho }^2-\mathrm{tr}(\mathbf{F})\rho +\det (\mathbf{F})=0.$$

$$|\mathrm{tr}(\mathbf{F})|-1\le \det (\mathbf{F})\text{and}\det (\mathbf{F})\le 1,$$

$$\mathrm{tr}\left(\mathbf{F}\right)=\mathrm{tr}({\mathbf{F}}_0)+\mathrm{tr}({\mathbf{F}}_1)+\mathrm{tr}({\mathbf{F}}_2)+\mathrm{tr}({\mathbf{F}}_3)+\mathrm{tr}({\mathbf{F}}_4)+\dots ,$$

$$\det (\mathbf{F})$$

$$\mathrm{tr}({\mathbf{A}}_j)=\frac{1}{T}{\int }_0^T\mathrm{tr}({\mathbf{J}}_j(\tau ))d\tau .$$

$$\displaystyle\det\!\left(\mathbf{F}\right)=e^{\int^{T}_{0}\operatorname{tr}(\mathbf{J}_{0}(\tau)){\rm\,d}\tau}\Big{(}1$$

Full text

1 Introduction

Let us study the stability of the equilibrium $y\equiv 0$ of a nonlinear system governed by ordinary differential equations $\dot{y}=\mathbf{J}(t)\,y+f(t,y)$, where for nonlinear vector-function $f$ exist constants $c$ and $\nu>1$, such that $|f(t,y)|\leq c\,|y|^{\nu}$ for all $t\geq 0$ and $y\in\mathbb{R}^{n}$. Vector-function $f(t,y)$ and matrix $\mathbf{J}(t)$ are piecewise continuous in $t$. Moreover $\mathbf{J}(t)$ is bounded and $T$-periodic.

According to Lyapunov’s theorem the trivial solution of such a nonlinear system is asymptotically stable if all Lyapunov exponents of the linear system $\dot{x}(t)=\mathbf{J}(t)\,x(t)$ are strictly negative and the solution is unstable if at least one Lyapunov exponent is strictly positive.111Lyapunov regularity condition of the linear system holds when its matrix is periodic. Asymptotic stability of the periodic linear system determines the asymptotic stability of the nonlinear system equilibrium and vice versa.222Periodic linear system is asymptotically stable if, and only if, it is exponentially stable. Exponential stability of $\dot{x}(t)=\mathbf{J}(t)\,x(t)$ results in exponential stability, and hence asymptotic stability, of the nonlinear system solution $y\equiv 0$. Stability of linear systems with time-periodic coefficients was also studied by Gaston Floquet [1]. Lyapunov exponent $\lambda$ of the linear periodic system can be expressed via its corresponding Floquet multiplier $\rho$ as $\lambda=\frac{1}{T}\ln|\rho|$. The theorem can also be reformulated to compare absolute values of Floquet multipliers with 1. Floquet multipliers are the eigenvalues of the monodromy matrix which is the fundamental matrix of the linear system taken at time $T$.

Monodromy matrix and Floquet multipliers can always be numerically calculated. Stability can also be checked by studying solutions of the system, see, e.g. [4, 5, 6] and references therein. But in practice it is often useful to have analytical approximations of stability regions in parameter space. To obtain a straightforward technique for deriving such analytical stability boundary approximations of any order we combine Floquet theory with asymptotic method of averaging, [6].

This technique yields same results as expansion of monodromy matrix in series in [2] up to the terms of higher order than the order of approximation. These terms are automatically eliminated in the averaging scheme making the technique more convenient in practice.

We demonstrate the proposed technique on the example of inverted pendulum, where destabilizing effect of damping (shift of the lower stability boundary) is manifested in the fourth approximation, see [7], though it would be natural to expect stabilization by damping, see, e.g. [4, 8]. But there is also destabilizing effect (shift of the upper boundary) demonstrated numerically in [9]. Here we obtain approximations of both stability boundaries and study analytically both effects of stabilization and destabilization by damping, Moreover, the fourth approximation yields in addition boundaries of another stability domain.

2 Statement of the problem

Consider linearization, $\dot{x}(t)=\mathbf{J}(t)\,x(t)$, of a nonlinear system about its equilibrium, where $x(t)$ is the vector of state variable perturbations and $\mathbf{J}(t)$ is piecewise continuous, $T$-periodic and thus integrable Jacobian matrix of the original nonlinear system. Solution of the matrix differential equation with the initial value being the identity matrix $\mathbf{I}$

[TABLE]

yields fundamental matrix and the monodromy matrix $\mathbf{F}=\mathbf{X}(T)$. If all eigenvalues of the monodromy matrix, Floquet multipliers, have absolute values smaller than one, then the equilibrium of the nonlinear system is asymptotically stable, and if at least one eigenvalue has absolute value grater than one, then the equilibrium is unstable, see, e.g. [2]. To have analytical approximations of stability regions in parameter space we apply the following.

3 Averaging scheme

Let the Jacobian matrix $\mathbf{J}(t)$ be expended into the series

[TABLE]

where the lower index denotes the order of smallness. Suppose we know solution $\mathbf{X}_{0}(t)$ of the matrix initial value problem $\dot{\mathbf{X}}_{0}(t)=\mathbf{J}_{0}(t)\cdot\mathbf{X}_{0}(t)$, where $\mathbf{X}_{0}(0)=\mathbf{I}$. Then the change of variable $\mathbf{X}(t)=\mathbf{X}_{0}(t)\cdot\mathbf{Y}(t)$ converts (1) to the standard form:

[TABLE]

where matrix $\mathbf{H}(t):=\mathbf{X}_{0}^{-1}(t)\cdot\left(\mathbf{J}(t)-\mathbf{J}_{0}(t)\right)\cdot\mathbf{X}_{0}(t)$ is small for $t\in[0,T]$. Approximate solution of (3) can be found with averaging method as follows. Let

[TABLE]

where $\mathbf{H}_{j}(t):=\mathbf{X}_{0}^{-1}(t)\cdot\mathbf{J}_{j}(t)\cdot\mathbf{X}_{0}(t)$ for all $j\geq 1$. We will find solution as

[TABLE]

where $\mathbf{U}_{j}(t)$ are matrix-functions, such that $\mathbf{U}_{j}(0)=\mathbf{U}_{j}(T)=0$ and $\mathbf{Z}(t)$ is the solution of the autonomous averaged differential equation:

[TABLE]

where $\mathbf{A}=\mathbf{A}_{1}+\mathbf{A}_{2}+\mathbf{A}_{3}+\ldots$, which can be written via the matrix exponential:

[TABLE]

The matrices $\mathbf{A}_{j}$ and matrix-functions $\mathbf{U}_{j}(t)$ can be found one by one the following expressions.333The expressions are derived by differentiating (5) w.r.t. time

$\displaystyle\dot{\mathbf{Y}}(t)=$ $\displaystyle\left(\dot{\mathbf{U}}_{1}(t)+\dot{\mathbf{U}}_{2}(t)+\ldots\right)\cdot\mathbf{Z}(t)+\left(\mathbf{I}+\mathbf{U}_{1}(t)+\mathbf{U}_{2}(t)+\ldots\right)\cdot\dot{\mathbf{Z}}(t),$

substituting there expressions for time derivatives from (3) and (6)

$\displaystyle\left(\mathbf{H}_{1}(t)+\mathbf{H}_{2}(t)+\ldots\right)\cdot\left(\mathbf{I}+\mathbf{U}_{1}(t)+\mathbf{U}_{2}(t)+\ldots\right)\cdot\mathbf{Z}$

$\displaystyle=\left(\dot{\mathbf{U}}_{1}(t)+\dot{\mathbf{U}}_{2}(t)+\ldots\right)\cdot\mathbf{Z}+\left(\mathbf{I}+\mathbf{U}_{1}(t)+\mathbf{U}_{2}(t)+\ldots\right)\cdot\left(\mathbf{A}_{1}+\mathbf{A}_{2}+\ldots\right)\cdot\mathbf{Z},$

collecting there terms of the same order, and canceling non-degenerate matrix $\mathbf{Z}$, which yield the following equalities. First order: $\mathbf{H}_{1}(t)=\dot{\mathbf{U}}_{1}(t)+\mathbf{A}_{1}.$Second order: $\mathbf{H}_{2}(t)+\mathbf{H}_{1}(t)\cdot\mathbf{U}_{1}(t)=\dot{\mathbf{U}}_{2}(t)+\mathbf{U}_{1}(t)\cdot\mathbf{A}_{1}+\mathbf{A}_{2},$Third order: $\mathbf{H}_{3}(t)+\mathbf{H}_{1}(t)\cdot\mathbf{U}_{2}(t)+\mathbf{H}_{2}(t)\cdot\mathbf{U}_{1}(t)=\dot{\mathbf{U}}_{3}(t)+\mathbf{U}_{2}(t)\cdot\mathbf{A}_{1}+\mathbf{U}_{1}(t)\cdot\mathbf{A}_{2}+\mathbf{A}_{3},$ Fourth order: $\mathbf{H}_{4}(t)+\mathbf{H}_{1}(t)\cdot\mathbf{U}_{3}(t)+\mathbf{H}_{2}(t)\cdot\mathbf{U}_{2}(t)+\mathbf{H}_{3}(t)\cdot\mathbf{U}_{1}(t)=\dot{\mathbf{U}}_{4}(t)+\mathbf{U}_{3}(t)\cdot\mathbf{A}_{1}+\mathbf{U}_{2}(t)\cdot\mathbf{A}_{2}+\mathbf{U}_{1}(t)\cdot\mathbf{A}_{3}+\mathbf{A}_{4}$, and so on…

For the first order approximation we calculate $\mathbf{A}_{1}$ as the average of $\mathbf{H}_{1}(t)$

[TABLE]

under assumption that $\mathbf{A}_{1}$ and $\mathbf{H}_{1}(t)$ are of the same order of smallness. In particular we assume that $\mathbf{H}_{1}(t)$ does not contain periodic functions with small frequencies, which would appear in the denominator during integration and could cause high value of $\mathbf{A}_{1}$, thus violating the assumption of its smallness.

For the second order approximation we have to calculate

[TABLE]

using matrix $\mathbf{A}_{1}$ already obtained in (8).

For the third order approximation we have

[TABLE]

and so on … for the $n+1$-th order approximation we calculate

[TABLE]

4 Monodromy matrix approximation

Due to $\mathbf{Y}(T)=\mathbf{Z}(T)$ we have

[TABLE]

where we denote $\mathbf{F}_{0}:=\mathbf{X}_{0}(T)$ as the zero order approximation of monodromy matrix, $\mathbf{F}\approx\mathbf{F}_{0}$. Hence, we can write expressions to find terms of the expansion $\mathbf{F}=\mathbf{F}_{0}+\mathbf{F}_{1}+\mathbf{F}_{2}+\mathbf{F}_{3}+\mathbf{F}_{4}+\ldots$, where $\mathbf{F}_{j}=\mathbf{F}_{0}\cdot\mathbf{Z}_{j}(T)$ with $\mathbf{Z}_{0}(T)=\mathbf{I}$. Expansion of the matrix exponential in (7) yields expressions for $\mathbf{F}_{j}$ via $\mathbf{A}_{k}$, where $k\leq j$.

For the first order approximation we have $\mathbf{Z}_{1}(T)=\mathbf{A}_{1}T$, so that

[TABLE]

For the second order approximation the expansion of the matrix exponential in (7) up to the second order terms yields $\mathbf{Z}_{2}(T)=\mathbf{A}_{2}T+\frac{1}{2}\mathbf{A}_{1}^{2}T^{2}$ and

[TABLE]

For the third order approximation we have

[TABLE]

For the fourth order approximation we calculate

[TABLE]

and so on …

5 Stability conditions in 2-dimensional case

The eigenvalues of the monodromy matrix, Floquet multipliers, determine the stability of the solution $x=(0,0)^{\prime}$ of the linearized system. Since $\mathbf{F}$ is the $2\!\times 2$ matrix its eigenvalues can be found analytically as roots $\rho_{1}$ and $\rho_{2}$ of the characteristic polynomial:

[TABLE]

Stability conditions ($|\rho_{1}|\leq 1$ and $|\rho_{2}|\leq 1$) written in the case of real roots as $\rho\in[-1,1]$ and in the case of complex conjugate roots as $\rho_{1}\rho_{2}\leq 1$, with the use of (20) and Vieta’s formula $\rho_{1}\rho_{2}=\det(\mathbf{F})$ correspondingly, take the form

[TABLE]

where for asymptotic stability all inequalities should be strict, see, e.g. [12], p. 213. For instability, it is sufficient that at least one of the conditions in (21) is violated.

Let us find approximations of stability conditions. Trace can be written as

[TABLE]

with $\mathbf{F}_{1}$, $\mathbf{F}_{2}$, $\mathbf{F}_{3}$, and $\mathbf{F}_{4}$ calculated by (16), (17), (18), and (19).

Notice that due to Liouville’s formula, see, e.g. [3], we have

[TABLE]

and since $\mathbf{F}=\mathbf{F}_{0}\cdot\mathbf{Z}(T)$ and (7) determinant can be written as

[TABLE]

so that for any dimension for all $j\geq 1$ we have

[TABLE]

Then expansion of the matrix exponential yields

[TABLE]

The first order approximation of stability conditions $|\operatorname{tr}(\mathbf{F})|-1\leq\det(\mathbf{F})\leq 1$ can be written with the use of (22) and (24) as

[TABLE]

and so on …

6 Two DOF system with impulse parametric excitation

6.1 High frequency stabilization of inverted pendulum

Let’s show that motion of an inverted pendulum can be stable if we supply to the suspension point rather high frequency vibrations in vertical direction $l\ddot{\varphi}+bl\dot{\varphi}+\left(g\pm c\right)\sin(\varphi)=0$, where $l$ is the length of the pendulum, $a\ll l$ is the amplitude of the vibrations of the pivot. The period of the pivot vibrations we normalize to be $T=2\pi$, besides in any semi-period acceleration of the pivot is constant and is equal $\pm c$, which sign changes each semi-period. We assume linear viscous friction with coefficient $b$. It turns out that for rather low relative eigenfrequency $\omega\ll 1$ the inverted position becomes stable. The equation of motion can be written in the form

[TABLE]

where $\omega^{2}=g/l$ is the relative eigenfrequency, $\varepsilon=c/l$ is the relative excitation acceleration, and $\beta=b/g$ is the new damping coefficient. Stability condition for this problem without damping can be found in [10]. The linearized case without damping coincides with the Meissner equation, [11].

This system, linearized about inverted vertical position $(\phi,s)=(\pi,0)$, has the form $\dot{x}(t)=\mathbf{J}(t)\,x(t)$, with vector $x(t)$ corresponding to the perturbation of vector $(\varphi(t),s(t))^{\prime}$ and $\mathbf{J}(t)$ being the piecewise constant Jacobian matrix of the original system: $\mathbf{J}(t)=\mathbf{J}^{+}$ if $t\in\left[0,\pi\right)$ and $\mathbf{J}(t)=\mathbf{J}^{-}$ if $t\in\left[\pi,2\pi\right)$, where

[TABLE]

It is easy to find exact stability conditions to which approximate conditions converge as we shall demonstrate.

6.2 Exact stability conditions

We have the following explicit expression of the monodromy matrix via matrix exponentials

[TABLE]

Since the determinant of a matrix product equals the product of the determinants, we have from (28)

[TABLE]

where we take into account that $\det(\exp\!\left(\pi\,\mathbf{J}^{\pm}\right))=e^{2\pi\operatorname{tr}(\mathbf{J}^{\pm})}.$ The same expression as (29), $\det(\mathbf{F})=e^{-2\pi\beta\omega}$, can be obtained by Liouville’s formula for any piecewise continuous integrable $2\pi$-periodic modulation function, see [3]. So with positive damping coefficient, $\beta>0$, asymptotic stability can only be lost when the first condition in (21) is violated, i.e. when

[TABLE]

We compare exact stability borders, determined by (27), (28), and (30) as

[TABLE]

with approximate stability boundaries obtained for this example.

6.3 Approximate stability conditions

We expend this matrix in the series $\mathbf{J}^{\pm}=\mathbf{J}_{0}+\mathbf{J}_{1}^{\pm}+\mathbf{J}_{2}$, where

[TABLE]

assuming that $\varepsilon$, $\omega$, and $\beta$ have the same order of smallness.

[TABLE]

According the the formula $\mathbf{H}_{j}(t):=\mathbf{X}_{0}^{-1}(t)\cdot\mathbf{J}_{j}(t)\cdot\mathbf{X}_{0}(t)$ we have

[TABLE]

Formula (8) reads as

[TABLE]

The zero order approximation of monodromy matrix is the following

[TABLE]

With (16) we find the first order adjustment of monodromy matrix

[TABLE]

Thus $\operatorname{tr}(\mathbf{F}_{0})=2$, $\operatorname{tr}(\mathbf{F}_{1})=0$, and $\operatorname{tr}(\mathbf{A}_{1})=0$.

For the second order approximation we take \mathbf{H}_{2}(t)=\left(\begin{array}[]{c@{\quad}c}-t\omega^{2}&-t\left(\omega^{2}t-\beta\omega\right)\\ \omega^{2}&\omega^{2}t-\beta\omega\\ \end{array}\right) and obtain with (9) and (10) the matrix

[TABLE]

We calculate according (17) the second order adjustment of monodromy matrix

[TABLE]

Thus $\operatorname{tr}(\mathbf{F}_{2})=-\frac{1}{3}\,{\pi}^{4}{\varepsilon}^{2}+4{\pi}^{2}{\omega}^{2}-2\pi\,\beta\,\omega$ and $\operatorname{tr}(\mathbf{A}_{2})=-2\pi\,\beta\,\omega$.

Second order approximation of stability border written from (30) as

[TABLE]

yields $4\,{\pi}^{2}{\omega}^{2}-\varepsilon_{p}^{2}{\pi}^{4}/3=0$ and $4-4\pi\,\beta\,w+4{\pi}^{2}{\omega}^{2}-\varepsilon_{n}^{2}{\pi}^{4}/3=0$ in cases of positive and negative value of the sum $\operatorname{tr}(\mathbf{F}_{0})+\operatorname{tr}(\mathbf{F}_{1})+\operatorname{tr}(\mathbf{F}_{2})$ correspondingly. Hence we have corresponding stability borders denoted by indexes $p$ and $n$

[TABLE]

same as in the third approximation, see dashed lines in the Figure on the left, because (11)–(18) yield $\operatorname{tr}(\mathbf{F}_{3})=0$ and $\operatorname{tr}(\mathbf{A}_{3})=0$.

Fourth order approximation, where stability boundary equation reads as

[TABLE]

yields in cases of positive and negative sums of traces the following two equations:

$\footnotesize\frac{{\pi}^{8}{\varepsilon}^{4}_{p}}{1260}-\frac{{\pi}^{4}}{3}\!\left(1+{\frac{4{\pi}^{2}{\omega}^{2}}{15}}-{\pi}\beta\omega\right){\varepsilon}^{2}_{p}+4{\pi}^{2}{\omega}^{2}\left(1+\frac{{\pi}^{2}{\omega}^{2}}{3}-\beta{\omega}{\pi}\right)=0,$

$\footnotesize\frac{{\pi}^{8}{\varepsilon}^{4}_{n}}{1260}-\frac{{\pi}^{4}}{3}\!\left(1+{\frac{4{\pi}^{2}{\omega}^{2}}{15}}-{\pi}\beta\omega\right){\varepsilon}^{2}_{n}+4\left(1-\beta\pi\omega+{\pi}^{2}{\omega}^{2}\!\left(1+\frac{{\pi}^{2}{\omega}^{2}}{3}-\beta{\omega}{\pi}+{\beta}^{2}\right)\right)\!=\!0.$

Solutions of these two equations with respect to $\varepsilon^{2}$ give us four borders, drawn in the Figure with solid lines, which approximate two exact stability domains determined by (31) and marked in gray.

7 Conclusion

We develop convenient algorithm for obtaining approximate stability boundaries of parametrically excited systems. We demonstrate how this algorithm can be applied to the particular case of parametric pendulum for obtaining approximate stability domains even in the case of damping and impulse parametric excitation. Stabilizing and destabilizing effects of damping on inverted equilibrium of the parametric pendulum are revealed with the use of the fourth order approximation of stability boundaries.

In the right figure we draw approximate stability boundaries of inverted vertical pendulum position. Addition of small linear viscous friction $\beta$ shifts both stability boundaries upward. Thus, at the lower boundary additional friction destabilizes the inverted pendulum while at the upper boundary friction stabilizes the pendulum position.

8 Acknowledgments

A.O. Belyakov received funding from the Russian Science Foundation grant 19-11-00223.