Mathematical Model of Van der Pol-Airy Fractional Oscillator

Full Text

Введение

Разработка математических моделей дробных колебательных систем (осцилляторов) является актуальной задачей в связи с широким спектром приложений: от физических [1] до экономических [2]. Это связано с тем, что колебания могут возникать в различных системах, в том числе обладающих наследственностью, и могут иметь важные свойства, которые необходимо изучать для установления закономерностей между различными параметрами, рассматриваемой системы. Дробные осцилляторы исследуются в рамках теории дробной динамики [3], а инструментом исследования дробных осцилляторов, как правило, является аппарат математического моделирования с привлечение дробного исчисления [4].

В настоящей статье была исследована нелинейная колебательная система Ван дер Поля-Эйри с эффектами наследственности. Особенностью этой колебательной системы заключается в наличии в модельном уравнении нелинейного трения, характерного для релаксационных колебаний осциллятора Ван дер Поля [5], линейной зависимости частоты собственных колебаний от времени, характерной для колебаний с плохо затухающей амплитудой в осцилляторе Эйри [6], а также в учете эффектов наследственности, которые характерны для вязкоупругих и пластичных сред [7]. Такие эффекты можно описать с помощью интегро-дифференциальных уравнений или производных дробных порядков [4, 8, 9].

В работе с помощью явной конечно-разностной схемы первого порядка точности в рамках теории конечно-разностных схем строится численное решение. Разработана компьютерная программа в среде Maple, которая реализует численный алгоритм. С помощью, разработанной компьютерной программы были получены осциллограммы и фазовые траектории при различных значениях параметров модели. Дана интерпретация результатов моделирования. Показано, что дробный осциллятор Ван дер Поля-Эйри обладает богатой динамикой, которая заслуживает дальнейшего изучения.

Постановка задачи и методика решения

Рассмотрим следующую задачу Коши:

${\partial }_{0t}^{\alpha }x\left(t\right)+\lambda \left(a{x}^2\left(t\right)-b\right){\partial }_{0t}^{\beta }x\left(t\right)+tx\left(t\right)=\mathrm{0,}x\left(0\right)=x_0,\dot{x}\left(0\right)=y_0,$ (1)

где $x\left(t\right)\in C^2\left[\mathrm{0,}T\right]$ — функция смещения; $t\in \left[\mathrm{0,}T\right]$ — текущее время процесса; $T>0$ - время моделирования; $\lambda$ - коэффициент трения, $a\mathrm{,}b$ - константы, $x_0\mathrm{,}{y}_0$ — константы, которые определяют начальные условия.

Операторы дробного дифференцирования в уравнении (1) понимаются в смысле Герасимова-Капуто порядков $1<\alpha <2$ и $0<\beta <1$ [10, 11]:

${\partial }_{0t}^{\alpha }x\left(t\right)=\frac{1}{\Gamma \left(2-\alpha \right)}\underset{0}{\overset{t}{\int }}\frac{\ddot{x}\left(\tau \right)d\tau }{{\left(t-\tau \right)}^{\alpha -1}},{\partial }_{0t}^{\beta }x\left(t\right)=\frac{1}{\Gamma \left(1-\beta \right)}\underset{0}{\overset{t}{\int }}\frac{\dot{x}\left(\tau \right)d\tau }{{\left(t-\tau \right)}^{\beta }},$ (2)

где $\Gamma \left(\cdot \right)$ — гамма-функция Эйлера.

Замечание 1. Если в модельном уравнении (1) положить $a=\mathrm{0,}b=-1$, то мы получаем уравнение дробного осциллятора Эйри, которое было рассмотрено в статьях [12–14].

Определение. Задачу Коши (1) будем называть дробным осциллятором Ван дер Поля-Эйри.

Замечание. Дробный осциллятор Ван дер Поля-Эйри содержит не только особенности осцилляторов Ван дер Поля и Эйри, а также дробные производные (2), которые характеризуют эффекты наследственности и при значении порядков $\alpha =\mathrm{2,}\beta =1$ переходят целочисленные.

