Исследование вынужденных колебаний линейной системы с двумя степенями свободы и малым параметром методом ляпунова–шмидта

Аннотация

Работа посвящена нахождению периодических решений линейных систем с двумя степенями свободы и малым параметром методом Ляпунова-Шмидта. На основе метода Ляпунова-Шмидта разработан алгоритм в пакете Maple, построены графики компонент периодических решений и фазовых траекторий возмущенной системы.

Полный текст

Введение.

Математические модели в линейном приближении могут быть описаны как неоднородные линейные системы обыкновенных дифференциальных уравнений с малым параметром. Общая методика исследования периодических решений таких систем описана в работе [1]. Она основана на представлении исследуемой системы в виде операторного уравнения в банаховом пространстве и применении метода Ляпунова-Шмидта [2–4].

В работе [5] приведено исследование периодических решений одной линейной неоднородной системы второго порядка с малым линейным возмущением. В работе [6] проведено исследование вынужденных колебаний одной линейной системы двух связанных осцилляторов с малым параметром.

В настоящей работе исследуются вынужденные колебания линейной системы с двумя степенями свободы и малым параметром при условии, что на систему действует внешняя периодическая сила с двумя соизмеримыми частотами.

Математическая модель вынужденных колебаний системы с двумя степенями свободы и малым параметром.

Рассмотрим математическую модель вынужденных колебания системы с двумя степенями свободы в виде [7]

