CC BY
35
УДК 517.925.42 Б01: 10.14529/ттр170109
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ЭРЕДИТАРНОГО ОСЦИЛЛЯТОРА ЭЙРИ С ТРЕНИЕМ
Р. И. Паровик
Работа посвящена вопросам математического моделирования эредитарных колебательных систем с помощью математического аппарата дробного исчисления на примере осциллятора Эйри с трением. Модельное уравнение Эйри было записано в терминах дробных производных Герасимова - Капуто. Далее для этого обобщенного уравнения предложена конечно-разностная схема для численного счета. Рассмотрены вопросы аппроксимации, устойчивости и сходимости такой численной схемы. Приведены результаты моделирования, на основе численного решения построены осциллограммы и фазовые траектории в зависимости от различных значений управляющих параметров.
Ключевые слова: осциллятор Эйри; эредитарность; производная Герасимова -Капуто, конечно-разностная схема; фазовая траектория.
Введение
В настоящее время широко исследуются эредитариые процессы или процессы с «памятью:». В книге В.В. Учайкина [1] эредитарным процессам посвящена целая глава. Понятие эредитарности было введено Вито Вольтерра [2] для установления 3 с1П эздыв&ютце и (эредитарной) временной причинно-следственной связи. Примером такой эредитарной связи является, например, дипольный момент атома, наведенный пролетающей вблизи него заряженной частицей. Так же в работе [2] был рассмотрен эредитарный осциллятор, для которого выведен закон сохранения полной механической энергии. Необходимо отметить, что в математической записи этого закона содержится положительное дополнительное слагаемое, отражающее диссипацию энергии в эредитарном процессе.
С точки зрения математики, понятие эредитарности можно описать интегро-дифференциальными уравнениями, где ядра в подынтегральных выражениях являются функциями памяти. Рассмотрим следующее эредитарное уравнение:
г г
J Кг (г — т) х (т) Лт + xJ К2 (г — т) х (т) + гх (г) = / (г), (1)
0 0
где Кг (г — т), К2 (г — т) - ядра, функции памяти, А - коэффициент трения, / (г) -функция, х (г) - функция смещения, г Е [0,Т] -
время процесса, Т - время моделирования, точки над функцией х (г) означают про-из водные первого и второго порядков соответственно.
Замечание 1. Уравнение (1) описывает эредитарный колебательный процесс с лиг / (г)
уравнение можно считать обобщением известного уравнения Эйри, которое было исследовано английским астрономом Дж. Эйри в 1838 году [3] для описания звезд, как точечного источника света в телескопе. Также колебания Эйри рассматриваются в рамках теории дифракционной оптики для получения лазерных пучков [4].
Замечание 2. Используя аппарат дробного исчисления, который достаточно хорошо разработан [5,6], можно уравнение (1) привести к уравнению с производными дробных порядков, если функции памяти имеют вид:
и~т )1-в и-т) —
К (I — т) = [Г(22-р)-К<( - т) = Г(гЛ)'1 2 0 <7< 1 (2)
Действительно, если в уравнение (1) подставить соотношения (2), то мы приходим к следующему уравнению:
двх <Т) + Лд^ <Т) + гх (г) — ! (г), (3)
йв ( \ — 1 г х (т) ^т Ж ( \ _ 1 г х (т) &т
где догХ (т) — . . J в 1 п догХ (т) — . \ } / производные
г (2 — р) о (г — т) г(1 — 7) о (г — т)
дробных порядков в и 1 в смысле Герасимова - Капуто [7].
Замечание 3. Отметим, что в случае, если в — 2 и 7 — 1, уравнение (3) переходит в классическое уравнение Эйри [3].
Замечание 4. В работе автора [8], в случае, когда Л — 0 и f (г) — 0, было получено точное решение в терминах трехпараметрической функцией Миттаг - Леффлера [9], а также были найдены аналоги функций Эйри первого и второго рода.
Замечание 5. В литературе [10,11] осцилляторы, которые описываются производными дробных порядков, иногда называются фрактальными, так как функции памяти являются степенными функциями вида (2), что характеризует одно из фрактальных свойств среды - степенную зависимость смещения от времени («тяжелые хвосты>). Мы будем называть такие осцилляторы эредитарными, так как фракталь-ность - частный случай эредитарности и функции памяти могут не ограничиваются зависимостями вида (2).
Теория эредитарных колебательных процессов получила широкое развитие в связи с различными приложениями [10-14], что подтверждает актуальность исследований. Поэтому исследуем эредитарный колебательный процесс Эйри более детально.
1. Постановка задачи и методика ее решения
Рассмотрим эредитарное уравнение (3) с начальными условиями:
х (0) — Хо, Х (0) — уо, (4)
где хо и Уо заданные константы.
Замечание 6. Начальные условия (4) согласуются с дробными производными Герасимова - Капуто и являются локальными. Если для описания эредитарных колебательных систем выбрать дробные производные Римана - Лиувилля, то необходимо ставить нелокальные начальные условия как пределы от интегральных операторов [9]. Например, в книге [12] рассмотрены некоторые линейные и нелинейные эредитарные колебательные системы в терминах производной Римана - Лиувилля, однако там были выбраны начальные условия типа (4).
Задача Коши (3) и (4) может быть решена численными методами, например с помощью теории конечно-разностных схем [15]. Разобьем временной отрезок Ь Е [0, Т] на N точек с постоянным шагом т = Т/М и введем сеточную функцию х(Ь^) = х^. Далее с помощью соответствующих аппроксимаций дробных операторов Герасимова-Капуто по формулам [7]:
т-в г ]
двх (т) « г(3_ о^ [(* + !)2-в - к2-в (Х]-к+1 - 2х— + х--),
k=0
j-i
dYtX (т) « к=0 ^ + ^ - k1-^ (xj-к + Xj-k-i)
мы приходим к явной конечно-разностной схеме [16]:
xi = туо + Xo, j = 0, j-i
Xj+i = AjXj - Bxj-i - BYs'Pk (Xj-k+i - 2xj-k + Xj-k-i) -
к=1
j-i
-CY, Qk (Xj-k+i - Xj-k) + Dfj, k=i
A = 2Ai + Bi - jT = Ai = Bi D = 1
j Ai + Bi ' Ai + Bi' Ai + Bi' Ai + Bi Ai = , Bi ^
(5)
Г(3 - ß)' Г (2 - 7)' Pk = (к + 1)2-ß - k2-ß, Qk = (к + 1)i-Y - к-, j = 1, ...,N - 1.
2. Аппроксимация, устойчивость и сходимость схемы
Исследуем вопросы аппроксимации, устойчивости и сходимости схемы (5). Аппроксимация дифференциальной задачи (3),(4) разностной (5) во внутренних узл&х расчетной сетки имеет второй порядок, это следуют из аппроксимации дробных про-ИЗ ВОДНЫХ. Однако за счет аппроксимации в граничных узлах общий порядок аппроксимации схемы является первый. Порядок аппроксимации мо^жно повысить, выбрав специальный вид аппроксимации в граничных узловых точках, например метод фиктивного узла. Сформулируем следующую лемму.
Лемма 1. Коэффициенты, явной конечно-разностной схемы (5) Aj,B,Bi} C,pk,qk обладают следующими свойствами:
Aj > 0, Bi > 0, 0 <qk < 1, 0 <pk < 1, 0 <B < 1, 0 <C < 1. (6)
Рассмотрим структуру явной конечно-разностной схемы (5), которую можно переписать в матричном виде следующим образом:
M • X = F, (7)
10 0 0
-2Ai - Bi Ai + Bi 0 0
M = Aipo т - 2Ai - Bi Ai + Bi 0
Aipi Ai (1 - 2pi) - Bq 2т + Ai (pi - 2) + Bi (qi - 1) Ai + Bi
X = [xo, туо + Xo,... , Xn ]T ,F = [fo,fi,..., fN ]T .
В матричном уравнении (7) треугольная матрица M обратима, ее определитель, с учетом Леммы 1, больше нуля. Поэтому матричное уравнение (7) имеет решение. Далее перепишем схему (5) в виде:
xj+i = (Aj - BB - CC) xj - DDxj-i-
j-i j-i (8) -BY, Pk (xj-k+i - 2xj-k + Xj-k-i) - CYsQk (xj-k+i - Xj-k) + Dfj, k=2 k=2
где BB = piB, CC = qiC ъ DD = B(1 - pi) - Cqi.
Следуя методике работы [17] и с учетом соотношения (8), составим ошибку £j+i = Xj+i - Xj+i между численным xj+i и точным Xj+i решениями. Подставим ошибку £j+i в уравнение (8). После простых преобразований мы приходим к следующему
HGpaBGHCTBy!
|£j+i|<|Aj - BB - CC - DD\\£j\. (9)
Для определения устойчивости и сходимости явной конечно-разностной схемы (5) с учетом неравенства (9) сформулируем теорему.
Теорема 1. Явная схем,а (5) является устойчивой и сходится, если выполняется условие:
0< j Г(3 - ^)Г(2 - 7) тm < 2 . = 01 N 1 ПО)
0 < Г (2 - 7) + АГ (3 - в) т- < 2,j = 0,1,-,N - 1- ^
Доказательство. Условие Теоремы 1 с учетом Леммы 1 должно удовлетворять неравенству [151:
\Aj - BB - CC - DD\ < 1. Откуда мы приходим к условию Теоремы 1. □
Выполнение условия (10) Теоремы 1 зависит от параметров в-, Ъ А т и не зависит от значений индекса j, так как увеличение количества узлов N в расчетной сетке
т
Рассмотрим работу явной конечно-разностной схемы (5) на следующем тестовом примере.
Пример 1. Рассмотрим задачу Коши:
двх (т) + AdYtX (т) + tx (t) = t4 + рЦ-g) + ^ТТЪ), X (0) = X (0) = 0. (11)
Вестник ЮУрГУ. Серия «Математическое моделирование i ал
и программирование» (Вестник ЮУрГУ ММП). 2017. Т. 10, № 1. С. 138-148
Простой проверкой можно убедиться, что точным решением задачи Копти (11) является элементарная функция [7]:
Пусть значения управляющих параметров: N = 100, т = 0,1, Л = 0,02, x0 = y0 =
1,в = 1, 8,7 = 0, 7.
Рис. 1. Конечно-разностная явная схема (кривая 1) и точное решение (кривая 2)
На рис. 1 можно отметить, что численное решение (5) практически совпадает с точным решением (12) задачи Копти (11). Рассмотрим как изменяется абсолютная ошибка и расчетный порядок точности p = ln(| е \) / 1п(т) схемы (5) при различных значениях дробных параметров в и 7-
Из таблицы следует, что абсолютная ошибка между точным решением и численным стремится к нулю при уменьшении шага т. Порядок точности p для схемы (5) можно считать близким к единице. Причем уменьшение значений дробных параметров в и 7 незначительно влияет на точность схемы для этого примера.
Используя теорему Лакса [15], мьт можем сделать вывод, что скорость сходимости численного решения к точному имеет первый порядок. Поэтому можно сделать вывод о том, что схема (5) может быть использована в численном моделировании эредитарного осциллятора Эйри.
В следующих примерах полученную схему (5) мьт будем использовать при построении фазовых траекторий в зависимости от различных значений дробных управляющих параметров.
3. Результаты моделирования
Пример 2. С учетом разностной схемы (5) моделирование проведем со следующими управляющими параметрами: N = 2000, т = 0,05, ф = = 0,3, Л = 0,85,5 = 12 • шв,7 = 0, 75, x (0) = 0, 02, X (0) = 0, f (t) = 5cos (pi).
x (t) = t3.
(12)
и*(0
Таблица
Исследование схемы (5)
в = 1, 8 и 7 = 0, 7
N т Абсолютная ошибка Порядок точности р
10 1/10 0,0558 1,25
20 1/20 0,0339 1,12
40 1/40 0,0195 1,006
80 1/80 0,0108 1,032
100 1/100 0,0089 1,024
320 1/320 0,003 0,999
в = 1, 5 и 7 = 0, 5
10 1/10 0,0899 1,095
20 1/20 0,0528 0,981
40 1/40 0,0289 0,96
80 1/80 0,0152 0,954
100 1/100 0,0123 0,953
320 1/320 0,004 0,955
в = 1, 3 и 7 = 0, 3
10 1/10 0,0989 1,004
20 1/20 0,0562 0,9609
40 1/40 0,0299 0,9508
80 1/80 0,0154 0,9509
100 1/100 0,0124 0,9517
320 1/320 0,0039 0,9577
Визуализация результатов моделирования приведена на рис. 2.
хт
Рис. 2. Осциллограммы эредитарного уравнения Эйри с трением и внешним периодическим воздействием: 1) в = 1, 3; 2) в = 1, 5; 3) в = 1, 7; 4) в = 2
На рис. 2 приведены осциллограммы при различных значениях дробного параметра в и фиксированного значения дробного параметра 7 = 0, 75, которое соответствует фрактальному трению.
в=2
по-видимому, проявляется эффект резонанса (совпадение соответствующих частот 1.
Далее можно отметить, что с некоторого момента времени происходит уменьшение амплитуды и потом она устанавливается.
С каждым уменьшением значений дробного параметра в ^ 1 амплитуды колебаний возрастают, однако с течением времени они устанавливаются примерно на одинаковом уровне. Поэтому можно сделать вывод о том. что расчетным кривым в рассматриваемом случае могут соответствовать траектории на фазовой плоскости, построенные в зависимости от значений дробного параметра в> которые выходят на предельные циклы практически похожие друг на друга (рис. 3).
Фазовые траектории на рис. 3 с течением времени выходят на предельный цикл, причем при в ^ 1 такой выход траекторий на этот режим происходит с каждым разом все медленнее и медленнее.
Пример 3. Рассмотрим другие параметры моделирования: N = 2000, г = 0, 05, ш = 0, 3, ф = 1,5, Л = 2, 5,5 = 12ш^,7 = 0, 9,x (0) = 0, 2,Х(0) = 0,f (t) = 5 cos^t).
На рис. 4 приведены результаты моделирования в виде осциллограмм. На рис. 4 мы можем заметить, что расчетные кривые синфазных колебаний имеют затухающий характер. Однако как в предыдущем случае амплитуда колебаний с течением времени устанавливается. Сохраняется роль дробного параметра
Уменьшение значений параметра в ^ 1 приводит к увеличению амплитуды колебаний, однако со временем амплитуды колебаний можно считать одинаковыми. Такое поведение расчетных кривых указывает на существование одного более общего предельного цикла на фазовой плоскости (рис. 5).
Фазовый портрет на рис. 5 подтверждает, что все траектории выходят на один и тот же предельный цикл. При в ^ 1 происходит более медленный выход траекторий на предельный цикл.
т
Рис. 4. Осциллограммы эредитарного уравнения Эйри с трением и внешним периодическим воздействием: 1) в = 1, 3; 2) в = 1, 5; 3) в = 1, 7; 4) в = 2
т
Рис. 5. Фазовые траектории при различных значения дробного параметра в: 1) в = 2;
2) в = 1, 7; 3) в = 1, 5; 4) в = 1, 3
Заключение
Введение дробных производных Герасимова - Капуто порядков в и 7 в модельное уравнение Эйри расширяют его свойства за счет дополнительных управляющих параметров Это позволяет более гибко проводить исследование эредитарных
колебательных процессов, которое найдет свое отражение в различных приложениях.
Результаты моделирования эредитарного колебательного процесса Эйри показали, что при значениях дробного параметра в ^ 1с учетом внешнего периодического воздействия, что фазовые траектории выходят с замедлением на предельный цикл. Значения дробного параметра 7 определяют степень вязкости трения, что приводит к еще более медленным колебательным процессам.
Определенный интерес, в рамках математического моделирования, представляет учет зависимости дробных параметров в и 1 от времени Ь или от функции решения х Ц).
Литература
1. Uchaikin, V.V. Fractional Derivatives for Physicists and Engineers. Vol. I. Background and Theory / V.V. Uchaikin. - Beijing, Berlin: Higher Education Press, Springer-Verlag, 2013. -373 p.
2. Вольтерра, В. Теория функционалов, интегральных и интегро-дифференциальных уравнений / В. Вольтерра. - М.: Наука, 1982. - 304 с.
3. Airy, G.В. On the Intensity of Light in the Neighbourhood of a Caustic / G.B. Airy // Transactions of the Cambridge Philosophical Society. - 1838. - V. 6. - P. 379-402.
4. Хонина, С.H. Зеркальные лазерные пучки Эйри / С.Н. Хонина, СТ. Волотовский // Компьютерная оптика. - 2014. - Т. 34, № 2. - С. 203-213.
5. Oldham, К.В., Spanier J. The Fractional Calculus. Theory and Applications of Differentiation and Integration to Arbitrary Order / K.B. Oldham, J. Spanier. - London: Academic Press, 1974. - 240 p.
6. Miller, K.S. An Introduction to the Fractional Calculus and Fractional Differntial Equations / K.S. Miller, B. Ross. - N.-Y.: A Wiley-Interscience Publication, 1993. - 384 p.
7. Паровик, Р.И. Математическое моделирование линейных эредитарных осцилляторов / Р.И. Паровик. - Петропавловск-Камчатский: изд-во Камчатского государственного университета имени Витуса Беринга, 2015. - 178 с.
8. Паровик, Р.И. О решении обобщенного уравнения Эйри / Р.И. Паровик // Доклады Адыгской (Черкесской) международной академии наук. - 2014. - Т. 16, № 3. - С. 64-69.
9. Kilbas, A.A. Theory and Applications of Fractional Differential Equations / A.A. Kilbas, H.M. Srivastava, J.J. Trujillo. - Amsterdam: Elsevier, 2006. - 523 p.
10. Mainardi, F. Fractional Relaxation-Oscillation and Fractional Diffusion-Wave Phenomena / F. Mainardi // Chaos, Solitons к Fractals. - 1996. - V. 7, № 9. - P. 1461-1477.
11. Афанасьев, В.В. Стабилизация фрактального осциллятора инерциальными воздействиями / В.В. Афанасьев, М.П. Данилаев, Ю.Е. Польский // Письма в журнал технической физики. - 2010. - Т. 36, № 7. - С. 1-6.
12. Petras, I. Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation / I. Petras.
- Beijing, Berlin: Higher Education Press, Springer-Verlag, 2011. - 218 p.
13. Tavazoei, M.S. Chaotic Attractors in Incommensurate Fractional Order Systems / M.S. Tavazoei, M. Haeri // Physica D: Nonlinear Phenomena. - 2008. - V. 237, № 20. -P. 2628-2637.
14. Rossikhin, Y.A. Application of Fractional Calculus for Dynamic Problems of Solid Mechanics: Novel Trends and Recent Results / Y.A. Rossikhin, M.V. Shitikova // Applied Mechanics Reviews. - 2010. - V. 63, № 1. - P. 010801.
15. Самарский, А.А. Устойчивость разностных схем / А.А. Самарский, А.В. Гулин. - M.: Наука, 1973. - 415 с.
16. Паровик, Р.И. Численный анализ некоторых осцилляционных уравнений с производной дробного порядка / Р.И. Паровик // Вестник КРАУНЦ. Физико-математические науки.
- 2014. - Т. 9, № 2. - С. 30-35.
17. Xu, Y. A Finite Difference Technique for Solving Variable-Order Fractional Integro-Differential Equations / Y. Xu, V. Suat Ertiirk // Bulletin of the Iranian Mathematical Society. - 2014. - V. 40, № 3. - P. 699-712.
Роман Иванович Паровик, кандидат физико-математических наук, доцент кафедры «Математика и физика>, Камчатский государственный университет имени Виту-са Беринга (г. Петропавловск-Камчатский, Российская Федерация); старший научный сотрудник лаборатории моделирования физических процессов, Институт кос-мофизических исследований и распространения радиоволн ДВО РАН (Камчатский край, с. Паратунка), romanparovik@gmail.com.
Поступила в редакцию 15 мая 2016 г. MSC 37С70 DOI: 10.14529/ттр 170109
MATHEMATICAL MODELLING OF HEREDITARITY AIRY OSCILLATOR WITH FRICTION
R.I. Parovik, Kamchatka State University by Vitus Bering, Petropavlovsk-Kamchatsky; Institute Space Physics Research and Radio Wave Propagation FEB RAS, Kamchatka region, Paratunka, Russian Federation, romanparovik@gmail.com
Work is devoted to mathematical modelling hereditarity oscillatory systems with the help of the mathematical apparatus of fractional calculus on the example of an Airy oscillator with friction. Model Airy equation was written in terms of Gerasimov -Caputo fractional derivatives. Next a finite-difference scheme to this generalized equation for numerical computation was proposed. The problems of approximation, stability and convergence of a numerical scheme are considered. The results of simulations are presented based on numerical solutions waveforms and phase trajectories depending on different values of the control parameters are built.
Keywords: Airy oscillator; hereditarity; Gerasimov - Caputo derivative; finite-difference scheme; the phase trajectory.
References
1. Uchaikin V.V. Fractional Derivatives for Physicists and Engineers. Vol. I. Background and Theory. Beijing, Berlin, Higher Education Press, Springer-Verlag, 2013. 373 p. DOI: 10.1007/978-3-642-33911-0
2. Volterra V. Theory of Functional and of Integral and Integra-Differential Equations. N.Y., DOVER, 1959. 304 p.
3. Airy G.B. On the Intensity of Eight in the Neighbourhood of a Caustic. Transactions of the Cambridge Philosophical Society, 1838, vol. 6, pp. 379-402.
4. Honina S.N., Volotovskij S.G. Mirror Laser Airy Beams. Computer Optics, 2014, vol. 34, no. 2, pp. 203-213. (in Russian)
5. Oldham K.B., Spanier J. The Fractional Calculus. Theory and Applications of Differentiation and Integration to Arbitrary Order. London, Academic Press, 1974. 240 p.
6. Miller K.S., Ross B. An Introduction to the Fractional Calculus and Fractional Differntial Equations. N.Y., Wiley-Interscience Publication, 1993. 384 p.
7. Parovik R.I. [Mathematical Modelling of Linear Oscillators Hereditarity]. Petropavlovsk-Kamchatskiy, Kamchatskiy Gosudarstvennyy Universitet Imeni Vitusa Beringa, 2015. 178 p.
8. Parovik R.I. Cauchy Problem of Generalized Airy Equation. Doklady Adygskoy (Cherkesskoy) mezhdunarodnoy akademii nauk [Reports Adyghe (Circassian) International Academy of Sciences], 2014, vol. 16, no 3, pp. 64-69. (in Russian)
9. Kilbas A.A., Srivastava H.M., Trujillo. J.J. Theory and Applications of Fractional Differential Equations. N.Y., Elsevier Science, 2006. 523 p. DOI: 10.1016/0960-0779(95)00125-5
10. Mainardi F. Fractional Relaxation-Oscillation and Fractional Diffusion-Wave Phenomena. Chaos, Solitons & Fractals, 1996, vol. 7, no. 9, pp. 1461-1477.
11. Afanas'ev V.V., Danilaev M.P., Pol'skij Yu.E. Stabilizatsiya fraktal'nogo ostsillyatora inertsiaPnymi vozdeystviyami [Stabilization of Fractal Oscillator Inertial Effects]. Technical Physics Letters, 2010, vol. 36, no. 7, pp. 1-6. (in Russian)
12. Petras I. Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation. Beijing, Berlin, Heidelberg, Higher Education Press, Springer, 2011. 218 p. DOI: 10.1007/978-3-64218101-6
13. Tavazoei M.S., Haeri M. Chaotic Attractors in Incommensurate Fractional Order Systems. Physica D: Nonlinear Phenomena, 2008, vol. 237, no. 20, pp. 2628-2637. DOI: 10.1016/j.physd.2008.03.037
14. Rossikhin Y.A., Shitikova M.V. Application of Fractional Calculus for Dynamic Problems of Solid Mechanics: Novel Trends and Recent Results. Applied Mechanics Reviews, 2010, vol. 63, no. 1, pp. 010801. DOI: 10.1115/1.4000563
15. Samarskij A.A., Gulin A.V. Ustojchivost' Raznostnyh Shem [The Stability of Difference Schemes]. Moscow, Nauka, 1973. 415 p.
16. Parovik R.I. Numerical Analysis Some Oscillation Equations with Fractional Order Derivatives. Bulletin KRASEC. Physical and Mathematical Sciences, 2014, vol. 9, no. 2, pp. 34-38. DOI: 10.18454/2313-0156-2014-9-2-34-38
17. Xu Y., Suat Erturk V. A Finite Difference Technique for Solving Variable-Order Fractional Integro-Differential Equations. Bulletin of the Iranian Mathematical Society, 2014, vol. 40, no. 3, pp. 699-712.
Received May 15, 2016