Укажем некоторые особенности, рассматриваемого осциллятора (1). Наличие нелинейного трения, характерного для осциллятора Ван дер Поля приводит к релаксационным колебаниям или автоколебаниям, которые на фазовой плоскости соответствуют предельным циклам. Релаксационные колебания Ван дер Поля находят свое приложение в физике, биологии и других науках [15]. Осциллятор Эйри обладает следующими особенностями: при отрицательных временах существует апериодический режим и происходит экспоненциальное затухание, а при положительных – наблюдаются колебания с плохо затухающей амплитудой подобно бесселевым функциям. Так как время в задаче Коши (1) положительно, то апериодических режимы будут отсутствовать. Осциллятор Эйри применяется в оптике [6], например, вблизи монохроматических каустик наблюдаются интерференционные полосы, интенсивность, которых описывается функциями Эйри, в лазерной оптике известны также пучки Эйри [16]. Дробный осциллятор Ван дер Поля и дробный осциллятор Эйри рассматривались независимо друг от друга в статьях [12–14, 17] соответственно. Было показано, что в основном при наличии дробной инерции, первый член в задаче Коши (1), происходит затухание колебаний, т.е. дробная инерция играет роль диссипативного члена как составляющей трения. Отметим, что в работах [18–21] было показано, что порядок дробной производной в инерциальном члене, а также и в диссипативном связан с качественной характеристикой колебательной системы – добротностью.

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

Введем равномерную сетку, для этого разобьём временной интервал на  равных частей с шагом $\tau =T/N$. Пусть функция решения задачи Коши (1) обладает нужной гладкостью. Тогда мы можем ввести сеточную функцию $x\left(t_k\right)=x_k,t_k=k\tau ,k=\mathrm{1,....,}N-1$, которая аппроксимирует искомую функцию $x\left(t\right)$. Отметим, что аппроксимация производных порядков дается следующими формулами [22, 23]:

${\partial }_{0t}^{\alpha }x\left(t\right)\approx A\sum _{j=0}^{k-1}{\omega }_j^{\alpha }\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right),{\partial }_{0t}^{\beta }x\left(t\right)\approx B\sum _{j=0}^{k-1}{\omega }_j^{\beta }\left(x_{k-j+1}-x_{k-j}\right),$ (3)

где $A=\frac{{\tau }^{-\alpha }}{\Gamma \left(3-\alpha \right)},B=\frac{{\tau }^{-\beta }}{\Gamma \left(2-\beta \right)},{\omega }_j^{\alpha }={\left(j+1\right)}^{2-\alpha }-j^{2-\alpha },{\omega }_j^{\beta }={\left(j+1\right)}^{1-\beta }-j^{1-\beta }$.

С учетом аппроксимаций (3) дифференциальную задачу Коши (1) мы можем записать в дискретной подстановке:

