CC BY
8
Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 41. №4. C. 146-166. ISSN 2079-6641
УДК 517.925 Научная статья
Исследование дробной динамической системы Селькова
Р. И. Паровик
Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, п. Паратунка, ул. Мирная, 7, Россия E-mail: [email protected]
Предложена дробная нелинейная динамическая система Селькова, для описания микросейсмических явлений. Эта система известна наличием автоколебательных режимов и применяется в биологии для описания гликолитических колебаний субстрата и продукта. Динамическая система Селькова также может по аналогии описать взаимодействие двух видов трещин в упруго-хрупкой среде. Первый вид - затравочные трещины с меньшей энергией, которые не регистрируются сейсмической аппаратурой, а второй тип - крупные трещины, которые порождают микросейсмы. Первый вид трещин является триггерами для трещин второго вида. Однако возможен и обратной переход. Например, когда крупные трещины теряют свою энергию и частично становятся затравочными. Далее после увеличения их концентрации процесс повторяется, обеспечивая автоколебательный характер источников микросейсм. Дробная динамическая система Селькова учитывает эффект наследственности (эредитарности) и описывается с помощью производных дробных порядков. Эредитарность колебательных систем исследуется в рамках наследственной механики и указывает на то, что динамическая система может «помнить» некоторое время, оказанное на нее воздействие, что характерно для вязкоупругих и пластичных сред. Порядки дробных производных связаны с эредитарностью системы и отвечают за интенсивность диссипации энергии, испускаемую трещинами первого и второго видов. В работе исследуется дробная динамическая модель Селькова с помощью численного метода Адамса-Башфорта-Моултона, построены осциллограммы и фазовые траектории, исследованы точки покоя. Показано, что дробная динамическая модель может обладать релаксационными и затухающими колебаниями, а также хаотическими режимами.
Ключевые слова: динамическая система Селькова, автоколебательный режим, осциллограммы, фазовые траектории, бифуркационные диаграммы, метод Адамса-Башфорта-Мултона
d DOI: 10.26117/2079-6641-2022-41-4-146-166
Поступила в редакцию: 07.12.2022 В окончательном варианте: 14.12.2022
Для цитирования. Паровик Р. И. Исследование дробной динамической системы Селькова // Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 41. № 4. C. 146-166. d DOI: 10.26117/2079-6641-2022-41-4-146-166
Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.Org/licenses/by/4.0/deed.ru)
© Паровик Р. И., 2022
Финансирование. Работа выполнена за счет средств РНФ (проект № 22-11-00064).
Введение
Микросейсмы - колебания земной поверхности малой амплитуды, источником которых являются процессы природного и техногенного характера. В частности, процессы в атмосфере (циклоны), воздействие морей и океанов на берега, процессы, связанные с человеческой деятельностью (строительство зданий и сооружений). Для возникновения микросейсм в зонах земной коры необходимо наличие сред, которые, с одной стороны, могут накапливать напряжения, что ассоциируется с упругостью, а с другой - обладают свойством хрупкости, т.е. способностью разрушаться под воздействием усилий, уровень которых заметно ниже предела текучести [1]. Хрупкость среды реализуется на макроскопическом уровне за счет возникновения и роста длины трещин. Процесс множественного трещинообразования заканчивается разрушением (утерей цельности) исследуемого объекта. В настоящей работе мы будем исследовать процесс формирования трещинообразования на ранней стадии разрушения.
Необходимо отметить, что трещины могут со временем не развиваться, а с учетом вязкости среды, некоторые трещины даже могут затягиваться. Такой процесс по мнению авторов статьи [2] может происходить в земной коре на глубинах, где существует высокое давление и температура, которые способствуют процессу диффузии у кончиков трещин, что приводит затягиванию трещин.
Микросейсмы можно разделить на два типа: микросейсмы первого и второго типа. Микросейсмы первого типа - регулярные слабые колебания с периодом от 2 до 10 с. Микросейсмы второго типа - менее регулярные с более длительным периодом колебаний до 30 с [2].
Микросейсмы первого рода возбуждаются трещинами малой длины, которые не регистрируются сейсмической аппаратурой. Такие трещины мы будем называть 1г-трещинами, они являются триггерами для более крупных трещин. Микросейсмы второго рода возбуждаются трещинами большей длины, которые уже регистрируются сейсмической аппаратурой. Эти трещины мы будем называть gs-трещинами.
Рис. 1. Автоколебательный процесс при взаимодействии между собой tr-трещин и gs-трещин
[Figure 1. Self-oscillatory process in the interaction between tr-cracks and gs-cracks]
Microbism5»
Warm
В настоящей работе мы будем исследовать механизм автоколебаний источников микросейсм или колебания концентрации gs-трещин по аналогии с работой [2]. Автоколебательный процесс здесь заключается во взаимодействии между собой ^-трещин и gs-трещин (рис. 1).
Первый тип трещин является затравочными с меньшей энергией и при достижении критического уровня концентрации переходит во второй тип. Трещины второго типа являются источником микросейсмических явлений (колебаний) и после отдачи своей энергии они частично исчезают и частично переходят в затравочные трещины (принцип Ле-Шателье-Брауна, [3]). Далее, этот автоколебательный процесс повторяется.
В статье [2] авторами был предложен интересный подход к описанию автоколебательного процесса при взаимодействии трещин в упруго-хрупкой среде, который основывается на применении нелинейной динамической системы Селькова, которая исследуется в рамках биологии [4].
В настоящей работе мы рассмотрим обобщение динамической системы Селькова на случай учета эффекта наследственности (эредитарности). Эредитарность означает, что система может некоторое время «помнить» оказанное на нее воздействие. Эффекты эредитарности впервые были рассмотрены в рамках наследственной механики для описания вязкоупругих или пластичных сред [5].
Математическое описание эредитарности основано на интегро-дифференциальных уравнениях вольтерровского типа с разностными ядрами в подынтегральных выражениях, которые иногда называют функциями памяти [6]. Функции памяти могут быть выбраны, исходя из экспериментальных данных, свойств среды или рассматриваемого процесса.
Однако, в первом приближении, обычно функции памяти выбирают степенными для того, чтобы описать эффект эредитарности, влияние которого уменьшается со временем. Степенные функции памяти дают возможность использовать математический аппарат дробного исчисления и перейти от интегро-дифференциальных уравнений к уравнениям с производными дробных порядков [7]-[9]. Динамические системы, которые описываются с помощью производных дробных порядков, называют дробными динамическими системами или динамическими системами дробного порядка [10].
В силу того, что необходимо изучить автоколебательный характер взаимодействия трещин двух видов, то главной целью настоящей работы является - установления возможности существования автоколебательных режимов в рамках дробной динамической системы Селькова (ДДСС). Более полный качественный анализ классической динамической системы Селькова (КДСС) рассмотрен в статье [11].
Некоторые сведение из теории дробного исчисления
В этом пункте мы рассмотрим основные определения из теории дробного исчисления, более детально его аспекты можно изучить в книгах [7]-[9].
Определение 1. Дробный интеграл Римана-Лиувилля порядка а
I 0V (t) =
г ( а)
x (т) dT л ——, а > 0,t > 0, (t - т)1—а
(1)
где Г(■) - гамма-функция. Дробный интеграл Римана-Лиувилля (1) имеет следующие свойства: ^х(т) = х(1), I оуовх(т) = 10а+рх(т), 1оуовх(т) = 10в10«х(т).
Определение 2. Дробная производная Герасимова-Капуто порядка а имеет вид:
1
э0> (t) = г (m-a)J
x(m) (т) dT Л -п—,0 < m — 1 < а < m,
(t — т)а+1 —m, < ,
(2)
dmx (t) dtm
,m e N.
Оператор (2) имеет следующие свойства
Э№ (1) = х (1), 10«30> (1) -
i—1 x(k) (0) tk
k=0
k!
,t>0.
Дадим определение обобщения оператора дробного дифференцирования (2) на случай переменного дробного порядка.
Определение 3. Дробная производная типа Герасимова-Капуто переменного порядка (1) имеет вид:
1
30at(t)x (t) = <|г (m—a(t))
x(m) (т) dT
(t — т)
a(t)+1—m
,0 < m — 1 < a(t) < m,
(3)
dmx (t) dtm
,m e N.
1
Замечание 1. Необходимо отметить, что существует много определений производной дробного переменного порядка [12]. Мы же выбрали определение как прямое обобщение оператора (2), поэтому мы его будем называть дробной производной типа Герасимова-Капуто. Также отметим, что для оператора (3) нарушаются свойства композиции операторов дробного интегрирования и дифференцирования переменных порядков.
Дробная динамическая система Селькова Случай постоянной а
Рассмотрим следующую нелинейную динамическую систему:
за1 х (1) =-х (1) + ау (1)+ Ьх2 М у (г) ,х (0) = х0
302у (1) = V- ау (1) -Ьх2 (1)у (1) ,у (0) = Уо.
149
где х (1) - функция, определяющая концентрацию затравочных трещин первого типа; у (1) - функция, определяющая концентрацию трещин второго вида, которые генерируют микросейсмы, 1 € [0, Т] - координата, отвечающая за текущее время процесса, Т>0 - константа, время моделирования; хо,уо,^,а,Ъ - заданные положительные константы; операторы дробного дифференцирования понимаются в смысле Герасимова-Капуто порядков 0 < а|,а2 < 1 и определятся согласно (2).
Замечание 2. Динамическая система (4) является обобщением известной динамической системы Селькова [11], которая используется в биологии для описания гликолитических колебаний субстрата и продукта, причем при значении параметров а = а2 = 1 с ней совпадает.
Замечание 3. Необходим отметить, что в ДДСС (4) для сохранения размерности обычно вводят дополнительный параметр 0, обладающий размерностью времени:
^ 1 а 1 дд2
йг 01-«1 01-«2
Мы будет полагать в дальнейших наших исследованиях 0 = 1 .
Целью исследований в этом пункте является нахождение решения ДДСС (4) при заданных значениях параметров Х0,У0,^,а,Ъ и определение условий существования автоколебаний. Прежде чем, перейти к методике решения ДДСС (4), введем некоторые определения из теории качественного анализа динамических систем для исследования асимптотической устойчивости ее точек равновесия.
Асимптотическая устойчивость точек равновесия
Определение 4. Дробная динамическая система Селькова называется соизмеримой, если порядки дробных производных равны между собой а = а2 = а, в противном случае - несоизмеримой.
Определение 5. ДДСС (4) имеет одну точку равновесия Е = (х*,у*), которая является решениям следующей системы:
- х (г) + ау (г) + Ъх2 (г) у (г) = 0, V - ау (1)- Ъх2 (1) у (1) = 0.
и имеет вид:
v
x = v'y = a+w • (5)
Определение 6. Матрицей Якоби ДДСС (4) называется матрица вида:
'-(1 - 2bxy) а + bx2 —2bxy — (а + bx2
j (x,y)=( —2bxy) a+• (6)
Заметим, что с учетом точки равновесия (5) Якобиан (6) имеет вид:
/ ^2-а а + ъv2 \
| (х*,у*)= , ^ . (7)
1( ,у ) Ш - (а+ъ^)' ()
Исследуем точку равновесия (5) на асимптотическую устойчивость для этого воспользуемся следующими теоремами, предложенными в работе [10].
Пусть т - наименьшее общее кратное знаменателей дробных порядков а = Р1/о"1 ,а2 = ^2/^2 € 7,г = 1,2. Для соизмеримой ДДСС (4) справедлива сле-
дующая теорема.
Теорема 1. Точка равновесия (5) является асимптотически устойчивой для соизмеримой ДДСС (4), если все собственные значения Лг ее матрицы Якоби (7), удовлетворяют условию.
а п
|ащ(Лг)|>—, ,1 = 1,...,2, (8)
и вычисляются с помощью характеристического уравнения:
(агад [Лта 1 ,Лта2] -1 (х*,у*)) = 0 (9)
Теорема 2. Точка равновесия (5) является асимптотически устойчивой для несоизмеримой системы (4), где аг = т ,г = 1,2, если все собственные значения Лг ее матрицы Якоби |, удовлетворяют условию.
|аг§(Л)| = 1/т, (10)
где Л вычисляются из характеристического уравнения (9).
Теорема 3. (Критерий существования хаотических режимов). Если выполняется условие:
п
М = ---|Лг1) <0, (11)
2т г
то в динамической системе отсутствуют хаотические режимы.
Методика решения
В качестве метода решения ДДСС (4) выберем численный метод Адамса-Башфорта-Моултона (АБМ), который относится к типу численных методов предиктор-корректор. Метод АБМ подробно изучен и обсужден в работах [13]-[15]. Адаптируем этот метод для решения ДДСС (4). Для этого мы воспользуемся определениями (1), (2), а также их свойствами и на равномерной сетке с шагом т = T/N введем функции ,уП+1 , n = 0,...,N — 1, которые будут определяться по формуле Адамса-Башфорта (предиктор):
хП+1 = x0 +
т
«1
Г (ai + 1)f
j=0
т a2 ' / \
уП+I = yo + Г (a2+eJn+1 (v——Wj,
[9}>n+i =(n — j +1)a — (n — j)a ,i = 1,2,
а также функции хп+1 ,уп+, которые будут определяться по формуле Адамса-Моултона для корректора:
та1 / / / \ 2 \ П ,
Т I I р р ' ' Р 1 Р 1 I ^ 1
Xn+1 = xo + г ^ + 2) ( -хП+1 + ауП+1 + Ъ (хП+1) уП+1 J + Pj,u+1 (-xj + аУ) + bx?yj
Уп+1 = yo + г (J + 2) (v- аУП+1 -b(*í+i) уП+i + ¿ Р2,п+1 (v- аУ)-Ц^) | ,
(13)
Г (а2 + 2К ^
где весовые коэффициенты в (13) определяются по формуле:
'па+1 _(п - а.)(п + 1)а1 ,3 = 0,
1 _ (п-3 + 2)а+1 +(п-])а1+1 -2(п-3 + 1)а+1,1 < 3 < п,
1, ,3 = п +1,
1 = 1,2. V '
Известно, что для метода АБМ (х1 = х (1) ,Х2 = у (1) ,1 = 1,2) справедлива оценка
I ( ) | /1 +т1 па^
погрешности: тах |х^ (3 - ху | = О ^т
Замечание 3. Заметим, что в классическом случае а! = 1, мы получаем классический метод АБМ второго порядка точности.
Результаты исследований
Пример 1. Если ввести в FSDS (4) ai = a, то мы приходим к модели КСДД [2]. Значения параметров системы ДДСС взяты из статьи [2]: а = 0,116, Ъ = 0,602,t е [0,100]. Значения параметра v = 0,64421 выбраны исходя из того, что след матрицы Якоби равен нулю, чтобы существовал предельный цикл.
В таком случае, характеристическое уравнение (9) имеет чисто мнимые корни: Á1,2 = ±0.60484371641, поэтому точка равновесия может быть центром или фокусом.
Следовательно, точка равновесия E = (0.6442125706,1.760933066) является центром. Точка равновесия типа центр устойчива, но не асимптотически устойчива. Действительно, согласно Теореме 1 для соизмеримой системы (4) при a = a = 1 условие (8) не выполняется. Однако, фазовые траектории в этом случае достигают устойчивого предельного цикла (рис. 2), как ив [2].
На рис. 2 показаны расчетные кривые осциллограмм и фазовых траекторий, полученный методом (12), (13) с учетом того, что t е [0,100] и количество вычислительных узлов N = 2000. Начальные условия (%0,У0) были выбраны следующим образом: (1,0.5) - черная кривая и (0.6,1.2) - красная кривая.
Это было сделано для того, чтобы показать устойчивость предельного цикла (рис. 2c) для Примера 1. Действительно, мы видим, что выделенная красным цветом фазовая траектория, раскручиваясь, стремится к некоторой предельной траектории.
b
10 20 30 40 50 60 70 80 90
10 20 30 40 50 60 70 80 90
X(t)
- (1,0.5) - <0.6Д.2)|
Рис. 2. Осциллограммы концентраций gs-трещин (а) и tr-трещин (b); c) фазовые траектории, построенные при различных начальных условиях для Примера 1
[Figure 2. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectories constructed under different initial conditions for Example 1]
c
С другой стороны, фазовая траектория, выделенная черным цветом, также стремится к этой предельной траектории. Эта предельная траектория является устойчивым циклом.
Рассмотрим пример в случае соизмеримой ДДСС.
Пример 2. Рассмотрим случай соизмеримый ДДСС (4) с а = а2 = 0.8,а = 1. Мы выбрали значения параметров а = 0.03,Ь = 1.3^ = 0.6, I € [0,200], количество узлов сетки N = 2000. Точка равновесия в этом случае будет иметь координаты Е = (0.6,1.204819277).
Согласно Теореме 1 характеристическое уравнение для этой точки имеет вид форма:
Л16 - 0.3815180723Л8 + 0.498 = 0,
корни которого удовлетворяют условию (8). Действительно, в Таблице 1 показаны все корни характеристического уравнения с проверкой условия (8).
Следовательно, точка равновесия асимптотически устойчива. Мера существование хаотических режимов отрицательна (10): М = —0.0050548591 < 0.
На основании вышеизложенного можно сделать вывод, что система имеет замкнутую фазовую траектории и не имеет хаотических режимов. В связи с тем, что точка покоя асимптотически устойчивая, замкнутая траектория представляет собой устойчивый предельный цикл (рис. 3с).
Таблица 1
Корни характеристического уравнения и их проверка для условия (8) [Roots of the characteristic equation and their verification for the condition (8)].
№ Корни Условие (8)
1 0.9448075810 - 0.1545424593I True
2 0.7773578684 + 0.5588018265I True
3 0.1545424593 + 0.9448075810I True
4 -0.5588018265 + 0.777357868I True
5 -0.9448075810 + 0.154542459I True
6 -0.7773578684 - 0.558801827I True
7 -0.1545424593 - 0.94480758I True
8 0.5588018265 - 0.7773578684I True
9 0.9448075810 + 0.154542459I True
10 0.5588018265 + 0.7773578684I True
11 -0.1545424593 + 0.94480758I True
12 -0.7773578684 + 0.55880183I True
13 -0.9448075810 - 0.15454246I True
14 -0.5588018265 - 0.77735787I True
15 0.1545424593 - 0.9448075810I True
16 0.7773578684 - 0.5588018265I True
_xll)
- (1,0.5) - (0.6,1.2);
Рис. 3. Осциллограммы концентраций gs-трещин (a) и tr-трещин (b); c) фазовые траектории, построенные при различных начальных условиях для Примера 2
[Figure 3. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectories constructed under different initial conditions for Example 2]
Как и в предыдущем примере, мы видим, что фазовая траектория, выделенная красным, по мере раскручивания стремится к некоторой предельной траектории. В то же время фазовая траектория, выделенная черным цветом, также стремится к этой предельной траектории.
Пример 3. Рассмотрим несоизмеримый FSDS (4) а = 1,СС2 = 0,9, остальные значения параметров выберем следующим образом а = 0,12,Ь = 1.1^ = 0.6.
Тогда точка равновесия будет иметь координаты: Е = (0.6,1.162790698) Исследуем эту точку на асимптотическую устойчивость.
Согласно Теореме 2 характеристическое уравнение (10) с учетом того, что а.1 = Ю,сс2 = 10, т = 10, имеет вид:
Л19 - 0,5348837209Л10 + 0,516Л9 + 0,516 = 0.
Все корни этого уравнения удовлетворяют условию (11) Теоремы 2, поэтому точка равновесия асимптотически устойчива (табл. 2).
Таблица 2
Корни характеристического уравнения и их проверка для условия (11) [Roots of the characteristic equation and their verification for the condition (11)].
№ Корни Условие (11)
1 0.9587107281 + 0.1596989864I True
2 0.8338678747 + 0.4505677490I True
3 0.6751185224 + 0.7318032772I True
4 0.3733655051 + 0.8484303397I True
5 0.08517791582 + 1.011916745I True
6 -0.222741495 + 0.8839001155I True
7 -0.562141402 + 0.8623046746I True
8 -0.711638255 + 0.5546773068I True
9 -0.980118554 + 0.3368792296I True
10 -0.8992016786 True
11 -0.980118554 - 0.3368792296I True
12 -0.711638255 - 0.5546773068I True
13 -0.562141402 - 0.8623046746I True
14 -0.222741495 - 0.8839001155I True
15 0.08517791582 - 1.011916745I True
16 0.3733655051 - 0.8484303397I True
17 0.6751185224 - 0.7318032772I True
18 0.8338678747 - 0.4505677490I True
19 0.9587107281 - 0.1596989864I True
На рис. 4 показаны расчетные кривые осциллограмм и фазовых траекторий. построен методом АБМ (12), (13) при N = 2000^ е [0,200]. На рис. 4а и рис. 4Ь мы видим, что амплитуда колебаний увеличивается до определенного постоянного уровня, что указывает на существование предельного цикла. Однако этот цикл, как показано на рис. 4с, неустойчив, так как существуют траектории, отталкивающиеся от него (красная кривая).
а
О 20 40 60
100 120 140 160
0.4 0.5 0.6 0.7 0.8 X(t)_
- (1,0.5) - (0.6,1.2);
0.9 1.0
Рис. 4. Осциллограммы концентраций gs-трещин (a) и tr-трещин (b); c) фазовые траектории, построенные при различных начальных условиях для Примера 3
[Figure 4. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectories constructed under different initial conditions for Example 3]
c
Максимальные показатели Ляпунова для дробного осциллятора Селькова
Анализ точек равновесия для ДДСС (4) показал, что фазовые траектории могут достигать устойчивого предельного цикла. Однако важен и вопрос о существовании хаотических режимов, который мы попытаемся исследовать, используя спектры максимальных показателей Ляпунова (МПЛ), по аналогии с работами автора [16]-[18]. С помощью алгоритма Вольфа-Беннетина (АВБ) [20],[21] и численного метода (12), (13) были построены максимальные показатели Ляпунова для различных параметров ДДСС (4).
Пример 4. Построим бифуркационные диаграммы методом АБМ с учетом значений параметров: v = 0.5, а = 0.1,Ъ = а = 1, а2 = 1,t е [0,100] ,x(0) = 0.2,y (0) = 0.3, с шагом т = 0.05.
На рис. 5 мы видим, что в спектре МПЛ присутствуют положительные значения, которые указывают на наличие хаотического режима. Проверим, действительно ли этот факт имеет место в Примере 4.
Спектр МПЛ Лтах (а-|) состоит из положительных значений в диапазоне параметра а2 € [0,0,325) и отрицательных значений из диапазона а € [0.325,1] (рис. 5).
Рис. 5. Спектр МПЛ Лтах [Figure 5. MLE spectrum Лтах (ат)|
Обратим внимание, что на рис. 5 показаны фазовые траектории для различных значений а. При а = 0,5 фазовая траектория имеет вид закручивающейся спирали (устойчивый фокус) (рис. 5а), при а = 0,9 фазовая траектория представляет собой предельный цикл (рис. 5б). Нас интересует значение а = 0.1, при котором спектр МПЛ положителен. Здесь мы имеем открытую фазовую траекторию (рис. 5с). Соответствует ли такая фазовая траектория хаотическому режиму? Воспользуемся критерием (11). В этом случае М < -0.1514163724. Это означает, что этот режим не является хаотичным. Фазовая траектория соответствует апериодическому режиму, т.е. процессу затухания без колебаний.
Здесь мы наблюдаем бифуркацию Андронова-Хопфа - переход от устойчивого фокуса к предельному циклу при а —> 1. Бифуркация Адронова-Хопфа и переход к апериодическому режиму указывают на то, что фазовые траектории при а —> 0 характеризуют диссипацию энергии в динамической системе. Здесь а порядок дробной производной связан с добротностью системы по аналогии с результатами в [19], т.е. при уменьшении значения порядка дробной производной предельный цикл переходит в устойчивый фокус, и тогда мы видим затухание без колебаний. Аналогично можно показать, что для порядка а2 эти выводы остаются в силе.
Пример 5. Приведем другие спектры МПЛ: ЛШах (а2), Лтах (Ь), Лтах (V) для первого МПЛ по методу АВБ для системы ДДСС (4):
Лтах (а): V = 0.5,а = а2 = Ь = 1,1 е [0,100] ,х(0) = 0.2,у (0) = 0.3, а е (0,10] с шагом т = 0.1.
Лтах (ог): V = 0.5, а = 0.1, а1 = Ь = 1,1 е [0,100] ,х (0) = 0.2,у (0) = 0.3, а2 е (0,1] с шагом т = 0.05.
Лтах (Ь) : V = 0.5, а = 0.1,Ь е (0,3] = аг = 1,1 е [0,100] ,х(0) = 0.2,у (0) = 0.3 с шагом т = 0.1.
Лтах (V): а = 2, а = а2 = Ь = 1,1 е [0,100] ,х (0) = 0.2,у (0) = 0.3, V е (0,4] с шагом т = 0.1.
a
b
а-
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
<
X -2.0
ее =
d
.5 2 2.5 3
О
С
X -1-4 се
3
< -И-
Рис. 6. Спектры МПЛ: а-Лтах(a); b - Лтах(а2); с-Лтах(b); d -Лтах(v) [Figure 6. Spectra MLE: а-Лтах(а); b - Лтах (a2); c-^и (b); d -Л max (v)]
С
Дробная динамическая система Селькова Случай непостоянной а
В случае непостоянной а ДДСС (4), с учетом оператора (3), перепишется в виде:
(1)х (1) = -х (1) + ау (1) + Ьх2 (1) у (1),х (0) = Х0 (1)= V- ау (1)-Ьх2 (1)у (1) ,у (0) = ус.
158
где 0 < а ("Ь),а/| (1) < 1 заданные непрерывные функции.
Замечание 4. В силу того, что операторы дробной производной и интеграла переменного порядка не обладают композицией, то применение описанной выше теории в пункте 3 затруднительно. Поэтому остановимся мы здесь на методе решения системы (14), а также на построении с его помощью осциллограмм и фазовых траекторий в зависимости от вида функций а (Ь),а1(Ь).
Методика решения
Методику решения (12), (13) согласно методу АБМ обобщим на случай переменных а (1), а-1 ("Ь).
xp -Xn+1 —
x0 + 0j,n+i (-xj + ay + bxjfyjj ,
j—0
(15)
уП+1 = У0 + к2п x 02,п+1 (V - - , 1=0
[9}>п+1 =(п -1 + 1ГП-(п - ,1 = 1,2,
а также функции хп+1 ,уп+, которые будут определяться по формуле Адамса-Моултона для корректора:
Xn+1 - X0 + sin
хП+1 + ауП+1 + b (хП+1) уП+1 ) + ü jn+i (-xj + аУ + bx2y
j—0
Уп+1 - У0 + s2n I v - ауП+1 - b (хП+1 ) УП+1 + Y- Р}>+1 (v - аУ - bxfyj
j—0
где весовые коэффициенты в (15) определяются по формуле:
(16)
Р>+1
—
fnÄin+1 -(n- oiu)(n + 1)ain - 0,
(n- j + 2)ain+1 +(n-j)ain+1 -2(n-j + 1)ain+1,1 < j < n,
к1п =
T
«ii
i - 1,2. V '
-, k2n -
T
«2 n
Г (а1п + 1)' Г (а2п +1)' 51п г (а1п + 2)' "2п Г (а2п + 2)'
Для метода АБМ (х1 = х (Ь) ,Х2 = у (Ь) ,1= 1,2) справедлива оценка погрешности
T
«Ii
,j - n + 1
S2n -
T
«2 n
minI 2,1 +mina
[22]: max |xi(t^ - xi,j| - 0(t V t ) ] ,a - {a1,a2}.
Результаты моделирования
Пример 6. С помощью метода АБМ по формулам (15), (16) построим осциллограммы и фазовые траектории. Возьмем значения параметров следующие:
v = 0.6, a = 0.03,b = 1.3,T = 300,т = 0.15. Функции порядков дробных производных:
,, ло 0.01t ,, 0.02t
a1 (t) = 0.8--t—,a2(t) = 0.9--—.
Рис. 7. Осциллограммы концентраций gs-трещин (а) и tr-трещин (b); c) фазовая траектория, построенная при начальных условиях x(0) = 1,y (0) = 0.5 для Примера 6
[Figure 7. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectory constructed under initial conditions x(0) = 1,y (0) = 0.5 for Example 6]
Из рис. 7, что в случае, если функции а1 (1) и а^И монотонно убывают, то может возникнуть регулярный режим, например, предельный цикл. Этот факт имеет важное значение при моделировании автоколебательных режимов микро-сейсм с помощью ДДСС.
Пример 7. С помощью метода АБМ по формулам (15), (16) построим осциллограммы и фазовые траектории. Возьмем значения параметров следующие: V = 0.6, а = 0.03,Ь = 1.3,Т = 100,т = 0.05. Функции порядков дробных производных:
а1 (1)= 0.08 + ^а^Н 0.19 + ^
а
b
О 10 20 30 40 50 60 70 80 90
О 10 20 30 40 50 60 70 80 90
Рис. 8. Осциллограммы концентраций gs-трещин (а) и tr-трещин (b); c) фазовая траектория при начальных условиях x(0) = 0.1,y (0) = 0.3 для Примера 7 [Figure 8. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectory constructed under initial conditions x(0) = 0.1,y (0) = 0.3 for Example 7]
a
b
10 20 30 40 50 60 70 80 90
Рис. 9. Осциллограммы концентраций gs-трещин (a) и tr-трещин (b); c) фазовая траектория при начальных условиях x(0) = 0.01,y (0) = 0.05 для Примера 8 [Figure 9. Oscillograms of concentrations of gs-cracks (a) and tr-cracks (b); c) phase trajectory constructed under initial conditions x(0) = 0.01,y (0) = 0.05 for Example 8]
c
c
Из рис. 8 видно, что в случае, когда функции a (t) и ct2(t) монотонно возрастают, то может возникнуть апериодический режим - отсутствие колебаний.
Пример 8. С помощью метода АБМ по формулам (15), (16) построим осциллограммы и фазовые траектории. Возьмем значения параметров следующие: v = 0.6, a = 0.03,b = 1.3,T = 100,т = 0.05. Функции порядков дробных производных:
1 5exp(cos(nt)) , . 1 3exp(sin(nt))
= 1--т-,a2(t) = 1--t-.
В случае, если функции a (t) и ct2(t) имеют периодический характер, то могут возникать нерегулярные режимы близкие хаотическим рис. 8.
Заключение
В настоящей работы мы показали с помощью элементов качественного анализа дробный динамических систем и численного моделирования, что ДДСС может обладать автоколебательными режимами [24]. Также мы исследовали точки покоя, регулярные и хаотические режимы ДДСС для постоянных порядков дробных производных. Показано, что может существовать устойчивый предельный цикл, также с помощью бифуркационных диаграмм МПЛ показана бифуркация Андронова-Хопфа - уничтожения предельного цикла. Хаотические режима по методу АВБ не были обнаружены, хотя МПЛ для а были положительные значения спектра. Здесь возможно, необходимо уточнить метод АВБ. Это можно сделать, заменив процедуру ортогонализации Грама-Шмидта более точной процедурой Хаусхолде-ра [23].
Хаотические режимы могут возникнуть, если порядки дробных производных зависят от времени. Однако необходимо более детальное исследование. Также можно отметить, что предельный цикл, также может возникнуть. Вопрос его устойчивости необходимо рассматривать отдельно.
Конкурирующие интересы. Конфликтов интересов в отношении авторства и публикации нет.
Авторский вклад и ответственность. Автор участвовал в написании статьи и полностью несет ответственность за предоставление окончательной версии статьи в печать.
Список литературы
1. Kearey Ph. The Encyclopedia of Solid Earth Sciences: Blackwell Sci., 1993. 722 pp.
2. Маковецкий В. И., Дудченко И. П., Закупин А. С. Автоколебательная модель источников мик-росейсм, Геосистемы переходных зон, 2017. №4(1), С. 37-46.
3. Shpielberg O., Akkermans E. LeChatelier principle for out-of-equilibrium and boundary-driven systems: Application to dynamical phase transitions, Physical review letters, 2016. vol.116, no. 24, 240603.
4. Selkov E. E. Self-oscillations in glycolysis. I. A simple kinetic model, Eur. J. Biochem., 1968. vol.4, pp. 79-86.
5. Работнов Ю. Н. Элементы наследственной механики твердого тела, 1980. 392 с.
a (t)
6. Volterra V. Sur les 'equations int 'egro-differentielles et leurs applications, Acta Mathematica, 1912. vol.35, no. 1, pp. 295-356.
7. Kilbas A. A., Srivastava H.M., Trujillo J.J. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier, 2006. 523 pp.
8. Oldham K. B., Spanier J. The fractional calculus. Theory and applications of differentiation and integration to arbitrary order. London: Academic Press, 1974. 240 pp.
9. Miller K. S., Ross B. An introduction to the fractional calculus and fractional differntial equations. New York: A Wiley-Interscience publication, 1993. 384 pp.
10. Petras I. Fractional Order Nonlinear Systems. Modeling, Analysis and Simulation. Beijing-Springer-Verlag Berlin Heidelberg: Springer, 2011..
11. Brechmann P., Rendall A. D. Dynamics of the Selkov oscillator, Mathematical Biosciences, 2018. vol.306, pp. 152-159 DOI: 10.1016/j.mbs.2018.09.012..
12. Patnaik S., Hollkamp J.P., Semperlotti F. Applications of variable-order fractional operators: A review, Proc. R. Soc. A R. Soc. Publ, 2020. vol.476, 20190498 DOI: 10.1098/rspa.2019.0498..
13. Garrappa R. Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial, Mathematics, 2018. vol.6, no. 16 D0I:10.3390/math6020016..
14. Yang C., Liu F.A computationally effective predictor-corrector method for simulating fractional-order dynamical control system, ANZIAM J.,2006. vol.47, pp. 168-184 DOI: 10.21914/anzi-amj.v47i0.1037.
15. Diethelm K, Ford N.J., Freed A.D. A predictor-corrector approach for the numerical solution of fractional differential equations, Nonlinear Dyn, 2002. vol.29, pp. 3-22 DOI: 10.1023/A:1016592219341.
16. Parovik R., Rakhmonov Z., Zunnunov R. Modeling of fracture concentration by Sel'kov fractional dynamic system, E3S Web of Conferences, 2020. vol. 196, 02018.
17. Parovik R. I. Research of the stability of some hereditary dynamic systems, Journal of Physics: Conference Series, 2018. vol. 1141, no. 1, 012079.
18. Parovik R. I. Chaotic modes of a non-linear fractional oscillator, IOP Conference Series: Materials Science and Engineering, 2020. vol. 919, no. 5, 052040.
19. Parovik R. I. Quality factor of forced oscillations of a linear fractional oscillator, Technical Physics, 2020. vol.65, no. 7, pp. 1015-1019.
20. Benettin G., Galgani L., Giorgilli A., Strelcyn J. M.Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory, Meccanica, 1980. vol. 15, no. 1, pp. 9-20.
21. Wolf A., Swift, J.B., Swinney, H.L., Vastano, J. A. Determining Lyapunov exponents from a time series, Physica D: nonlinear phenomena, 1985. vol. 16, no. 3, pp. 285-317.
22. Ma S., Xu Y., Yue W. Numerical solutions of a variable-order fractional financial system, Journal of Applied Mathematics, 2012. vol.. 2012, 417942 DOI: 10.1155/2012/417942..
23. Geist K., Parlitz U., Lauterborn W. Comparision of different methods for computing Lyapunov exponents, Prog. Theor. Phys., 1990. vol. 83, no. 5.
24. Parovik R.I. Studies of the Fractional Selkov Dynamical System for Describing the Self-Oscillatory Regime of Microseisms, Mathematics, 2022. vol. 10(22), 4208 DOI: 10.3390/math10224208..
Паровик Роман Иванович А - доктор физико-математических наук, доцент, ведущий научный сотрудник лаборатории моделирования физических процессов института космофизических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия, ОЯСГО 0000-0002-1576-1860.
Vestnik KRAUNC. Fiz.-Mat. Nauki. 2022. vol. 41. no. 4. pp. 146-166. ISSN 2079-6641
MSC 86-10 Research Article
Investigation of the Selkov fractional dynamical system
R. I. Parovik
Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, 684034, Paratunka, Mirnaya str., 7, Russia E-mail: [email protected]
A fractional nonlinear Selkov dynamic system is proposed to describe microseismic phenomena. This system is known for the presence of self-oscillatory regimes and is used in biology to describe glycolytic oscillations of the substrate and product. The Selkov dynamic system can also, by analogy, describe the interaction of two types of cracks in an elastic-brittle medium. The first type is seed cracks with less energy, which are not recorded by seismic equipment, and the second type is large cracks that generate microseisms. The first type of cracks are triggers for cracks of the second type. However, the reverse transition is also possible. For example, when large cracks lose their energy and partially become seed cracks. Further, after increasing their concentration, the process is repeated, providing the self-oscillating nature of microseismic sources. The Selkov fractional dynamical system takes into account the effect of hereditarity and is described using derivative fractional orders. The heredity of oscillatory systems is studied within the framework of hereditary mechanics and indicates that a dynamic system can "remember" some time, the impact on it, which is typical for viscoelastic and plastic media. The orders of fractional derivatives are related to the hereditarity of the system and are responsible for the intensity of energy dissipation emitted by cracks of the first and second types. In this paper, the Selkov fractional dynamic model is investigated using the Adams-Bashforth-Moulton numerical method, oscillograms and phase trajectories are constructed, and rest points are investigated. It is shown that a fractional dynamic model can have relaxation and damped oscillations, as well as chaotic modes.
Key words: Selkov dynamic system, self-oscillating mode, oscillograms, phase trajectories, bifurcation diagrams, Adams-Bashforth-Multon method
d DOI: 10.26117/2079-6641-2022-41-4-146-166
Original article submitted: 07.12.2022 Revision submitted: 14.12.2022
For citation.Parovik R.I. Investigation of the Selkov fractional dynamical systems. Vestnik KRAUNC. Fiz.-mat. nauki. 2022,41: 4,146-166. d DOI: 10.26117/2079-6641-2022-41-4-146166
Competing interests. The author declare that there is no conflicts of interest regarding authorship and publication.
Funding. The work was carried out with the financial support of the Russian Science Foundation (project No. 22-11-00064).
Contribution and Responsibility. The author is solely responsible for providing the final version of the article in print. The final version of the manuscript was approved by the author. The content is published under the terms of the Creative Commons Attribution 4-0 International License (https://creativecommons.org/licenses/by/4.0/deed.ru)
© Parovik R.I., 2022
References
Kearey Ph. The Encyclopedia of Solid Earth Sciences. Blackwell Sci., 1993. 722 p. Makovetsky V. I., Dudchenko I. P., Zakupin A. S. Auto oscillation model of microseism's sources, Geosistemy perehodnykh zon. 2017, 4 (1), 37-46 (In Russian). Shpielberg O., Akkermans E. Le Chatelier principle for out-of-equilibrium and boundary-driven systems: Application to dynamical phase transitions, Physical review letters. 2016, 116:24, 240603.
Selkov E.E. Self-oscillations in glycolysis. I. A simple kinetic model, Eur. J. Biochem. 1968, 4, 79-86.
Rabotnov Yu. N. Elementy nasledstvennoy mekhaniki tverdogo tela [Elements of hereditary solid mechanics]. Moscow: Mir, 1980, 392 (In Russian).
Volterra V. Sur les 'equations int'egro-differentielles et leurs applications, Acta Mathematica, 1912, 35:1, 295-356.
Kilbas A.A., Srivastava H.M., Trujillo J.J. Theory and Applications of Fractional Differential Equations. Amsterdam, Elsevier, 2006. 523.
Oldham K. B., Spanier J. The fractional calculus. Theory and applications of differentiation and integration to arbitrary order. London: Academic Press, 1974, 240. Miller K. S., Ross B. An introduction to the fractional calculus and fractional differntial equations. New York, A Wiley-Interscience publication, 1993. 384. Petras I. Fractional Order Nonlinear Systems. Modeling, Analysis and Simulation. - Beijing, Springer-Verlag Berlin Heidelberg, Springer, 2011.
Brechmann P., Rendall A.D. Dynamics of the Selkov oscillator, Mathematical Biosciences, 2018, 306, 152-159, DOI: 10.1016/j.mbs.2018.09.012.
Patnaik S., Hollkamp J. P., Semperlotti F. Applications of variable-order fractional operators: A review, Proc. R. Soc. A R. Soc. Publ. 2020, 476, 20190498,DOI: 10.1098/rspa.2019.0498.
Garrappa R. Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial, Mathematics, 2018, 6:16, DOI:10.3390/math6020016. Yang C., Liu F. A computationally effective predictor-corrector method for simulating fractional-order dynamical control system, ANZIAM J. 2006, 47, 168-184, DOI: 10.21914/anziamj.v47i0.1037
Diethelm K, Ford N. J., Freed A. D. A predictor-corrector approach for the numerical solution of fractional differential equations, Nonlinear Dyn., 2002, 29, 3-22, DOI: 10.1023/A:1016592219341
Parovik R., Rakhmonov Z., Zunnunov R. Modeling of fracture concentration by Sel'kov fractional dynamic system, E3S Web of Conferences. EDP Sciences, 2020, 196, 02018. Parovik R.I. Research of the stability of some hereditary dynamic systems, Journal of Physics: Conference Series. IOP Publishing, 2018, 1141:1, 012079. Parovik R.I. Chaotic modes of a non-linear fractional oscillator, IOP Conference Series: Materials Science and Engineering. IOP Publishing, 2020, 919:5, 052040.
3
4
5
6
7
8
9
[19] Parovik R. I. Quality factor of forced oscillations of a linear fractional oscillator, Technical Physics, 2020, 65:7, 1015-1019.
[20] Benettin G., Galgani L., Giorgilli A., Strelcyn J. M., Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory, Meccanica, 1980, 15:1, 9-20.
[21] Wolf A., Swift, J. B., Swinney, H. L., Vastano, J. A. Determining Lyapunov exponents from a time series, Physica D: nonlinear phenomena, 1985, 16:3, 285-317.
[22] Ma S., Xu Y., Yue W. Numerical solutions of a variable-order fractional financial system, Journal of Applied Mathematics, 2012, 2012, Article ID 417942, DOI: 10.1155/2012/417942.
[23] Geist K., Parlitz U., Lauterborn W. Comparision of different methods for computing Lyapunov exponents, Prog. Theor. Phys., 1990, 83:5.
[24] Parovik R. I. Studies of the Fractional Selkov Dynamical System for Describing the Self-Oscillatory Regime of Microseisms, Mathematics. 2022, 10(22), 4208, DOI: 10.3390/math10224208.
Parovik Roman IvanovichA - D. Sci. (Phys. & Math.), Associate Professor, Leading researcher laboratory of modeling physical processes Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, ORCID 0000-0002-1576-1860.