$\left\{\begin{array}{l}{\ddot{Q}}_1+n_1^2{Q}_1-K{Q}_2=F_1\left(t\right),\\ {\ddot{Q}}_2+n_2^2{Q}_2-K{Q}_1=F_2\left(t\right),\end{array}\right.$(1)

где $n_1^2иn_2^2$ − парциальные частоты, 𝐾 > 0, функции $F_1\left(t\right)иF_2\left(t\right)$ имеют вид

$\begin{array}{l}{F}_1\left(t\right)=r_{11}\sin ({\omega }_{1}t+{\theta }_{11})+r_{12}\sin ({\omega }_{2}t+{\theta }_{12}),\\ F_2\left(t\right)=r_{21}\sin ({\omega }_{1}t+{\theta }_{21})+r_{22}\sin ({\omega }_{2}t+{\theta }_{22})\end{array}$ (2)

где $r_{ks},{\theta }_{ks},{\omega }_k\in R,k,s=\mathrm{1,2,}{\omega }_2=\alpha {\omega }_1,\alpha \in Z$. Обозначим $T=\frac{2\pi }{{\omega }_1}$.

Исследуем вынужденные колебания системы (1), когда параметры системы мало отклоняются от заданных значений и две частоты вынужденных колебаний совпадают с двумя частотами собственных колебаний (случай резонанса). В этом случае, система (1) с учетом малого параметра будет иметь вид

$\left\{\begin{array}{l}{\ddot{q}}_1+\left(n_1^2+\epsilon d_{11}\right)q_1+\left(-K+\epsilon d_{12}\right)q_2=F_1\left(t\right),\\ {\ddot{q}}_2+\left(-K+\epsilon d_{21}\right)q_1+\left(n_2^2+\epsilon d_{22}\right)q_2=F_2\left(t\right),\end{array}\right.$ (3)

Где $d_{ks}\in R,k,s=\mathrm{1,2,}\epsilon$ – малый вещественный параметр. Заметим, что система (1) получается из системы (3), если положить 𝜀 = 0.

Сформулируем задачу для системы (3): при достаточно малых вещественных 𝜀 найти T–периодическое решение $q_1(t,\epsilon ),q_2(t,\epsilon )$ системы (3), удовлетворяющее условию $q_1\left(t\mathrm{,0}\right)=Q_1\left(t\right),q_2\left(t\mathrm{,0}\right)=Q_2\left(t\right),гдеQ_1\left(t\right),Q_2\left(t\right)$ есть T–периодическое решение системы (1).

Периодичность решения $Q_1\left(t\right),Q_2\left(t\right)$ системы (1) достигается за счет подбора амплитуд $r_{ks}$ и начальных фаз ${\theta }_{ks}$ в формуле (2) для периодических функций 𝐹1(t) и 𝐹2(t).

Представление математической модели в виде системы обыкновенных дифференциальных уравнений в нормальной форме.

Делая в системе (3) замену

$\begin{array}{l}{x}_1=q_1,x_2={\dot{q}}_1=>{\dot{x}}_1=x_2,\\ x_3=q_2,x_4={\dot{q}}_2=>{\dot{x}}_3=x_4\end{array}$ (4)

получим систему обыкновенных дифференциальных уравнений в нормальной форме

$\frac{dx}{dt}=(B_0-\epsilon B_1)x-f\left(t\right),$ (5)

где $x\in R^4$,

$B_0=\left(\begin{array}{l}0100\\ -n_1^{2}0K0\\ 0001\\ K0-n_2^{2}0\end{array}\right),B_1=\left(\begin{array}{l}0000\\ d_{\mathrm{1,1}}0d_{\mathrm{1,2}}0\\ 0000\\ d_{\mathrm{2,1}}0d_{\mathrm{2,2}}0\end{array}\right),f\left(t\right)=\left(\begin{array}{l}0\\ -F_1\left(t\right)\\ 0\\ -F_2\left(t\right)\end{array}\right).$

Системе (1) соответствует неоднородная линейная система обыкновенных дифференциальных уравнений в нормальной форме

$\frac{dz}{dt}=B_{0}z-f\left(t\right).$ (6)

В этом случае задача, сформулированная для систем (1) и (3), перейдет в задачу для систем (5) и (6) в следующем виде: при достаточно малых вещественных 𝜀 найти T−периодическое решение 𝑥(𝑡,𝜀) системы (5), удовлетворяющее условию 𝑥(𝑡,0)=𝑧(𝑡), где 𝑧(𝑡) является T−периодическим решением системы (6).

В дальнейшем систему (5) будем называть возмущенной системой, а систему (6) – невозмущенной.

Вычисление периодического решения возмущенной системы методом Ляпунова- Шмидта.

В работе [1] приведено решение поставленной задачи методом Ляпунова-Шмидта. Периодическое решение возмущенной системы будем вычислять с помощью разработанного на основе метода Ляпунова-Шмидта алгоритма в математическом пакете Maple.

Для иллюстрации разработанного алгоритма в пакете Maple выберем параметры системы (5) следующим образом:

${\omega }_1=\mathrm{3,}{\omega }_2=\mathrm{9,}{n}_1^2=n_2^2=\mathrm{45,}K=\mathrm{36,}{d}_{11}=-\mathrm{2,}{d}_{22}=\mathrm{2,}{d}_{12}=d_{21}=0.$

Собственные значения матрицы 𝐵₀при выбранных значениях параметров равны

${\lambda }_{\mathrm{1,2}}=\pm 3i,{\lambda }_{\mathrm{3,4}}=\pm 9i.$

Поскольку для системы (6) частоты собственных колебаний совпадают с частотами вынужденных колебаний, то для система (6), а значит и системы (1) имеет место случай резонанса. В этом случае $T=\frac{2\pi }{3}$ и задача сводится к нахождению $\frac{2\pi }{3}$ − периодического решения системы (5).

Вычислим обобщенный жорданов набор оператора

$B_0=B_{0}x-A\frac{dx}{dt}$

по формулам

$B_0{\phi }_k^{\left(1\right)}=\mathrm{0,}{B}_0{\phi }_k^{\left(j\right)}=B_1{\phi }_k^{(j-1)},$(7)

где $j=\overline{\mathrm{2,}{p}_k},k=\overline{\mathrm{1,}n}.$

Согласно методу Эйлера, решения уравнений (7), будем искать в виде

${\phi }_k^{\left(j\right)}=u_k^{\left(j\right)}{e}^{\lambda t},$ (8)

где $\lambda \in R,u_k^{\left(j\right)}\in R^4$.

Подставляя формулы (8) в уравнения (7), получим следующие элементы обобщенных жордановых цепочек оператора B₀

${\phi }_1^{\left(1\right)}\left(t\right)=e^{i3t}\left(\begin{array}{l}-\frac{i}{3}\\ 1\\ -\frac{i}{3}\\ 1\end{array}\right),{\phi }_1^{\left(2\right)}\left(t\right)=e^{3it}\left(\begin{array}{l}-\frac{i}{54}\\ \frac{1}{18}\\ 0\\ 0\end{array}\right),$ (9)

${\phi }_2^{\left(1\right)}\left(t\right)=e^{i9t}\left(\begin{array}{l}-\frac{i}{9}\\ -1\\ -\frac{i}{9}\\ 1\end{array}\right),{\phi }_2^{\left(2\right)}\left(t\right)=e^{i9t}\left(\begin{array}{l}-\frac{i}{54}\\ \frac{1}{18}\\ 0\\ 0\end{array}\right).$ (10)

Здесь число и длины жордановых цепочек равны =2, 𝑝₁=𝑝₂=2, соответственно.

Элементы обобщенных жордановых цепочек сопряженного оператора $B_0^{*}$, удовлетворяющих условию биортогонализации, имеют вид

Коэффициенты $c_{kj}$ в разложении искомого периодического решения вычислим по формулам

${\psi }_1^{\left(1\right)}=e^{3it}\left(\begin{array}{l}81\\ 27i\\ 81\\ 27i\end{array}\right),{\psi }_1^{\left(2\right)}=e^{3it}\left(\begin{array}{l}0\\ 0\\ -\frac{9}{2}\\ -\frac{3i}{2}\end{array}\right)$ (11)

${\psi }_2^{\left(1\right)}=e^{9it}\left(\begin{array}{l}729\\ 81i\\ -729\\ -81i\end{array}\right),{\psi }_2^{\left(2\right)}=e^{9it}\left(\begin{array}{l}0\\ 0\\ -\frac{81}{2}\\ -\frac{9i}{2}\end{array}\right).$ (12)

Учитывая формулы (11) и (12), получим

$\begin{array}{l}{c}_{11}=\frac{27r_{11}{e}^{i{\theta }_{11}}}{2}+\frac{27r_{12}{e}^{i{\theta }_{12}}}{2},c_{12}=\frac{3r_{12}{e}^{i{\theta }_{12}}}{4},\\ c_{21}=\frac{81r_{21}{e}^{i{\theta }_{21}}}{2}-\frac{81r_{22}{e}^{i{\theta }_{22}}}{2},c_{22}=\frac{9r_{22}{e}^{i{\theta }_{22}}}{4}.\end{array}$ (13)

Для того, чтобы система (6) имела T-периодическое решение, необходимо и достаточно, чтобы 𝑐₁₁=0 и 𝑐₂₁=0 [1], что эквивалентно следующим условиям, соответственно,

$r_{11}{e}^{i{\theta }_{11}}+r_{12}{e}^{i{\theta }_{12}}=\mathrm{0,}$ (14)

$r_{21}{e}^{i{\theta }_{21}}+r_{22}{e}^{i{\theta }_{22}}=0$ (15)

Решения уравнений (14) и (15) можно записать в параметрической форме, соответственно,

$r_{11}=r_{12}=r_1,{\theta }_{11}={\theta }_1{\theta }_{12}={\theta }_1+\pi ,$ (16)

$r_{21}=r_{22}=r_2,{\theta }_{21}={\theta }_2{\theta }_{22}={\theta }_2,$ (17)

где  $r_k,{\theta }_k,k=\mathrm{1,2,}$ – произвольные вещественные параметры.

После подстановки формул (16) и (17) в выражения (13), получим

$c_{11}=\mathrm{0,}{c}_{12}=\frac{3r_1{e}^{i{\theta }_1}}{4},c_{21}=\mathrm{0,}{c}_{22}=-\frac{9r_2{e}^{i{\theta }_2}}{4}.$ (18)

Тогда при 𝜀 ≠ 0 вычисляя величины ${\xi }_k$ по формулам (3.3) из работы [1], получим

${\xi }_1=-\frac{3r_1{e}^{i{\theta }_1}}{4\epsilon },{\xi }_2=-\frac{9r_2{e}^{i{\theta }_2}}{4\epsilon }$ (19)

Дополнительное слагаемое y (𝜀), входящее в  $\frac{2\pi }{3}$ −периодическое решение системы (5) и принадлежащие к дополнению корневого пространства, имеет вид

$y(t,\epsilon )=\left(\begin{array}{l}-\frac{(18\epsilon +1)}{36{\epsilon }^2-36}{r}_1\sin (3t+{\theta }_1)-\frac{(18\epsilon -1)}{36{\epsilon }^2-36}{r}_2\sin (9t+{\theta }_2)\\ -\frac{(18\epsilon +1)}{12{\epsilon }^2-12}{r}_1\cos (3t+{\theta }_1)-\frac{(18\epsilon -1)}{4{\epsilon }^2-4}{r}_2\cos (9t+{\theta }_2)\\ -\frac{\epsilon }{2{\epsilon }^2-2}{r}_1\sin (3t+{\theta }_1)+\frac{\epsilon }{2{\epsilon }^2-2}{r}_2\sin (9t+{\theta }_2)\\ -\frac{3\epsilon }{2{\epsilon }^2-2}{r}_1\cos (3t+{\theta }_1)-\frac{9\epsilon }{36{\epsilon }^2-36}{r}_2\cos (9t+{\theta }_2)\end{array}\right)$ (20)

Таким образом, подставляя формулы (9), (10), (18) и (19) в выражение (3.4) из работы [1], получим $\frac{2\pi }{3}$– периодическое решение системы (5)

$x(t,\epsilon )=\left(\begin{array}{l}-\frac{1}{2\epsilon }{r}_1\sin (3t+{\theta }_1)-\frac{1}{2\epsilon }{r}_2\sin (9t+{\theta }_2)\\ -\frac{3}{2\epsilon }{r}_1\cos (3t+{\theta }_1)-\frac{9}{2\epsilon }{r}_2\cos (9t+{\theta }_2)\\ -\frac{1}{2\epsilon }{r}_1\sin (3t+{\theta }_1)+\frac{1}{2\epsilon }{r}_2\sin (9t+{\theta }_2)\\ -\frac{3}{2\epsilon }{r}_1\cos (3t+{\theta }_1)-\frac{9}{2\epsilon }{r}_2\cos (9t+{\theta }_2)\end{array}\right)$ (21)

Построение графиков компонент периодического решения и фазовых траекторий системы (3).

Учитывая замену (4) из формулу (21), получим

$\begin{array}{l}{q}_1(t,\epsilon )=-\frac{1}{2\epsilon }{r}_1\sin (3t+{\theta }_1)-\frac{1}{2\epsilon }{r}_2\sin (9t+{\theta }_2),\\ q_2(t,\epsilon )=-\frac{1}{2\epsilon }{r}_1\sin (3t+{\theta }_1)+\frac{1}{2\epsilon }{r}_2\sin (9t+{\theta }_2).\end{array}$

Для построения графиков компонент решений и фазовых траекторий системы (3) выберем следующие наборы параметров:

$\begin{array}{l}1)r_1=\mathrm{0,5};r_2=\mathrm{0,3};{\theta }_1=0;{\theta }_2=0;\\ 2)r_1=\mathrm{0,5};r_2=\mathrm{0,3};{\theta }_1=0;{\theta }_2=\frac{2\pi }{5};\\ 3)r_1=\mathrm{0,5};r_2=\mathrm{0,3};{\theta }_1=0;{\theta }_2=\frac{4\pi }{5}.\end{array}$

