CC BY
4
Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 40. №3. C. 179-198. ISSN 2079-6641
УДК 517.938 Научная статья
Неявная конечно-разностная схема для осциллятора Дуффинга с производной переменного дробного порядка типа
Римана-Лиувилля
В. А. Ким1, Р. И. Паровик1'2
1 Камчатский государственный университете имени Витуса Беринга, 683032, г. Петропавловск-Камчатский, ул. Пограничная, 4, Россия
2 Институт космофизических исследований и распространения радиоволн ДВО РАН 684034, п. Паратунка, Камчатский край, ул. Мирная, 7, Россия
E-mail: valentinekim@mail.ru, romanparovik@gmail.com
В статье рассматривается неявная конечно-разностная схема для уравнения Дуффинга с производной дробного переменного порядка типа Римана-Лиувилля. Рассматриваются вопросы устойчивости и сходимости неявной конечно-разностной схемы. Для обоснования теоретических результатов приводятся тестовые примеры. С помощью правила Рунге сравниваются результаты работы неявной схемы с результатами явной схемы. Построены фазовые траектории и осциллограммы для осциллятора Дуффинга с дробной производной переменного порядка типа Римана-Лиувилля, с помощью спектра максимальных показателей Ляпунова и сечений Пуанкаре детектируются хаотические режимы. Построены поверхности добротности, амплитудно-частотной и фазо-частотной характеристик для исследования вынужденных колебаний. Результаты исследования показали, что неявная конечно-разностная схема показывает более точные результаты, чем явная.
Ключевые слова: осциллятор Дуффинга, правило Рунге, оператор Римана-Лиувилля, оператор Грюнвальда-Летникова, амплитудно-частотная характеристика, фазо-частотная характеристика, добротность, показатели Ляпунова, сечения Пуанкаре, фазовая траектория, осциллограмма.
d DOI: 10.26117/2079-6641-2022-40-3-179-198
Поступила в редакцию: 24.11.2022 В окончательном варианте: 05.12.2022
Для цитирования. Ким В. А., Паровик Р. И. Неявная конечно-разностная схема для осциллятора Дуффинга с производной переменного дробного порядка типа Римана-Лиувилля // Вестник КРАУНЦ. Физ.-мат. науки. 2022. Т. 40. № 3. C. 179-198. d DOI: 10.26117/2079-6641-2022-40-3-179-198
Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru)
© Ким В. А., Паровик Р. И., 2022
Финансирование. Финансовая поддержка выполнена в рамках гранта президента РФ "Развитие математических моделей дробной динамики с целью исследования колебательных процессов и процессов с насыщением № МД-758.2022.1.1
Введение
Осциллятор Дуффинга с дробной производной в диссипативном слагаемом имеет большое значение в решении прикладных задач математики [1], физики [2], теории хаоса [3]. Уравнение Дуффинга описывает нелинейные колебательные процессы, для которых характерны бистабильность и наличие хаотической динамики. В мировой практике бистабильность нелинейных колебаний представляет особый интерес в оптических технологиях [4], электросетях и т.д. Выявление хаотических режимов является одной из главных задач, например, в клинической медицине [5], исследовании турбулентности [3]. Существует множество методов решения дробного уравнения Дуффинга. В работе [6] рассматривался метод гомотопического анализа, метод конечно-разностных схем [7] в случае постоянного порядка дробной производной. Также в последнее время широкий интерес представляет качественный анализ дробного уравнения Дуффинга [8]-[14].
Введение дробной производной переменного порядка в уравнение Дуффинга, позволит еще гибче описывать нелинейные колебания с эффектами памяти и хаотическими режимами. Решение уравнения Дуффинга с дробной производной переменного порядка, которое было найдено с помощью явной конечно-разностной схемы, подробно рассмотрено в работе [15]. В [16] с помощью спектра максимальных показателей Ляпунова и сечений Пуанкаре исследованы хаотические и регулярные режимы. Для явной схемы теоретически обоснованы вопросы устойчивости и сходимости [17]. В работах [18] и [19] с помощью амплитудно-частотной (АЧХ), фазо-частотной характеристик (ФЧХ) и добротности исследуются свойства вынужденных колебаний осциллятора Дуффинга с дробной производной переменного порядка типа Римана-Лиувилля. Выяснилось, что порядок дробной производной влияет на скорость затухания колебаний.
Однако, точность вычислений по явной схеме не высокая. Для повышения точности вычислений и уменьшения погрешности в настоящей работе используется неявная конечно-разностная схема. Более того, в отличие от явной схемы, устойчивость и сходимость неявной не зависит от ограничений на шаг расчетной сетки. В данной статье, по аналогии с работами [15]-[19] исследуется неявная конечно-разностная схема для решения уравнения Дуффинга с дробной производной переменного порядка типа Римана-Лиувилля, обосновываются вопросы устойчивости и сходимости численной схемы, исследуются хаотические режимы и бистабильность колебаний.
Постановка задачи
Рассмотрим следующую задачу Коши для нелинейного осцилляционного уравнения по аналогии с работой [17]:
х М+ ЛЭ^^О^х м + ^0х (1)+ f (х,-Ь) = 0,х(0) = х0,х(0) = у0, (1)
где х (1) € С2 (0,Т) - функция смещения, Л - коэффициент вязкого трения, ^0 -собственная частота, х0,у0 - заданные константы, которые определяют начальные условия, 0 - параметр, отвечающий за масштаб по времени. В уравнении (1) оператор дробной производной берется в смысле Римана-Лиувилля.
Определение 1. Дробная производная Римана-Лиувилля переменного порядка 0 < < 1, функции х("Ь) € С1 [0, Т], Т > 0 имеет вид:
Dq«x (t)-_1_—
u0t х (t) = г (1 _ ^ (t)) dtj
x (т) dT (2)
(t _ т)
q(t)'
где Г (у +1) = |ху е-Мх - гамма-функция Эйлера.
0
Нелинейная функция ^х,"Ь) удовлетворяет условию Липшица относительно х.
№1 (1),1)- Пх2(-Ь),-Ь)|| <Ь|х1(1)-х2(1)|, (3)
где Ь > 0 - постоянная Липшица.
Замечание 1. Задача Коши (1) представляет собой математическую модель, описывающую широкий класс нелинейных дробных осцилляторов, вид и тип которых определяется функцией ^х,"Ь). Первое слагаемое описывает силу инерции, второе силу трения (демпфирование), а третье возвращающую силу осциллятора.
Замечание 2. Если в уравнении (1) ^х,"Ь) = Ьх3 — 6 соэ(ш"Ь), где Ь - коэффициент нелинейности, 6 и ш - амплитуда и частота внешней силы, то мы получаем уравнение осциллятора Дуффинга с дробной производной переменного порядка [17]. В дальнейшем в уравнении (1) будем считать 0 = 1.
Численный алгоритм
Так как уравнение (1) является нелинейным, то его решение ищется с помощью конечно-разностных схем. Введем равномерную расчетную сетку. Отрезок [0,Т] разобьём на N равных частей с шагом Н = N : 0 = "0 < "1 = Н<"2 = 2Н < ... < ^ = кН< ... < "м = Т. Функции q(t),x(t),f(x(t),t) перейдут в сеточные qk = q(tk),xk = х("к)Дк = ^х^),^). Аппроксимация второй производной дает:
•■ и\ хк+1 — 2хк + хк—1 , г>гг,2л (А\
хи) =-Н2-+ 0(Н ), (4)
Дробная производная переменного порядка аппроксимируется оператором Грюнвальда-Летникова (5).
Определение 2. Дробная производная Грюнвальда-Летникова переменного порядка 0 < qk < 1 имеет вид:
к+1
¿оЙ! Хк+1 = Г е?к+1 хк_1+1 + О(Н),
/ . 1=0 (5) с°°к+1 = ^ _ с0-!+1 ,Со°к+1 = 1.
Здесь с?к - весовые коэффициенты Грюнвальда-Летникова. В работе [17] была доказана следующая Лемма для весовых коэффициентов.
Лемма 1. Весовые коэффициенты с?к обладают следующими свойствами:
с0°к = 1,сОк =_Ок,с0к <0(1 = 1)
оо I (6)
х с?к = 0,VI = 1,2,..., £ с?к < 0 . (6)
1=0 1=1
В работе [17] для решения задачи коши 1 была построена нелокальная явная конечно-разностная схема (ЯКРС) для к = 1,2,...,N _ 1:
Хк+1 = АкХк _ Хк_1 _ Вк^ с°к Хк_1 _ (7)
>к 1=1
где Ак = 2 _ЛН2_°к _Н2^0,Вк = ЛН2_°к. Были доказаны следующие теоремы Теорема 1. Явная схема (7) устойчива, если выполнено условие
ЛЙ2_р + Н2^2 < 1.
Теорема 2. Явная схема (7) сходится к точному решению с первым порядком, если выполнено условие Лй2_р + Н2^0 < 1.
В теоремах 1 и 2 = тах(ок).
к
В настоящей статье, с учетом соотношений 4 и 5 построим нелокальную неявную конечно-разностную схему (НКРС) и рассмотрим вопросы ее устойчивости и сходимости
к+1
(1 + Н2^0)хк+1 _ 2хк + Хк_1 + Вк+1 ^ сг°к+1 Хк_1+1 _ нЧ+1 = 0. (8)
1=0
где Вк+1 = ЛН2_°к+1.
Замечание 3. Для построения конечно-разностной схемы функция смещения должна рассматриваться в третьем классе гладких функций х (1) € С3 [0, Т]
Устойчивость и сходимость НКРС
Определение 3. Разностная аппроксимация (8) устойчива, если для любого вектора ошибок между точным и численным решением Е0 существует положительное, число Рк : Пт Рк = 0, и выполнено условие:
к
Рк+1|1< Рк Р0Н. (9)
Теорема 3. Неявная конечно-разностная схема (8) безусловно устойчива.
Доказательство. Пусть ошибка е^ = х^ — = 0,...,М где х^ - приближенное решение задачи Коши (1). Тогда уравнение (8) в терминах ошибки примет вид:
(1 + Н2^2) £к+1 =
к+1
2ек — ек-1 — Бк+1 е^—3+1 — Н2 (f (Xk+1,tk+1) — f (xk+1,tk+1)),
1=0
k = 1,...,М — 1. (10)
Введем норму 1^+1^ = max|еk+1|. Перейдем в (10) к абсолютной величине,
k
получим:
(l + h2^0) kk+i I < 2 Ы _ |ek_il _
-г
k+1
-Бk+1 ^ | е^1 — Н2 (^ (Xk+1, tk+1) — f (xk+1, tk+1) |) <
3=0
k+1
< |ек| — Бk+1 ^ е^^1 |ем+11 — Н2Ь |ек+11
3=0
Перейдем к норме. В силу Леммы 1, соотношения (9), условия Липшица (3) и Ь > 0, мы получаем оценку
k+1
(1 + h2^2 + H2l) ||Ek+i II«, < ||Ek^_Bk+1 ^cqk+1 ||Ek_j+iL < IIEkll
j-0
! ^+1 < - 2 2 2 \ < --—-Г—к ||Е0||00.
О + Н2^0+Н2Ч (1 + н2^2+Н2^
т.е. при ^^ оо, 11~> 0. Теорема доказана. □
Пусть х - точное решение задачи Коши (1) в точке Определим = х^)— Xk и соответственно вектор Yk = ("Пь-^"^) . Заметим, что У0 нулевой вектор. Подставим Xk = х(tk)— Лk в уравнение (8), получим:
+ Лk+1 =
k+1
2лk—Лk—1 — Bk+1 Лk-j+1 — н2 ^ (х (tk+1) ,tk+1) — f (xk+1,tk+1)) +
3=0
+ Н^+ь (11)
Здесь
1^с+11 < СН,
где С - константа, не зависящая от шага Н расчетной сетки. Справедлива следующая теорема.
оо
Теорема 4. Неявная конечно-разностная схема (8) безусловно безусловно сходится к точному решению с первым порядком.
Доказательство. Перейдем в (11) к абсолютной величине, получим:
(l + h2w0J |Лк+1| < 2 |nkl - Ink—11 -
k+!
_Вк+^ с?к+1 |пк_э+11 _ й2 (|f (х (1к+1), 1к+1) _ f (хк+1, 1к+1) |) +
1=0
к+1
+Н2СН < 2 |Лк1 _ 1Пк_11 _ Вк+1 ^ с?к+1 |Лк_1+11 _ Н2Ь |Лк+11 + й3 .
1=0 к+1
(1 + Н2^2) |Ук+1^<|Ук^_Вк+1 ^с?к+1 ||Ук_5+1^ _Н21|Ук+1^ + СН3
1=0
Перейдем к норме. В силу Леммы 1, соотношения (9), условия Липшица (3) и Ь > 0, мы получаем оценку
11Ук+1 И«, < (^н2 V н2-п НУк^ + СЙ< --—2—--к||У0^+
(1 + +н2Ч (1 + н2^2+н2Тк
+ 1-1-к +--1-кт +... +1 ) СН,
\(1 + Н2^2 + Н2Ь)к (1 + Н2^2 + Н2Ь) )
т.е. при оо, ||Ук+1 СН. Теорема доказана. □
Результаты исследований
Для подтверждения теоретических результатов рассмотрим тестовые примеры. Пусть в модельном уравнении (1) нелинейная функция имеет вид ^х,"Ь) =
а ( Г(4)"4_о(") \
Ьх3 ("Ь) _ М9 _ ^0"Ь3 _ 6" _ Л— —-—т^ . Тогда точное решение новой задачи
а^Г (5 _ яи)) )
Коши находится в виде:
х(")= "3. (12)
Погрешность НКРС (8) находим с по формуле:
£ = тах(|хМш_хм[]]|) ,1 = 0,...,N, (13)
где хМШ - точное решение (14), хмШ - численное решение, полученное по схеме (8). Если точное решение неизвестно, то используем правило Рунге:
£ = тах(|хз_х^|),] = 0,...,^ (14)
где х2^ - численное решение на шаге Н, х^ - численное решение на шаге Н/2. Вычислительная точность (8) определяется по формуле:
In
а =
У ^
1п(2)
(15)
где ег - ошибка при шаге Н/2г, ег+1 - ошибка при шаге Н/2г+1, г = 0,1,..., М — 1.
Сравним результаты, полученные с помощью ЯКРС (7) и НКРС(8).
Пример 1. Рассмотрим случай когда для ЯКРС условия теорем 1 и 2 соблюдаются. В модельном уравнении (1) выберем следующие параметры: t € [0,2], х(0) = х(0) = 0,Л = 1,ш0 = Ь = 6 = 1,ш = 2, q(t) = 0.8соз(0.5"), ЛН2-р + Н2^2 = 0.908, р = тах^^к))
Рис. 1. Тестовый пример. Решения задачи Коши (1), полученные по схемам (7) и
(8), а также точное решение (12). [Figure 1. Test case. Solutions to the Cauchy problem (1) obtained by the schemes (7) and (8), as well as the exact solution (12).]
Рис. 2. а) погрешность e; б) вычислительная точность а. [Figure 2. a) e error; b) computational accuracy а.]
По рис.1 и 2 видно, что ЯКРС и НКРС достаточно хорошо аппроксимируют точное решение (12). Однако погрешность, изображенная на рис.2а, для НКРС меньше, чем погрешность для ЯКРС. Это значит, что НКРС (8) показывает более
точные результаты, чем ЯКРС (7). Вычислительная точность (рис.2б) для ЯКРС и НКРС принимает значения близкие к 1 при увеличении узлов сетки, что говорит о первом порядке сходимости ЯКРС (7) и НКРС (8) к точному решению (12).
Пример 2. Рассмотрим случай когда для ЯКРС условия теорем 1 и 2 нарушаются. В модельном уравнении (1) выберем следующие параметры: 1 € [0,2], х(0) = х(0) = 0,Л = 3,^0 = 10,Ь = 6 = 1,ш = 2,я(1") = 0.8соз(0.51"), ЛН2_Р + Н2^2 = 1.2, р = тах^"*)).
Рис. 3. Тестовый пример. Решения задачи Коши (1), полученные по схемам (7) и (8), а также точное решение (12). Для ЯКРС (7) нарушаются условия теорем 1 и 2.
[Figure 3. Test case. Solutions to the Cauchy problem (1) obtained by the schemes (7) and (8), as well as the exact solution (12). For EFDS (7), the conditions of Theorems 1 and 2 are violated.]
Рис. 4. а) погрешность e; б) вычислительная точность а. [Figure 4. a) e error; b) computational accuracy а.]
При нарушении условий теорем 1 и 2 ЯКРС (7) расходится (рис.3). В то же время НКРС (8) аппроксимирует точное решение с достаточно большой точностью. По рис.4а видно, что на промежутке нарушения условий теорем 1 и 2 погрешность ЯКРС принимает достаточно большие значения и меняется скачкообразно,
а вычислительная точность ЯКРС на этом промежутке принимает значения, сильно отличающиеся от 1. Для НКРС вычислительная точность принимает значения близкие к 1. Все это говорит о том, что устойчивость и сходимость НКРС (8) не зависит от ограничений на шаг и подтверждает справедливость теорем 3 и 4 о безусловной устойчивости и сходимости.
Сравним результаты ЯКРС (7) и НКРС (8) для дробного осциллятора Дуф-финга. Нелинейная функция ^хД) берется в форме, представленной в Замечании 1. В силу того, что осциллятор Дуффинга не имеет точного решения, то погрешность вычисляется по правилу Рунге (14).
Пример 3. В модельном уравнении (1) выберем следующие параметры: t € [0,50],х(0) = х(0) = 0,Л = 1,ш0 = Ь = 6 = 1,ш = 2^) = 0.8соз(0.5"), ЛН2-р + Н2^2 = 0.9796, р = тах^^к)).
Рис. 5. Фазовая траектория а) и осциллограмма б) для задачи Коши (1), полученные по явной (7) и неявной (8) схемам. [Figure 5. Phase trajectory a) and oscillogram b) for the Cauchy problem (1), obtained by explicit (7) and implicit (8) schemes.]
Погрешность aj Вычислительная точность
2D 4D ВО 160 320 t4L' N - Fffi -EFDS
Рис. 6. а) погрешность e; б) вычислительная точность а. [Figure 6. a) e error; b) computational accuracy а.]
Осциллятор Дуффинга обладает различными колебательными регулярными и хаотическими режимами. Регулярные режимы могут быть много периодическими. На рис. 5 приведен пример двух периодического режима, когда существуют
колебания с несколькими периодами. Осциллограмма на рис. 5б показывает, что со временем колебания выходят на установившийся двупериодический режим, а фазовая траектория рис. 6а имеет форму двух замкнутых петель, что характеризует несколько периодов колебаний.
Из рис. 6а и 6б мы видим, что увеличение узлов расчетной сетки в 2 раза приводит к сокращению ошибки в 2 раза, при этом вычислительная точность метода стремиться к единице.
Пример 4. Рассмотрим случай когда для ЯКРС условия теорем 1 и 2 нарушаются. В модельном уравнении (1) выберем следующие параметры: 1 € [0,50], х(0) = х(0) = 0,Л = 5,^0 = 8,Ь = 6 = = = 0.8соз(0.51:),
ЛН2-р + Н2^0 = 1.15,р = тах(я(1|с))
Рис. 7. а) погрешность е; б) вычислительная точность а. Условия теорем 1 и 2 нарушаются.
[Figure 7. a) е error; b) computational accuracy а. The conditions of Theorems 1 and 2 are violated.]
При нарушении условий теорем 1 и 2 для ЯКРС погрешность (рис. 7а) и вычислительная точность (рис. 7б) метода имеет ярко выраженный немонотонный характер в изменении своих значений, что подтверждает нарушение устойчивости и сходимости схемы (7). Устойчивость НКРС (8), в свою очередь не зависит от условий теорем 1 и 2.
Осциллятор Дуффинга. Хаотические режимы
При исследовании нелинейных систем одной из важных задач является определение типа колебаний - периодического, квазипериодического, случайного, хаотического [7]. Особенностью хаотических колебаний является их высокая чувствительность к малым изменениям начальных условий. Поэтому одним из наиболее надежных способов детектирования хаоса является определение скорости разбе-гания траекторий, которая оценивается с помощью спектра показателей Ляпунова. Спектр максимальных показателей Ляпунова строился по модифицированному алгоритму Вольфа-Беннетина с учетом процедуры ортогонализации Грамма-Шмидта, который был подробно рассмотрен в работе [16], а также с приминением неявной схемы.
Замечание 4. Наличие в спектре хотя бы одного положительного ляпуновско-го показателя означает наличие хаотического режима (асимптотической неустойчивости) рассматриваемой фазовой траектории [16]. Отрицательный показатель свидетельствует о регулярном режиме.
Замечание 5. Хаотические режимы можно определять с помощью сечений Пуанкаре. Если сечения Пуанкаре представляют собой облако, то наблюдается хаотический режим [7].
Пример 5. В задаче (1) выберем следующие параметры: q(t) = 0.8сов2(0.5г), г е [0,50] ,х(0) = х(0) = 0,^0 = 6 = Ь = 1, ш = 2
Рис. 8. Фазовые траектории при а) Л = 0.18; в) Л = 0.6; б) спектр показателей Ляпунова от Л. Точками обозначены сечения Пуанкаре. [Figure 8. Phase trajectories for a) Л = 0.18; c) Л = 0.6; b) the spectrum of Lyapunov exponents of Л. The dots denote the Poincare sections.]
Пример 6. Выберем следующие параметры: h = 0.025,
t e [0,50],x(0) = X(0) = 0,Л = = b = 1, w = 2,6 = 1.
Рис. 9. Фазовые траектории при а) q = 0.15; в) q = 0.6; б) Спектр показателей Ляпунова от q. Точками обозначены сечения Пуанкаре. [Figure 9. Phase trajectories for a) q = 0.15; c) q = 0.6; b) the spectrum of Lyapunov exponents of Л. The dots denote the Poincare sections.]
На рис. 8б и 9б изображены спектры максимальных показателей Ляпунова в зависимости от Л и q , соответственно. При положительных значениях ляпуновских показателей фазовые траектории выходят на хаотический режим (рис. 8а и 9а), а при отрицательных - на регулярный (рис. 8в и 9в).
Осциллятор Дуффинга. Вынужденные колебания
Большой интерес представляет исследование систем при воздействии на них различного рода переменных возмущающих нагрузок. Колебания в таких системах, вызываемые периодическими внешними силами называются вынужденными [18]. Различным по характеру возмущающим силам отвечают различные сценарии поведения колебательной системы, которые, таким образом, уже не определяются полностью собственными характеристиками системы, а отражают ее реакцию на возмущающую силу.
Для практики наибольшее значение, имеют те случаи, когда f(x,t) содержит периодическую функцию зависящую только от времени t, которая имеет свои собственные частоту и амплитуду, отличные от рассматриваемой системы. В таком случае возникает суперпозиция колебаний внешней периодической силы и собственных колебаний системы . Колебания, описываемые уравнением (1) выйдут на некоторый установившийся режим с течением времени. В таких системах, часто можно наблюдать такие явления резонанс и бистабильность [4], [19]. Для изучения этих явлений одной из важных задач является построение амплитудно-частотных (АЧХ), фазо-частотных характеристик (ФЧХ) и добротности.
Определение 4. АЧХ называется зависимость амплитуды установившихся колебаний выходного сигнала некоторой системы от частоты ее входного гармонического сигнала.
Определение 5. ФЧХ называется зависимость разности фаз между выходным и входным сигналами от частоты сигнала.
Определение 6. Добротностью называется количественная характеристика резонансных свойств колебательных систем, показывающая во сколько раз полная энергия системы больше затрачиваемой
В работах [18] и [19] с помощью метода гармонического баланса была доказана следующая теорема.
Теорема 5. Задача Коши (1) эквивалентна линейной задаче Коши с целочисленной производной
x (t)+ p(^,t)x (t)+ s2(^,t)x (t) = 6cos(^t) ,x(0) = Xo,x(0) = y0, В уравнении (16) коэфициенты p(^,t) и s2(^,t) ищутся в виде:
p(w,t) = -2Awq(t)-1 sin^
q(t)n 2
+
(16)
2А^«-2 dt
(InM-¥(1 -q(t)))cos ( ^^ I +
2
2Kt)= w2-2W(t)cos (
nsm
3A2b 4 '
(17)
(18)
s
dt
(ln(w)-¥(1 - q(t))) sin
q(t)n
2
neos
q(t)n
¥(1 - q(t)) = —y + Z--
1
1
=1 yn n +1 — q(t)
- дигамма-функция, где
у = 1 (^ГО — 1) - постоянная Эйлера.
Для линейной задачи Коши известны формулы расчета АЧХ, ФЧХ и добротности [19]
ЛМ =
б
v/(s2(^,t)- w2)2 + p2(w,t)w2' p(w,t)w
ф(^) = arctan
Q =
(s2(w,t)- w2) J ' s(w,t)
(19)
(20) (21)
Построим поверхности АЧХ, ФЧХ и добротности для немонотонной функции q(t) = 0.8соз(0.5^, используя формулы (19), (20), (21).
Пример 7. В модельном уравнении (1) выберем следующие парамет-
ры: q(t) = 0.8соз(0.5^Д е [0,100],Н = 0.05,х(0) = х(0) = О,Л = 1,ш0 = Ь = 6 = 1, ш е [0,3], Нш = 0.6.
Рис. 10. а) АЧХ и b) ФЧХ для оператора (2) с порядком q(t) = 0.8cos(0.5t). [Figure 10. a) AFC and b) PFC for the operator (2) with order q(t) = 0.8cos(0.5t).]
На рис. 10 даны поверхности АЧХ (рис. 10г) и ФЧХ (рис. 10Ь) для немонотонной функции q(t). На рис.11а дана поверхность добротности с учетом изменения показателя дробной производной по закону q(t) = 0.8соз(0.5^. На рисунке 11б изображена поверхность добротности, когда параметр q е [0,1] является независимой переменной. По рис.11б видно, что при уменьшении параметра q добротность увеличивается. Максимальной амплитуде соответствует максимум добротности, а при уменьшении частоты уменьшается добротность. Но сильнее всего добротность зависит от параметра q.
2
Рис. 11. Добротность а) q(t) = 0.8cos(0.5t), б) для q е [0,1]. [Figure 11. Q factor a) q(t) = 0.8cos(0.5t), b) for q е [0,1].]
Построим АЧХ на плоскости. Для этого мы проводим расчет по схемам (7) и (8) с достаточно большим временем моделирования, при котором вынужденные колебания выходят на установившийся режим. Далее фиксируется значения амплитуды при различных значениях частоты внешнего воздействия, получаем зависимость А(ш), которая строится по точкам. Аналитические АЧХ будем строить при фиксированном значении времени, соответствующему максимальному значению амплитуды установившихся колебаний по формулам (19), (20), (21), т.е. для постоянного значения дробной производной.
Пример 8. В модельном уравнении (1) выберем следующие параметры: 1 € [0,100],Н = 0.05,х(0) = х(0) = 0,Л = 1,ш0 = Ь = 6 = 1,ш € [0,3],Нш = 0.6.
Е IFDSAFC-EFDS AFC-AFC |-IFDS AFC-EFDS AFC-AFcl |-IDFS AFC-EFDS AFC-AFc]
Рис. 12. АЧХ дробного осциллятора Дуффинга при различных видах функции
q(t).
[Figure 12. Oscillograms for q(t) = 0.8cos(0.5t) for a) ш = 1.2, b) ш = 1.5, c) ш = 1.9.]
Результаты приведенные на рис. 12(а, б, в) подтверждают то, что НКРС (8) дает более точные результаты, чем ЯКРС (7). При приближении частоты внешней силы к резонансной = 1.2, разница между амплитудой колебаний, полученной по ЯКРС и амплитудой, полученной по НКРС увеличивается.
Рис. 13. Осциллограммы для q(t) = 0.8cos(0.5t) при а) w = 1.2, б) w = 1.5, в) w = 1.9.
[Figure 13. Oscillograms for q(t) = 0.8cos(0.5t) for a) w = 1.2, b) w = 1.5, c) w = 1.9.]
Это хорошо видно на рис. 13(а, б, в). Такое поведение связано с бистабильностью осциллятора Дуффинга.
А в
О-'" 1 I С
-— I ■ г ---г ■ . - Г---r7^1 W
0.5 1 1.5 2 2.5 3
Рис. 14. АЧХ. Петля гистерезиса.
[Figure 14. AFC. Hysteresis loop.]
Рис. 12 и 14 наглядно демонстрируют бистабильное поведение осциллятора Дуффинга с дробной производной переменного порядка. При стремлении частоты внешней силы к резонансной wr (участок AB) амплитуда колебаний начинает возрастать. Достигая порогового значения (точка B) колебания выходят на неустойчивый режим (участок BD), в результате чего происходит скачок амплитуды колебаний с одного устойчивого режима на другой (участок ВС), а затем происходит уменьшение амплитуды. При уменьшении частоты сначала происходит увеличение амплитуды (участок CD), затем при достижении частоты меньшей резонансной (точка C) w < wr происходит скачок с одного режима на другой (участок AD). Далее происходит уменьшение амплитуды. Многоугольник ABCD называется петлей гистерезиса.
Расчеты в статье, проводились с помощью комплекса программ VOFDDE 1.0, разработанного в среде Maple.
Заключение
В статье была рассмотрена неявная конечно-разностная схема (8) для уравнения Дуффинга с производной дробного переменного порядка типа Римана-Лиувилля. Были обоснованы вопросы устойчивости и сходимости неявной конечно-разностной схемы. Для обоснования теоретических результатов приводились тестовые примеры. С помощью правила Рунге сравнивались результаты работы неявной схемы (8) с результатами явной схемы (7). Построены фазовые траектории и осциллограммы для осциллятора Дуффинга с дробной производной переменного порядка типа Римана-Лиувилля, с помощью спектра максимальных показателей Ляпунова и сечений Пуанкаре исследовались хаотические и регулярные режимы. Были построены поверхности добротности, амплитудно-частотной и фазо-частотной характеристик для исследования вынужденных колебаний. Результаты исследования показали, что неявная конечно-разностная схема показывает более точные результаты, чем явная.
Конкурирующие интересы. Авторы заявляют об отсутствии конфликта интересов в отношении авторства и публикации.
Авторский вклад и ответственность. Все авторы внесли свой вклад в эту статью. Авторы несут полную ответственность за предоставление окончательной версии статьи в печать. Окончательный вариант рукописи был одобрен всеми авторами.
Список литературы
1. Нахушев А.М. Дробное исчисление и его приминение. М.: Физматлит, 2003.272 с.
2. Petras I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation. New York: Springer, 2010.180 pp.
3. Кузнецов С. П. Динамический хаос. М.: Физматлит, 2001. 296 с.
4. Зельдович Б. Я., Табирян Н. В. Оптическая бистабильность на ориентационной нелинейности жидких кристаллов, Квантовая электроника, 1984. Т. 11, №12, С. 2419-2426 DOI: 10.1070/QE1984v014n12ABEH006232.
5. Еськов В. В. и др. Хаотическая динамика миограмм, Вестник новых медицинских технологий. Электронное издание, 2016. №3 DOI: 12737/21668.
6. Ejikeme C.L., et al. Solution to nonlinear Duffing oscillator with fractional derivatives using ho-motopy analysis method (HAM), Global Journal of Pure and Applied Mathematics, 2018. vol. 14, no. 10, pp. 1363-1388 ISSN 0973-1768.
7. Syta A. Chaotic vibrations of the Duffing system with fractional damping, Chaos: An Interdisciplinary Journal of Nonlinear Science, 2014. vol.24, no. 1, pp. 10-16 DOI: 10.1063/1.4861942.
8. Xing W. Threshold for chaos of a duffing oscillator with fractional-order derivative, Shock Vib., 2019. vol. 2019, pp. 1-16.
9. Shen Y., Li H., Yang S., Peng M., Han Y. Primary and subharmonic simultaneous resonance of fractional-order Duffing oscillator, Nonlinear Dyn.,2020. vol. 102, pp. 1485-1497.
10. El-Dib Y. O. Stability approach of a fractional-delayed Duffing oscillator, Discontinuity Nonlinearity Complex, 2020. vol. 9, pp. 367-376.
11. Eze S. C. Analysis of fractional Duffing oscillator, Rev. Mex. Fisica, 2020. vol. 66, pp. 187-191.
12. Gouari Y., Dahmani Z., Jebril I. Application of fractional calculus on a new differential problem of Duffing type, Adv. Math. Sci. J., 2020. vol.9, pp. 10989-11002.
13. Syam M.I.The Modified Fractional Power Series Method for Solving Fractional Undamped Duffing Equation with Cubic Nonlinearity, Nonlinear Dyn. Syst. Theory, 2020. vol. 20, pp. 568-577.
14. Barba-Franco J. J., Gallegos A., Jaimes-Reátegui R., Pisarchik A.N. Dynamics of a ring of three fractional-order Duffing oscillators, Chaos Solitons Fractals, 2022. vol. 155, pp. 111-747.
15. Ким В. А. Осциллятор Дуффинга с внешним гармоническим воздействием и производной переменного дробного порядка Римана-Лиувилля, характеризующая вязкое трение, Вестник КРА-УНЦ. Физ.-мат. науки, 2016. Т. 13, №2, С. 46-49 DOI: 10.18454/2079-6641-2016-13-2-50-54.
16. Ким В. А., Паровик Р. И. Расчет максимальных показателей Ляпунова для колебательной системы Дуффинга со степенной памятью, Вестник КРАУНЦ. Физ.-мат. науки, 2018. Т. 23, №3, С. 98-105 DOI: 10.18454/2079-6641-2018-23-3-98-105.
17. Kim V. A., Parovik R. I. Application of the Explicit Euler Method for Numerical Analysis of a Nonlinear Fractional Oscillation Equation, Fractional and fractals, 2022. vol. 6, no. 5, pp. 274-293 DOI: 10.3390/fractalfract6050274.
18. Ким В. А., Паровик Р. И. Исследование вынужденных колебаний осциллятора Дуффинга с производной переменного дробного порядка Римана-Лиувилля, Известия Кабардино-Балкарского научного центра РАН, 2020. Т. 93, №1, С. 46-56 DOI: 10.35330/1991-6639-2020-193-46-56.
19. Kim V. A., Parovik R. I. Mathematical model of fractional Duffing oscillator with variable memory, Mathematics, 2020. vol.8, no. 11, pp. 20-34 DOI: 10.3390/math8112063.
Ким Валентин Александрович А - младший научный сотрудник интегративной научно-исследовательской лаборатории природных катастроф Камчатки - землетрясений и извержений вулканов, Камчатского государственного университета имени Витуса Беринга Петропавловск-Камчатский, Россия, СЖСГО 0000-00018895-6821.
Паровик Роман Иванович А - доктор физико-математических наук, доцент, профессор кафедры информатики и математики, Камчатского государственного университета имени Витуса Беринга, Петропавловск-Камчатский; ведущий научный сотрудник лаборатории моделирования физических процессов, Институт космофи-зических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия СЖСГО 0000-0002-1576-1860.
Vestnik KRAUNC. Fiz.-Mat. Nauki. 2022. vol. 40. no. 3. pp. 179-198. ISSN 2079-6641
MSC 26A33, 34C15 Research Article
Implicit finite-difference scheme for a Duffing oscillator with a
derivative of variable fractional order of the Riemann-Liouville
type
V. A. Kim1, R.I. Parovik1'2
1 Vitus Bering Kamchatka State University, 683032, Petropavlovsk-Kamchatskiy,
Pogranichnaya str., 4, Russia
2 Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS,
7, Mirnaya st., Kamchatka Krai, Yelizovsky district, 684034, c. Paratunka, Russia E-mail: valentinekim@mail.ru, romanparovik@gmail.com
The article considers an implicit finite-difference scheme for the Duffing equation with a derivative of a fractional variable order of the Riemann-Liouville type. The issues of stability and convergence of an implicit finite-difference scheme are considered. Test examples are given to substantiate the theoretical results. Using the Runge rule, the results of the implicit scheme are compared with the results of the explicit scheme. Phase trajectories and oscillograms for a Duffing oscillator with a fractional derivative of variable order of the Riemann-Liouville type are constructed, chaotic modes are detected using the spectrum of maximum Lyapunov exponents and Poincare sections. Q-factor surfaces, amplitude-frequency and phase-frequency characteristics are constructed for the study of forced oscillations. The results of the study showed that the implicit finite-difference scheme shows more accurate results than the explicit one.
Key words: Duffing oscillator, Runge rule, Riemann-Liouville operator, Grunwald-Letnikov operator, amplitude-frequency response, phase-frequency response, Q-factor, Lyapunov exponents, Poincare sections, oscillogram.
d DOI: 10.26117/2079-6641-2022-40-3-179-198
Original article submitted: 24.11.2022 Revision submitted: 05.12.2022
For citation. KimV. A., Parovik R. I. Implicit finite-difference scheme for a Duffing oscillator with a derivative of variable fractional order of the Riemann-Liouville type. Vestnik KRA UNC. Fiz.-mat. nauki. 2022,40: 3,179-198. d DOI: 10.26117/2079-6641-2022-40-3-179-198
Competing interests. The authors declare that there are no conflicts of interest regarding authorship and publication.
Contribution and Responsibility. All authors contributed to this article. Authors are solely responsible for providing the final version of the article in print. The final version of the manuscript was approved by all authors.
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)
© Kim V. A., Parovik R. I., 2022
Funding. Financial support was provided within the framework of the grant of the President of the Russian Federation, No. MD-758.2022.1.1
References
[1 [2
Nakhushev A. M. Drobnoye ischisleniye i yego primineniye [Fractional calculus and its application], Moscow. Fizmatlit, 2003, 272 (In Russian).
Petras I. Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation. New York. Springer, 2010, 180.
Kuznetsov S. P. Dinamicheskiy khaos [Dynamic chaos]. Moscow, Fizmatlit, 2001, 296 (In Russian).
Zel'dovich B.Ya., N.V. Tabiryan, Optical bistability associated with orientational nonlinearity of liquid crystals, Quantum Electron., 1984, 14:12, 1599-1604, DOI:10.1070/QE1984v014n12ABEH006232.
Eskov V. M., et al. Chaotic dynamics of the myograms, Journal of new medical technologies, eEdition, 2016, 3 DOI: 12737/21668.
Ejikeme C. L., et al. Solution to nonlinear Duffing oscillator with fractional derivatives using homotopy analysis method (HAM), Global Journal of Pure and Applied Mathematics, 2018, 14:10, 1363-1388, ISSN 0973-1768.
Syta A. Chaotic vibrations of the Duffing system with fractional damping, Chaos: An Interdisciplinary Journal of Nonlinear Science, 2014, 24:1, 10-16, DOI:10.1063/1.4861942. Xing W. Threshold for chaos of a duffing oscillator with fractional-order derivative, Shock Vib., 2019, 2019, 1-16, DOI:10.1155/2019/1230194.
Shen Y., Li H., Yang S., Peng M., Han Y. Primary and subharmonic simultaneous resonance of fractional-order Duffing oscillator, Nonlinear Dyn., 2020, 102, 1485-1497, DOI:10.1007/s11071-020-06048-w.
El-Dib Y. O. Stability approach of a fractional-delayed Duffing oscillator, Discontinuity
Nonlinearity Complex, 2020, 9, 367-376, DOI:10.5890/DNC.2020.09.003.
Eze S. C. Analysis of fractional Duffing oscillator, Rev. Mex. Fisica, 2020, 66, 187-191, DOI:
10.31349/revmexfis.66.187.
Gouari Y., Dahmani Z., Jebril I. Application of fractional calculus on a new differential problem of Duffing type, Adv. Math. Sci. J., 2020, 9, 10989-11002, DOI: 10.37418/amsj.9.12.82.
Syam M.I. The Modified Fractional Power Series Method for Solving Fractional Undamped Duffing Equation with Cubic Nonlinearity, Nonlinear Dyn. Syst. Theory, 2020, 20(5), 568-577.
Barba-Franco J. J., Gallegos A., Jaimes-Reategui R., Pisarchik A. N. Dynamics of a ring of three fractional-order Duffing oscillators, Chaos, Solitons & Fractals, 2022, 155, 111-747 DOI: 10.1016/j.chaos.2021.111747.
Kim V. A. Duffing oscillator with external harmonic action and variable fractional Riemann-Liouville derivative characterizing viscous friction, Bulletin KRASEC. Physical and Mathematical Sciences. 2016, 13:2, 46-49. DOI: 10.18454/2313-0156-2016-13-2-46-49. Kim V. А., Parovik R.I. Calculation the maximum Lyapunov exponent for the oscillatory system of Duffing with a degree memory, Vestnik KRAUNC. Fiz.-mat. nauki. 2018, 23: 3, 98-105. DOI: 10.18454/2079-6641-2018-23-3-98-105.
Kim V. A., Parovik R. I. Application of the Explicit Euler Method for Numerical Analysis of a Nonlinear Fractional Oscillation Equation, Fractal and Fractional, 2022, 6:5, 274-293, DOI: 10.3390/fractalfract6050274.
Ким В. А., Паровик Р. И. Исследование вынужденных колебаний осциллятора Дуффинга с производной переменного дробного порядка Римана-Лиувилля, Известия Кабардино-Балкарского научного центра РАН, 2020, 93:1, 46-56, DOI: 10.35330/19916639-2020-1-93-46-56.
ISSN 2079-6641
KHM B. A., napoBHK P. H.
[19] Kim V. A., Parovik R.I. Mathematical model of fractional Duffing oscillator with variable memory, Mathematics, 2020, 8:11, 20-34, DOI: 10.3390/math8112063.
A
if*
Rü
IP
ü»
/ IM »
Kim Valentine AleksandrovichA - Junior researcher of the integrative Scientific Research Laboratory of Natural Disasters of Kamchatka - earthquakes and volcanic eruptions, Petropavlovsk-Kamchatskiy, Russia, ORCID 0000-0001-8895-6821.
Parovik Roman IvanovichA - D. Sci. (Phys. & Math.), Associate Professor, Prof., Dep. of Informatics and Mathematics, Vitus Bering Kamchatka State University, Petropavlovsk-Kamchatsky; Leading Researcher, Laboratory for Simulation of Physical Processes, Institute for Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, ORCID 0000-0002-1576-1860.
