Study of bifurcation diagrams of Selkov’s fractional dynamic system to describe self-oscillatory modes of microseisms

Full Text

Введение

Динамические системы играют важную роль в различных областях знаний и зачастую бывает так, что одна и также динамическая система может описывать похожие процессы, но разной природы. Это свойство динамической системы иногда называют универсальностью. Не исключением является и динамическая система Селькова. Она часто встречается в биологии при исследовании гликолитических реакций, которые имеют автоколебательные режимы [1].

Далее в статье [2] было предложено использовать динамическую систему Селькова для описания автоколебательных режимов микросейсм $—$ колебаний земной поверхности малой амплитуды, источником которых являются природные и техногенные процессы.

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

Был проведен количественный и качественный анализ динамической системы Селькова с учетом наследственности, основные результаты которого были отражены не только в статье [3], но и в других работах автора [7, 9, 16].

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

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

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

$\left\{\begin{array}{l}{\partial }_{0t}^{{\alpha }_1\left(t\right)}x\left(t\right)=-v_1\left(t\right)x\left(t\right)+w_1\left(t\right)y\left(t\right)+h_1\left(t\right)x^2\left(t\right)y\left(t\right),x\left(0\right)=x_0,\\ {\partial }_{0t}^{{\alpha }_2\left(t\right)}y\left(t\right)=v_2\left(t\right)-w_2\left(t\right)y\left(t\right)-h_2\left(t\right)x^2\left(t\right)y\left(t\right),y\left(0\right)=y_0.\end{array}\right.$ (1)

где $x\left(t\right)\mathrm{,}y\left(t\right)\in C^1\left[\mathrm{0,}T\right]$ $—$ функции решения, $v_1\left(t\right)\mathrm{<mo>=</mo>}{\theta }^{1-{\alpha }_1\left(t\right)}$, $v_2\left(t\right)\mathrm{<mo>=</mo>}{v}_0{\theta }^{1-{\alpha }_2\left(t\right)}$, $w_1\left(t\right)\mathrm{<mo>=</mo>}{w}_0{\theta }^{1-{\alpha }_1\left(t\right)}\mathrm{,}{w}_2\left(t\right)\mathrm{<mo>=</mo>}{w}_0{\theta }^{1-{\alpha }_2\left(t\right)}\mathrm{,}{h}_1\left(t\right)\mathrm{<mo>=</mo>}{h}_0{\theta }^{1-{\alpha }_1\left(t\right)}\mathrm{,}{h}_2\left(t\right)\mathrm{<mo>=</mo>}{h}_0{\theta }^{1-{\alpha }_2\left(t\right)}$ $—$ функции из класса $C\left[\mathrm{0,}T\right]$, $\theta$ $—$ параметр имеющий размерность времени, $v_0\mathrm{,}{w}_0\mathrm{,}{h}_0$ $—$ заданные константы, $t\in \left[\mathrm{0,}T\right]$ $—$ текущее время процесса, $T\mathrm{<mo>></mo>0}$ $—$ время моделирования; $x_0\mathrm{,}{y}_0$ $—$ положительные константы, отвечающие за начальные условия; операторы дробных производных имееют вид:

${\partial }_{0t}^{{\alpha }_1\left(t\right)}x\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}-{\alpha }_1\left(t\right)\mathrm{<mo\; stretchy="false">)</mo>}}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{\dot{x}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}d\tau }{{\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}}^{{\alpha }_1\left(t\right)}}\mathrm{,}{\partial }_{0t}^{{\alpha }_2\left(t\right)}y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{1}{\Gamma \mathrm{<mo\; stretchy="false">(</mo>1}-{\alpha }_2\left(t\right)\mathrm{<mo\; stretchy="false">)</mo>}}{\displaystyle \underset{0}{\overset{t}{\int }}}\frac{\dot{y}\mathrm{<mo\; stretchy="false">(</mo>}\tau \mathrm{<mo\; stretchy="false">)</mo>}d\tau }{{\mathrm{<mo\; stretchy="false">(</mo>}t-\tau \mathrm{<mo\; stretchy="false">)</mo>}}^{{\alpha }_2\left(t\right)}}\mathrm{,}$