На рис. 1 и рис. 2. для случая 1) приведены графики компонент $\frac{2\pi }{3}$−периодического решения и фазовых траекторий системы (3) при различных значениях параметра 𝜀.

Рис. 1. Графики компонент а) ${\mathit{q}}_{\mathbf{1}}\mathbf{(}\mathit{t}\mathbf{,}\mathit{\epsilon }\mathbf{)}$ и б) ${\mathbf{q}}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{,}\mathbf{\epsilon }\mathbf{)}\frac{\mathbf{2}\mathbf{\pi }}{\mathbf{3}}$  — периодического решения системы (3) при различных значениях параметра 𝜀 в случае 1).

Рис. 2. График фазовой траектории системы (3) в конфигурационном пространстве 𝑂𝑞₁𝑞₂ при различных 𝜀 в случае 1).

На рис. 3 и рис. 4. для случая 2) приведены графики компонент $\frac{2\pi }{3}$−периодического решения и фазовых траекторий системы (3) при различных значениях параметра 𝜀.

Рис. 3. Графики компонент а) ${\mathit{q}}_{\mathbf{1}}\mathbf{(}\mathit{t}\mathbf{,}\mathit{\epsilon }\mathbf{)}$ и б) ${\mathit{q}}_{\mathbf{2}}\mathbf{(}\mathit{t}\mathbf{,}\mathit{\epsilon }\mathbf{)}\frac{\mathbf{2}\mathbf{\pi }}{\mathbf{3}}$ — периодического решения системы (3) при различных значениях параметра 𝜀 в случае 2).