$\left\{\begin{array}{l}{x}_0=c_1,\\ x_1=c_1+\tau c_2,y_0=c_2,\\ x_{k+1}=\frac{2A+\lambda \left(a{x}_k^2-b\right)B-\tau k}{A+\lambda \left(a{x}_k^2-b\right)B}{x}_k-\frac{A}{A+\lambda \left(a{x}_k^2-b\right)B}{x}_{k-1}-\\ -\frac{A{\displaystyle \sum _{j=1}^{k-1}}{\omega }_j^{\alpha }\left(x_{k-j+1}-2x_{k-j}+x_{k-j-1}\right)}{A+\lambda \left(a{x}_k^2-b\right)B}-\frac{\lambda \left(a{x}_k^2-b\right)B{\displaystyle \sum _{j=1}^{k-1}}{\omega }_j^{\beta }\left(x_{k-j+1}-x_{k-j}\right)}{A+\lambda \left(a{x}_k^2-b\right)B},\\ k=\mathrm{2,...,}N-1.\end{array}\right.$ (4)

Задача (4) представляет собой нелокальную явную конечно-разностную схему. Можно показать, что схема (4) условно устойчива и сходится с первым порядком. Покажем с помощью вычислительной точности, полученной по правилу Рунге, что при увеличении узлов расчетной сетки N вычислительная точность стремиться к теоретической, т.е. к единице. Вычислительная точность p определяется по формуле:

$p={\log }_2\left({\xi }_j/{\xi }_{j+1}\right),$ (5)

где $\xi =\underset{i}{max}\left|x_i-x_{2i}\right|$ — ошибка, $x_i$ — численное решение вычисленное на шаге $\tau$, $x_2$ — численное решение вычисленное на шаге $\tau /2$, j — индекс, число итераций.

Пример 1. Классический осциллятор Ван дер Поля-Эйри ($\alpha =\mathrm{2,}\beta =1$). Значения параметров модели (1): $\lambda =\mathrm{0.15,}a=b=\mathrm{1,}{x}_0=\mathrm{0.2,}{y}_0=0.3$, $T=1.$

Результаты вычисления вычислительной точности схемы (4) по формуле (5) приведены в табл. 1.

 

Таблица 1

Классический осциллятор Ван дер Поля-Эйри

[Van der Pol-Airy classical oscillator]

N	$\mathit{\tau }$	$\mathit{\xi }$	p
10	1/10	0.0006974210	—
20	1/20	0.0004640763	0.5876677840
40	1/40	0.0002628898	0.8199038548
80	1/80	0.0001392983	0.9162805155
160	1/160	0.0000712656	0.9668998922

 

Из табл. 1 видно, что для Примера 1 при увеличении узлов расчетной сетки N вычислительная точность стремиться к единице.

Пример 2. Дробный осциллятор Ван дер Пол-Эйри ($\alpha =\mathrm{1.9,}\beta =0.8$). Остальные значения параметров возьмем из предыдущего примера.

Из табл. 2 мы также видим, что при увеличении узлов расчетной сетки вычислительная точность стремиться к единице.

 

Таблица 2

Дробный осциллятор Ван дер Поля-Эйри

[Van der Pol-Airy fractional oscillator]

N	$\mathit{\tau }$	$\mathit{\xi }$	p
10	1/10	0.0003115359	—
20	1/20	0.0002240659	0.4754753168
40	1/40	0.0001633936	0.4555716297
80	1/80	0.0001101243	0.5692186263
160	1/160	0.0000733393	0.5864744480
320	1/320	0.0000379055	0.9521793008

 

Результаты исследования

Рассмотрим некоторые примеры применения нелокальной явной конечно-разностной схемы (4), которая была реализована в компьютерной программе на языке Maple.

Пример 3. Классический осциллятор Ван дер Поля-Эйри ($\alpha =\mathrm{2,}\beta =1$). Значения остальных параметров модели (1): $\lambda =\mathrm{0.15,}a=\mathrm{10,}b=\mathrm{20,}t\in \left[\mathrm{0,100}\right]$, $N=2500$.

Результаты моделирования по нелокальной явной конечно-разностной схеме (4) приведена на рис. 1 при различных значениях начальных условия: $\left(\mathrm{4,0}\right)$ и $\left(\mathrm{0.1,0.1}\right)$.

 

Рис. 1. Осциллограммы a) и фазовые траектории b) для классического осциллятора Ван дер Поля-Эйлера, построенные при различных начальных условиях: $\left(\mathrm{4,0}\right)$ — серая кривая и $\left(\mathrm{0.1,0.1}\right)$ — красная кривая; с) более детальное построение фазовых траекторий.

Figure 1. Oscillograms a) and phase trajectories b) for the classical Van der Pol-Airy oscillator, plotted under different initial conditions: $\left(\mathrm{4,0}\right)$ — gray curve and $\left(\mathrm{0.1,0.1}\right)$ — red curve; c) more detailed construction of phase trajectories.

 