понимются в смысле Герасимова-Капуто [11, 18], порядки которых $\mathrm{0<mo><</mo>}{\alpha }_1\left(t\right)\mathrm{,}{\alpha }_2\left(t\right)\mathrm{<mo><</mo>1}$ являются функциями из класса $C\left[\mathrm{0,}T\right]$.

Замечание 1. Отметим, что информацию о производных дробного переменного порядка можно найти в обзорной статье [12].

В настоящей работе мы будем использовать численный алгоритм, основанный на семействе методов предиктор-корректор (метод Адамса-Башфорта-Мултона) [13-16].

$\left\{\begin{array}{l}{x}_{k+1}^p\mathrm{<mo>=</mo>}{x}_0+\frac{{\tau }^{{\alpha }_{\mathrm{1,}k}}}{\Gamma \left({\alpha }_{\mathrm{1,}k}+1\right)}{\displaystyle \sum _{j\mathrm{<mo>=</mo>0}}^k}{\theta }_{j\mathrm{,}k+1}^1\left(-v_{\mathrm{1,}j}{x}_j+w_{\mathrm{1,}j}{y}_j+h_{\mathrm{1,}j}{x}_j^2{y}_j\right)\mathrm{,}\\ y_{k+1}^p\mathrm{<mo>=</mo>}{y}_0+\frac{{\tau }^{{\alpha }_{\mathrm{2,}k}}}{\Gamma \left({\alpha }_{\mathrm{2,}k}+1\right)}{\displaystyle \sum _{j\mathrm{<mo>=</mo>0}}^k}{\theta }_{j\mathrm{,}k+1}^2\left(v_{\mathrm{2,}j}-w_{\mathrm{2,}j}{y}_j-h_{\mathrm{2,}j}{x}_j^2{y}_j\right)\mathrm{,}\\ {\theta }_{j\mathrm{,}k+1}^i\mathrm{<mo>=</mo>}{\left(k-j+1\right)}^{{\alpha }_{i\mathrm{,}k}}-{\left(k-j\right)}^{{\alpha }_{i\mathrm{,}k}}\mathrm{,}i\mathrm{<mo>=</mo>1,2.}\end{array}\right.$ (2)

Для корректора (формула Адамса-Моултона) получим:

$\left\{\begin{array}{l}{x}_{k+1}=x_0+K_{\mathrm{1,}k}\left(-v_{\mathrm{1,}k+1}{x}_{k+1}^p+w_{\mathrm{1,}k+1}{y}_{k+1}^p+h_{\mathrm{1,}k+1}{x}_{k+1}^{p2}{y}_{k+1}^p\right)+\\ +K_{1.k}\left(\underset{j=0}{\overset{k}{}}{\rho }_{j,k+1}^1\left(-v_{\mathrm{1,}j}{x}_j+w_{\mathrm{1,}j}{y}_j+h_{\mathrm{1,}j}{x}_j^2{y}_j\right)\right),\\ y_{k+1=}{y}_0+K_{\mathrm{2,}k}\left(v_{\mathrm{2,}k+1}-w_{\mathrm{2,}k+1}{y}_{k+1}^p-h_{\mathrm{2,}k+1}{x}_{k+1}^{p2}{y}_{k+1}^p\right)+\\ +K_{\mathrm{2,}k}\underset{j=0}{\overset{k}{}}{\rho }_{j,k+1}^2\left(v_{\mathrm{2,}j}-w_{\mathrm{2,}j}{y}_j-h_{\mathrm{2,}j}{x}_j^2{y}_j\right)\end{array}\right.$ (3)

где $K_{\mathrm{1,}k}\mathrm{<mo>=</mo>}\frac{{\tau }^{{\alpha }_{\mathrm{1,}k}}}{\Gamma \left({\alpha }_{\mathrm{1,}k}+2\right)}\mathrm{,}{K}_{\mathrm{2,}k}\mathrm{<mo>=</mo>}\frac{{\tau }^{{\alpha }_{\mathrm{2,}k}}}{\Gamma \left({\alpha }_{\mathrm{2,}k}+2\right)}$, а весовые коэффициенты в (3) определяются по формуле:

  ${\rho }_{j,k+\text{1}}^i=\left\{\begin{array}{l}{k}^{{\alpha }_{i,k}+\text{1}}-\left(k-{\alpha }_{i,k}\right){\left(k+\text{1}\right)}^{{\alpha }_{i,k}},j=\text{0},\\ {\left(k-j+\text{2}\right)}^{{\alpha }_{i,k}+\text{1}}+{\left(k-j\right)}^{{\alpha }_{i,k}+\text{1}}-\text{2}{\left(k-j+\text{1}\right)}^{{\alpha }_{i,k}+\text{1}},\text{1}\le j\le k,\\ \text{1},j=k+\text{1,}\\ i=\text{1},\text{2}.\end{array}\right.$

Замечание 2. Исследование свойств метола Адамса-Башфорта-Мултона проводилось в статье автора [16].

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

Численный алгоритм (2), (3) был реализован на языке программирования Python [17] в среде PyCharm 2024.1 с возможностью визуализации результатов моделирования [18].

Пример 1. На рис.1 приведен график 3D поверхностей $x\left({\alpha }_1\mathrm{,}{\alpha }_2\right)$ и $y\left({\alpha }_1\mathrm{,}{\alpha }_2\right)$, где ${\alpha }_1\mathrm{,}{\alpha }_2\in \left[\mathrm{0.1,1}\right]$, $v\mathrm{<mo>=</mo>0.6,}w\mathrm{<mo>=</mo>0.03,}h\mathrm{<mo>=</mo>1.3,}\theta \mathrm{<mo>=</mo>1,}{x}_0\mathrm{<mo>=</mo>0.1}$, $y_0\mathrm{<mo>=</mo>0.1,}t\in \left[\mathrm{0,100}\right]$, $N\mathrm{<mo>=</mo>3000}$, ${\alpha }_1\mathrm{,}{\alpha }_2$ $–$ константы.

Рис. 1. Поверхности: a) x = x(α₁,α₂); b) y = y(α₁,α₂).

Figure 1. Surfaces: a) x = x(α₁,α₂); b) y = y(α₁,α₂).