Рис. 4. График фазовой траектории системы (3) в конфигурационном пространстве 𝑂𝑞₁𝑞₂ при различных 𝜀 в случае 2).

На рис. 5 и рис. 6. для случая 3) приведены графики компонент  − периодического решения и фазовых траекторий системы (3) при различных значениях параметра 𝜀.

Рис. 5. Графики компонент а)${\mathit{q}}_{\mathbf{1}}\mathbf{(}\mathit{t}\mathbf{,}\mathit{\epsilon }\mathbf{)}$ и б)${\mathit{q}}_{\mathbf{2}}\mathbf{(}\mathit{t}\mathbf{,}\mathit{\epsilon }\mathbf{)}\frac{\mathbf{2}\mathbf{\pi }}{\mathbf{3}}$  − периодического решения системы (3) при различных значениях параметра 𝜀 в случае 3).

Рис. 6. График фазовой траектории системы (3) в конфигурационном пространстве 𝑂𝑞₁𝑞₂ при различных 𝜀 в случае 3).

Заключение.

Из формулы (21) следует, что каждая компонента $\frac{2\pi }{3}$ – периодического решения системы (3) имеют полюс первого порядка в точке 𝜀 = 0 и зависит от тех же начальных фаз ${\theta }_{ks}$ и амплитуд $r_{ks}$, что и периодические функции 𝐹1(t) и 𝐹2(t).