Из рис. 1 мы видим, что в классическом случае для Примера 1 наблюдаются релаксационные колебания, на фазовая траектория выходит на предельный цикл, который является устойчивым аттрактором.

Определение 2. Предельный цикл называется устойчивым, если все фазовые траектории, начинающиеся в некоторой окрестности, асимптотически приближаются к предельному циклу при $t\to \infty$.

Определение 3. Предельные циклы, у которых близкие фазовые траектории неограниченно к ним приближаются называются аттракторами. Предельные циклы, у которых близкие фазовые траектории от них отталкиваются называются репеллерами.

Рассмотрим другой пример, когда $a=\mathrm{0.1,}b=1$ и значения начальных условий $\left(\mathrm{8,0}\right)$ и $\left(\mathrm{0.1,0.1}\right)$, а остальные значения параметров оставим без изменения. Результаты моделирования приведены на рис. 2.

 

Рис. 2. Осциллограммы a) и фазовые траектории b) для классического осциллятора Ван дер Поля-Эйлера, построенные при $a=\mathrm{0.1,}b=1$ и различных начальных условиях: $\left(\mathrm{8,0}\right)$ — серая кривая и $\left(\mathrm{0.1,0.1}\right)$ — красная кривая; с) более детальное построение фазовых траекторий.

Figure 2. Oscillograms a) and phase trajectories b) for the classical van der Pol-Euler oscillator, plotted for $a=\mathrm{0.1,}b=1$ and different initial conditions: $\left(\mathrm{8,0}\right)$ — gray curve and $\left(\mathrm{0.1,0.1}\right)$ — red curve; c) more detailed construction of phase trajectories.

 

Здесь мы видим, что предельные циклы при различных значениях являются различными, что говорит о неустойчивости предельного цикла.

На рис.2а может сложится впечатление, что при более больших временах все-таки предельный цикл станет устойчивым. Однако это не так, увеличим время моделирования $T=200$, остальные значения оставим без изменения (рис. 3).

 

Рис. 3. Осциллограммы a) и фазовые траектории b) для классического осциллятора Ван дер Поля-Эйлера, построенные при $a=\mathrm{0.1,}b=1$, $T=200$ и различных начальных условиях: $\left(\mathrm{8,0}\right)$ — серая кривая и $\left(\mathrm{0.1,0.1}\right)$ - красная кривая; с) более детальное построение фазовых траекторий.

Figure 3. Oscillograms a) and phase trajectories b) for the classical Van der Pol-Euler oscillator, constructed at $a=\mathrm{0.1,}b=1$, $T=200$ and different initial conditions: $\left(\mathrm{8,0}\right)$ — gray curve and $\left(\mathrm{0.1,0.1}\right)$ — red curve; c) more detailed construction of phase trajectories.

 

На рис. 3 мы видим детерминированный хаотический режим, который характерен для нелинейных колебательных систем и вызван чувствительностью к изменению начальных условий. Хаотические режимы заслуживают отдельного изучения, которое опирается на известные качественные методы нелинейной динамики [24].

Рассмотрим пример дробного осциллятора Ван дер поля-Эйри.

Пример 4. Дробный осциллятора Ван дер Поля-Эйри ($\alpha =\mathrm{1.8,}\beta =0.85$). Значения параметров: $\lambda =\mathrm{0.15,}a=\mathrm{10,}b=\mathrm{20,}t\in \left[\mathrm{0,100}\right],N=2500$, начальные условия $\left(\mathrm{0.2,0}\right)$ и $\left(\mathrm{4,0}\right)$.

Результаты моделирования приведены на рис.4.

 

Рис. 4. Осциллограммы a) при начальных условиях $\left(\mathrm{0.2,0}\right)$ и $\left(\mathrm{4,0}\right)$ и фазовая траектория для начального условия $\left(\mathrm{0.2,0}\right)$.