На рис. 1. приведены бифуркационные диаграммы в виде поверхностей искомого решения $x$ и $y$ от значений порядков дробных производных ${\alpha }_1$ и ${\alpha }_2$. Отметим, что на поверхностях рис.1a и рис. 1b есть области которые отвечают за регулярные режимы, например, затухающие колебания соответствуют областям без всплесков, в области всплесков могут формироваться предельные циклы, а также предхаотические или хаотические режимы. Кроме того, мы видим область рваные области, что, как мы покажем дальше, связано с сингулярностью.

На рис. 2 дается бифуркационная диаграмма $—$ сечение поверхности на рис. 1 при ${\alpha }_2\mathrm{<mo>=</mo>1}$ для решения $x$ (рис.2a) и при ${\alpha }_1\mathrm{<mo>=</mo>1}$ для решения $y$ (рис. 2b). Мы видим на этих бифуркационных диаграммах, например, на рис.2a три режима, сначала идет затухающий режим вплоть до ${\alpha }_1\mathrm{<mo>=</mo>0.6}$, причем прерывистая линия вначале указывает на сингулярность. Далее идут всплески, которые указывают предельный цикл. Причем на рис.2a всплески с возрастающей амплитудой указываю на то, что орбита предельного цикла увеличивается. Это подтверждается фазовыми траекториями на врезках к рис.2a и рис.2b.

Рис. 2. Расчетные кривые a) x(α₁), α₂ = 1; b) y(α₂), α₁ = 1.

Figure 2. Calculated curves a) x(α₁), α₂ = 1; b) y(α₂), α₁ = 1.

На рис. 3 приведены бифуркационные диаграммы, построенные при других значениях параметров с врезками фазовых траекторий для различных участков диаграмм. Здесь мы можем отметить, например, на рис. 3a всплески идут с уменьшающейся амплитудой, что указывает на уменьшении орбиты предельного цикла. Здесь сингулярности нет.

Рис. 3. Расчетные кривые a) x(α₁), α₂ = 0.8; b) y(α₂), α₁ = 0.8.

Figure 3. Calculated curves a) x(α₁), α₂ = 0.8; b) y(α₂), α₁ = 0.8.

На рис. 3b всплески идут сначала с возрастающей амплитудой, потом с убывающей и т.д. Однако если такое чередование будет непоследовательным или иметь хотический характер, то мы будет приходить к хаотическим или предхаотическим режимам.

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