Figure 4. Oscillograms a) for initial conditions $\left(\mathrm{0.2,0}\right)$ and $\left(\mathrm{4,0}\right)$ and phase trajectory for the initial condition $\left(\mathrm{0.2,0}\right)$.

 

На рис.4a приведены осциллограммы при начальных условиях $\left(\mathrm{0.2,0}\right)$ — кривая и $\left(\mathrm{4,0}\right)$ — красная кривая. Здесь мы видим, что предельный цикл на рис.4b является неустойчивым. Красная кривая представляет собой апериодический режим (колебания отсутствуют).

 

Рис. 5. Осциллограмма a) и фазовая траектория b) для начального условия $\left(\mathrm{0.2,0}\right)$.

Figure 5. Oscillogram a) and phase trajectory b) for the initial condition $\left(\mathrm{0.2,0}\right)$

 

На рис. 5 приведена фазовая траектория и осциллограмма, полученная по формуле (4) в зависимости от значений $\alpha =1.9$ и $\beta =0.1$, остальные значения параметров оставлены без изменения. Здесь релаксационных колебаний мы не видим, зато видим затухающие колебания.

В заключении нашего исследования посмотрим наличие хаотического режима для дробного осциллятора Ван дер Поля-Эйри. Для этого мы увеличим время моделирования, как в Примере 1, $T=160$, $\alpha =\mathrm{1.8,}\beta =0.85$ остальные параметры оставим без изменения. Результаты моделирования приведены на рис.6.

 

Рис. 6. Осциллограммы a) и фазовые траектории b) построенные по формуле (4) при различных начальных условиях $\left(\mathrm{4,0}\right)$ — красная кривая и $\left(\mathrm{0.2,0}\right)$ — черная кривая.

Figure 6. Oscillograms a) and phase trajectories b) constructed according to formula (4) under different initial conditions $\left(\mathrm{4,0}\right)$ — red curve and $\left(\mathrm{0.2,0}\right)$ — black curve.

 

На рис. 6 мы видим несколько интересных режимов. Хаотический режим наблюдается при начальном условии $\left(\mathrm{0.2,0}\right)$ мы видим здесь перемежаемость — чередования квазипериодических колебаний с хаотическими. Другой режим при начальном условии $\left(\mathrm{4,0}\right)$, характеризует бифуркацию Адронова-Хопфа — рождение и уничтожение предельного цикла.

Заключение

В статье была предложена новая математическая модель нелинейного дробного осциллятора Ван дер Поля-Эйри. Для которой был разработан численный алгоритм на основе нелокальной явной конечно-разностной схемы первого порядка точности. Численный алгоритм был реализован в компьютерной программе на языке Maple. С помощью компьютерной программы были построены осциллограммы и фазовые траектории при различных значениях параметров модели. Показано, наличие различных колебательных режимов: релаксационных, затухающих, хаотичных. Более детальный анализ режимов дробного осциллятора Ван дер Поля-Эйри в зависимости от значений его ключевых параметров может дать построение карт динамических режимов [25].

About the authors

Asal I. Salimova

National University of Uzbekistan named after Mirzo Ulugbek

Email: parovik@ikir.ru
ORCID iD: 0009-0003-9945-0991

1st year master’s student "Applied Mathematics"

Uzbekistan, 100174, Tashkent, Universitetskaya str., 4

Roman I. Parovik

Institute of Cosmophysical Research and Radio Wave Propagation, Far Eastern Branch of the Russian Academy of Sciences

Author for correspondence.
Email: parovik@ikir.ru
ORCID iD: 0000-0002-1576-1860

D. Sci. (Phys. & Math.), Associate Professor, Leading researcher, laboratory of modeling physical processes

Russian Federation, 684034, Paratunka, Mirnaya str., 7

References

Supplementary files