Рассмотрим теперь другой пример дробной динамической системы Селькова, когда ${\alpha }_1\left(t\right)$ и ${\alpha }_2\left(t\right)$ являются функциями от $t$.

Пример 2. Значения параметров выберем следующми: $N\mathrm{<mo>=</mo>10000,}t\in \left[\mathrm{0,1000}\right]$, остальные параметры возьмем из Примера 1. Порядки дробных производных изменяются во времени $t$ по следующим законам:

${\alpha }_1\left(t\right)\mathrm{<mo>=</mo>0.8}-\frac{1}{100}\cos \left(0.1\pi t\right)\mathrm{,}{\alpha }_2\left(t\right)\mathrm{<mo>=</mo>0.8}-\frac{9}{1000}\sin \left(0.1\pi t\right)\mathrm{.}$ (4)

Построим бифуркационные диаграммы в виде поверхностей для решений $x\left({\alpha }_1\mathrm{,}{\alpha }_2\right)$ и $y\left({\alpha }_1\mathrm{,}{\alpha }_2\right)$ (Рис. 4).

Рис. 4. 3D поверхности a) x(α₁,α₂); b) y(α₁,α₂).

Figure 4. 3D surfaces a) x(α₁,α₂); b) y(α₁,α₂).

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

На рис. 5 приведены расчетные кривые ${\alpha }_1\left(t\right)$ и ${\alpha }_2\left(t\right)$ по формулам (4) (рис.5a,b). Сечения поверхности плоскостями $x\left({\alpha }_1\right)$ и $y\left({\alpha }_2\right)$ (рис. 5,c,d), также фазовая траектория (рис. 5e).

Рис. 5. Расчетные кривые a) α₁ (t); b) α₂ (t); c) x(α₁); d) y(α₂); e) y = y(x).

Figure 5. Calculated curves a) α₁ (t); b) α₂ (t); c) x(α₁); d) y(α₂); e) y = y(x).

В заключение приведем бифуркационные диаграммы других ключевых параметров дробной динамической системы Селькова (рис.6).

Рис. 6. Бифуркационные диаграммы зависимостей решения x и y от различных значений параметров модели.

Figure 6. Bifurcation diagrams of the dependences of the solution x and y on various values of the model parameters.

Здесь мы также видим на бифуркационных диаграммах (рис. 6), что есть "спокойные" участки, а есть участки со всплесками. Все это указывает на наличие различных динамических режимов.

Заключение

Исследованы различные бифуркационные диаграммы для дробной динамической системы Селькова в случае, когда $\theta \mathrm{<mo>=</mo>1}$. Алгоритм построения бифуркационных диаграмм основан на численном алгоритме Адамса-Башфорта-Мултона (2), (3). Фазовые траектории и осциллограммы были получены с помощью программного комплекса ABMSelkovFracSim[16] написанного на языке программирования Python на в среде PyCharm 2014.1.

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

Дальнейшее изучение бифуркационных диаграмм связано с построением карт динамических режимов [19, 20], а также в случае когда $\theta \ne 1$. Для этих целей необходимо привлекать более мощные вычислительные ресурсы, например, вычислительные серверы с возможностью использования процессоров CPU или GPU.

About the authors

Roman I. Parovik

Author for correspondence.
Email: parovik@ikir.ru
ORCID iD: 0000-0002-1576-1860
Russian Federation, 684034, v. Paratunka, Mirnaya st., 7

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Figure 1. Surfaces: a) x = x(α1,α2); b) y = y(α1,α2).

Download (293KB)

Indexing metadata

3. Figure 2. Calculated curves a) x(α1) ,α2 = 1; b) y(α2) ,α1 = 1.

Download (123KB)

Indexing metadata

4. Figure 3. Calculated curves a) x(α1) ,α2 = 0.8; b) y(α2) ,α1 = 0.8.

Download (97KB)

Indexing metadata

5. Figure 4. 3D surfaces a) x(α1,α2); b) y(α1,α2).

Download (329KB)

Indexing metadata

6. Figure 5. Calculated curves a) α1 (t); b) α2 (t); c) x(α1); d) y(α2); e) y = y(x).

Download (87KB)

Indexing metadata

7. Figure 6. Bifurcation diagrams of the dependences of the solution x and y on various values of the model parameters.

Download (78KB)

Indexing metadata