Численный анализ задачи Коши для широкого класса фрактальных осцилляторов Текст научной статьи по специальности «Математика»

Вестник КРАУНЦ. Физ.-мат. науки. 2018. № 1(21). C. 93-116. ISSN 2079-6641

DOI: 10.18454/2079-6641-2018-21-1-93-116

УДК 517.938

ЧИСЛЕННЫЙ АНАЛИЗ ЗАДАЧИ КОШИ ДЛЯ ШИРОКОГО КЛАССА ФРАКТАЛЬНЫХ ОСЦИЛЛЯТОРОВ *

Р. И. Паровик1,2

1 Камчатский государственный университет имени Витуса Беринга, 683032 Россия, г. Петропавловск-Камчатский, ул. Пограничная, 4

2 Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034 Россия, Камчатский край, с. Паратунка, ул. Мирная, 7

E-mail: romanparovik@gmail.com, parovik@ikir.ru

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

Ключевые слова: задача Коши, фрактальные осцилляторы, эредитарность, оператор Герасимова-Капуто, численная схема, устойчивость, сходимость, правило Рунге.

© Паровик Р. И., 2018

*Работа выполнена при финансовой поддержке гранта президента РФ № МК-1152.2018.1 и по теме НИР КамГУ имени Витуса Беринга «Применение дробного исчисления в теории колебательных процессов» №АААА-А17-117031050058-9.

Введение

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

Математический аппарат описания наследственной механики — аппарат интегро-дифференциальных уравнений с разностными ядрами, которые называются функциями памяти [11]. В случае, если функции памяти степенные, то мы можем перейти к другим типам уравнений, которые основаны на производных дробных порядков, свойства которых рассмотрены в настольных книгах по дробному исчислению [13]-[15]. Процессы и системы, которые описываются с помощью производных дробных порядков называются фрактальными, например, фрактальными называются осцилляторы в теории эредитарных колебаний [16]. В работе предложена математическая модель широкого класса фрактальных осцилляторов и проведен ее численный анализ.

Постановка задачи

Задача. Рассмотрим следующее эредитарное уравнение: г г

У К (г — п) X (п) Ап + Я (х (г), г) ^К (г — п) X (п) Ап = / (х (г), г), (1)

о о

где К1 (г — п), К2 (г — п) — разностные ядра, которые определяют функции памяти или эредитарности, х (г) — функция смещения, Я (х (г), г) — непрерывная нелинейная функция, которая характеризует трение, / (х (г), г) — непрерывная нелинейная функция, которая характеризует внешнее воздействие, х (г) = Ах(г)/Аг,х (г) = А2х(г)/Аг2 — производная первого и второго порядков.

Интегро-дифференциальное уравнение (1) описывает широкий класс эредитарных осцилляторов, причем, если функция f (х(г),г) линейная, а функция Я (х(г),г) зависит только от г или является постоянной, то мы получаем класс линейных эредитарных осцилляторов, а если функции f (х (г),г) и Я (х (г),г) нелинейные, то мы приходим к классу нелинейных эредитарных осцилляторов. Например, если функция Я (х (г),г) четная, а функция f (х(г),г) нечетная функция по переменной х, то мы приходим к нелинейным фрактальным автогенераторам типа Ван дер Поля [ 17]-[ 19].

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

п1—в п—г

К1 (г — п) = Г (2 — в), К2 (г — п) = щ!—;), 1 < в < 2,0 < у < 1, (2)

Замечание 1. Заметим, что в работе [20] отмечено, что если функцию памяти выбрать в виде К(г — п) = ®8 (г — п), где ® - положительная константа, 8 (•) -функция Дирака, то мы получим отсутствие памяти по аналогии с марковским процессом. Случай, когда в динамической системе функция памяти имеет вид К (г — п) = г — Н (г — п), где Н (•) - функция Хевисайда, приводит к системам с полным наличием памяти. Степенной характер функции памяти в виде (2) приводит к динамическим системам с "частичной потерей памяти".

Замечание 2. Выбор степенной функции в виде (2) обусловлен возможностью перехода к производным дробного порядка и не является ограниченным. Так, если подставить функции (2) в уравнение (1) и принять во внимание определение производной дробного порядка в смысле Герасимова-Капуто [3],[21], то приходим к следующему уравнению:

д0х (п) + Я (х (г), г) д** (п) = / (х (г), г), (3)

которое описывает один из классов эредитарных осцилляторов — фрактальные осцилляторы.

Для уравнения (3) справедливы традиционные начальные условия:

х (0) = «1, х (0) = а2, (4)

где «1, «2 — заданные константы.

Задача (3) и (4) является задачей Коши, которая описывает класс фрактальных осцилляторов с постоянной эредитарностью или памятью. Существование и единственность решения задачи Коши (3) и (4) исследованы в работе [22].

Отметим, что дробные параметры в и у в уравнении (3) могут быть связаны со свойствами среды, например, с фрактальной размерностью и определяются с помощью экспериментальных данных.

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

п 1—в (г) п —у(г)

К1 (г — п) = пТГлл,К2 (г — п)= '

г (2 - в (t))' 24 '7 г (1 - у(0)'

где 1 < в (t) < 2,0 < у (t) < 1.

п 1-в (п) п-у(п)

Ki (t - П )= П ot ^, K (t - п )= '

г (2 - в (п))' ^ '7 г(1 - Y(п))' где 1 < в (п) < 2,0 < у (п) < 1.

п 1-в (^п) п-У^-п)

K1 (t - п )= J -vT, K2 (t - п )= '

г (2 — в (г — п))' ^ 'у г (1 — у (г — п))'

где 1 < в (г — п) < 2,0 < у (г — п) < 1.

Которые при подстановке в уравнение (1) приводят нас к следующим трем задачам Коши:

д0(г }х (п) + Я (х (г), г) д0гг(г }х (п) = / (х (г), г), х (0) = «1, х (0) = «2, (5)

д0вг (п }х (п) + Я (х (г п), г) д0уг(п }х (п) = / (х (г), г), х (0) = «1, х (0) = «2, (6)

д0в{(—п]х (п) + Я (х (г), г) д?-пп]х (п)= f (х (г), г),х (0) = «1,х (0) = «2, (7)

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

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

На рис.1 приведена условная классификация эредитарных осцилляторов, которые мы будет исследовать в нашей статье.

Рис. 1. Классификация эредитарных осцилляторов

Условность классификации на рис.1 обусловлено тем, что существуют другие фрактальные осцилляторы, например, активные фрактальные осцилляторы (АФО), которые не исследуются в настоящей статье.

Активные фрактальные осцилляторы представляют собой радиоэлектронный автогенератор с дробной цепью обратной связи. Математическое описание АФО дается уравнением с классическими целочисленными производными, но с дробным интегралом Римана-Лиувилля в его правой части [25], [26].

Нелокальная явная конечно-разностная схема. Устойчивость и сходимость

Пусть для достижения нужной гладкости при построении явной конечно-разностной схемы функция решения х(г) е С3 (0, Т), функция, отвечающая за трение Я (х (г),г) е С (0, Т), а функция f (х (г), г) удовлетворяет условию Липшица с константой Ь > 0 по переменной х:

|f (хх (г), г) - f (х2 (г), г) | < Ь|хх (г) - х2 (г) |. (8)

Напомним, что условие Липшица (8) означает, что функция f принадлежит классу абсолютно непрерывных функций на интервале (0,Т), т.е. f е АС(0, Т).

Разобьем временной интервал [0, Т] на N равных частей с постоянным шагом т = Т^. Функция решения х (г) перейдет в сеточную функцию х (гк) = хк, где гк = кт,к = 1,...,N — 1. Аппроксимации производных дробного порядка в уравнении (3) согласно работе [11], имеют вид:

к—1

(<§) « А £ о/ (хк—/+1 — 2*к-; + хк—у—0 ,о/ = (7 + 1)2—в — /2—в, (9)

7=0

где А =

Г (3 - в)'

£-1

\!-Г_ ,-1-у

dOtx (<§) « BAk £ bj (x£-j+1 - X£-j-0 ,bj = (j + 1)1-Y- j1-r (10)

j=o

т—^

где 5 = 2Г(2 — у) и функция Як = Я (х ,¿к) .

Подставляя аппроксимации (9) и (10) в исходное уравнение (3) приходим к следующей дискретной задаче Коши:

к—1 к—1 А £ о/ (хк—/+1 — 2хк—/ + хк—/—1) + ВЯк £ Ь/ (хк—/+1 — хк—7—1) = /к, (11)

/=0 /=0

х0 = а1, х1 = а2 + та1 (12)

Для дискретной задаче Коши (11) и (12) справедлива следующая явная конечно-разностная схема:

хк+1 = А+^В (2Ахк — (А — ЯкВ) хк—1 + fk) — (13)

А £-1

- А + hВ £ aj ^j+1 - 2xk-j + Xk-j-^ -

BAk

k- 1

—А+Ъ £ /хк—/+1 — хк—10 .

Явная конечно-разностная схема (13) нуждается в исследовании на устойчивость и сходимость, потому что существует ограничение на шаг т. Заметим, что коэффициенты

схемы (13) обладают свойствами, которые можно оформить в виде следующей вспомогательной леммы.

Лемма 1. Коэффициенты явной конечно-разностной схемы (13) обладают следующими свойствами:

k-1

к-1

1) £ aj = к2-в, £ bj = k1-Y , j=0 j=0

2) 1 = a0 > ai > ■ ■ ■ > 0,1 = b0 > bi > ■ ■ ■ > 0,

3) A > 0, B > 0.

Доказательство. Первое свойство вытекает из определения:

к-1 к-1 г I aj = £ (j + 1) j=o j=o

2-ß - j2-ß

= 1-0 + 22-ß — 1 + 32-ß- 22-ß + ••• +

+ (к- 1)2-в + к2-в - (к- 1)2-в = к2-в Е Ъ) = Е [(7 + 1)1-Г - >1-Г1 = 1 - 0 + 21-Г - 1 + 31-Г - 21-у ••• +

7=0 ]=0

+ (к - 1)1-г+к1-г - (к - 1)1-г = к1-г

Второе условие докажем следующим образом. Введем в рассмотрении функции: ф (г) = (г + 1)2-в - г2-в и п (г) = (г + 1)1-г - г1-г, где г > 0. Эти функции являются монотонно убывающими. Действительно, возьмем производные по х от этих функций. Поучим:

ф (x) = (2 - ß) [(x + 1)1-ß - x1-ß П' (x) = (1 - y) Г(х + 1)1-Y- x1-y

< 0, < 0.

Третье свойство вытекает из свойства гамма-функции. Известно, что гамма функция г (г) монотонно убывающая функция на интервале 0 < г < 1, следовательно функция

монотонно возрастающая функция, причем 0 < т_/ ч < 1. Так как т > 0, то мы

Г (z)

приходим к тому что A > 0, B > 0. □

Г (z)

Исследуем теперь порядок аппроксимации дробных операторов дЦ{х) и d^tx ).

Пусть <d0tx ) и d^x ) - операторы аппроксимации. Тогда справедлива следующая лемма.

Лемма 2. Аппроксимации дЦ{х ) и d0tx (%) операторов Герасимова-Капуто дЦ{х ) and d0tx ) удовлетворяет следующей оценке:

dßtx(S) -dßtx(S) <C1T, ^(S) -doYtx(S)| <C2T,

где С1 и С2 - константы, которые не зависят от шага т.

Доказательство. Используя свойства Леммы 1 и определений [11]:

dßx (S) =

Г (3 - ß)

I aj [x(t - jT) + O (t2)]

j=0

Т2-в £2-в

0 т 2 =

Т 2-в £-1

Т2-в (2-в 2 T(3-5) j=0j(t"jT)+ г(3-в)O(Т ' = Т 2-в £-1 = г(з-в) =j (t- jT) + 0 (т 2' ■

1 k-1 <j+/)T г(2-в) = J П 1-вX(t - П) =

j jT

дов(Х (*)

1 £-1

г (2-в) Ej ((- j' *jg [jT >(j+1) т]

дов(Х (* ) - ^Х (*)

Т2-в £-1

f^-g) £a ^((- jT) -^((- *^ + 0 (Т2

Т 2-в £-1 гС3-7в) j-0 ■0 (Т) + 0 (т 2

Т2-в£2-в

г (3 - в)

0 (т ) + 0 (т 2

0 (т ) + 0 (т 2) = 0 (т ).

Аналогично, мы можем показать вторую оценку в лемме (14).

Т1 Y £-1

доГ(Х (*) = £ bj [x (t- jT) + 0 (Т)]

Г (2 - Y)

j=0

Т1-Y£1-Y

t1-y £-1

f(2-Y) £ j ((-jT) + f(2-Y)

0 (T ) =

.1-7 £-1

t1-Y

Г (2 - Y) j=, j(t-jT) + Г (2 - Y)■

T1-Y £-1

= г(2-Y) j=0 j (t- jT)+0 (T) ■

0 т 2 =

do'x (*) =

1

(j+1)T

Г (1 - Y) |0 J П ^ ('- П) =

JT

£1

Г (1 - Y)j=0 ^x(*) - dOtx(*)|

£bjx (t - *j), *j G [jT, (j + 1) T]

T1-Y £-1

£ bj [x (t - jT) - x (t - *j)] + 0 (T)

T1-Y £-1

гёГ-YÖ £- j0 (T)+0 (T)

т 1-Y£1-Y

г (2 - y)

0 (т ) + 0 (Т )

1

= О (т) + О (т ) = О (т).

Лемма доказана. □

Замечание 4. Можно показать с помощью Леммы 2, что явная конечно-разностная (13) имеем ошибку е = О(т). Мы будет этот факт учитывать при оценке вычислительной точности численного метода.

Формулу (13) удобно представить в матричном виде:

Хк+1 = МХк + (15)

Хк+1 = (XI, Х2,..., Хк )Т, Хк = (хо, ХЬ..., Х-)Т,

/у у \Т

¥к =(У°, А^,..., А^^) ,/о = та2, где матрица (17) М = (ш^), I = 1,...,N — 1,] = 1,...,N — 1:

0, ] > I + 1,

А(2 — а1) — **—!*!, . = . = 3,...,N — 1,

mij =

A + Xi-1B ' (16)

-A (at- j+1 - 2a- j + a- j_i) - BAt-1 (bi- j+1 - hi-j—i)

A + А; 1B

j < i - 1,

2A BXi-ibi-2 -Aa— . 0 1

mi,i = 1, m2,2 = , , a „, mi,i =-——---, i = 2,...,N - 1,

A + A + Ai—1B

A (2ai-2 - a—з) + ¿¿-^¿-3 .

mi,2 =--г---5-, i = 3,...,N - 1.

A + Ai—1B

Для матричного уравнения (15) справедлива следующая теорема.

Теорема 1. Явная конечно-разностная схема (13) сходится с первым порядком \xk — xk\ = O (тX если выполнено условие:

т < То = min

( , , 1 \

1 / 2г (2 - у) \ ß-у

Л На ||г (3 - ß)J

(17)

где ||А || = max\Ak\ - векторная норма. k

Доказательство. Пусть Xk = (x0,...,xN—2) — точное решение системы (18) и ветор ошибок ek+1 = Xk+1 — Xk+1, eo = 0. Тогда система (16), с учетом Леммы 1, можем записать следующим образом:

ek+1 = Mek + Fe,k + O (т), (18)

1 т

Fe,k = A + Ак B (0, \ f (X1, tk—1) - f (X1, tk— 1) \,..., \ f (xk-1, tk-1) - f (Xk-1, tk— 1) \)

< A +АкB (Loeo,L1e1,...,Lkek) = aFtek,

1 T AFk = a + XkB diag(Lo, L1,..., Lk) .

Заметим, что для любого к, выполняется условие |Ьк| < Ь (8). Далее введем в

( 1,

||М|| = max

KKN-1

A1 ВЬо-Aao + |||

A+A^

ВЙ1A2- Aa1 +

A+A2B

Bb2Aз—Aa2 || + ||

A+A^

+

/ £-1 \

max £ mij

' v j=1 )

2A

A+A^

A(2a1 -ao)+BA2bo

A+A^

| A(2a2- al)+BAзbl

A+A3B

А(а2-2а1+ао)-BA3(bo-¿2) A+A3B

+

A(2-g1)-BA2b1 A+A2B

+ +

A(2-q1)-BA3b1 A+A3B

BAi-1 bi-2 -Ag;-2

+

V

+

A(2ai-2-ai-3)+Ai-1Bbi-3

A+A;_ 1B

A+Ai-1B

A(ai-2-2ai-3+ai-4)-Ai-1B(bi-4-bi-2) A+A;_ 1B

A(2-g1)-BAi-1b1

A+Ai-1B

+ ■■■

+

+ ■■

(19)

/

Заметим, что норма матрицы ||М|| < 3, это условие можно получить непосредственно с помощью несложных вычислений (19) и с учетом Леммы 1, а также, что все диагональные элементы матрицы М положительны, если выполнено условие А >

В. Например, для второй строки в норме (19):

2A

A - A1B

2A

3A - A1B

А + Я1В А + Я1В А + Я1В Заметим, что согласно неравенству:

=3

A + A1B 4A1B

A + A1B

3

> 0. Пусть А > Я1В, тогда

4||Я ||В А +

B

< 3.

A1B - A

+

2A

<

А + Я1В А + Я1В мы приходит к другому неравенству:

Я1В — А

A1B - A

A + A1B

+

2A

A + A1B

1

A + A1B

+

2A

A + A1B

Для третьей строки из нормы (19) с учетом А > Я2В:

Ао1 — ВЬ1Я2 + А (1 — 2о1) — ВЯ2 + А (2 — о1) — ВЯ2Ь1

A + BA2

A + A2B

A (3 - 2fl1) - BA2 (2b1 + 1) 3A - В

A + A2B С другой стороны:

A + ||A ||B

A + A2B

(2b1 + 1) < 3_ В

(4 + 2b1)

A+

В

<3.

Bb1A2 - Afl^ A (2fl1 - 1) + BA2 + A (2 - m) - BA2b

<

A + BA2

Bb1A2 - Aa1 A + BA2

+

A + A2B

A (2a1 - 1) + BA2 A + A2 В

<

+

A + A2B

A (2 - a1) - BA2b1 A + A2B

или

2BX2

3--— <

A + X2B <

BbiXi — Aai

A + BX2

+

A (2ai — 1) + BX2

A + X2 B

+

A (2 — ai) — BXibi

A + X2B

при X2 ^ ^ или B ^ те, будем иметь:

i <

BbiXi — Aai

A + BbiX2

+

A (2ai — i) + BX2

A + X2 B

+

A (2 — ai) — BXi bi

A + X2B

Аналогично, суммы остальных строк не превосходят 3 и больше единицы. Из уравнения (18) для любой константы С > 0, независимой от шага т и для ошибки получаем оценку:

lk+i|| < ИД** + M|| \\ekII + CT < ( 3 + j+XBB) Ы + СТ'

L

(20)

Введем для соотношения (18) следующие обозначения: 5к = 3 + -—г—, 5 = Ст. Тогда

А + Х^В

мы получим, следующую оценку:

1к+1 II < 5к\\ек II + 5 < 5к (^-1 ||бк-1 II + 5) + 5 =

= ОДк-1II ек -1II + 5 + 1) < ВДк-1 (5к-211 ек-21 + 5) + 5 + 1) = ВДк-1Лк-2Цек-2 II + ^ (5к5к-1 + 5к + 1) <

< SkSk-1 Sk-2 (sk-3 II ek-31| + s) + s (sksk-1 + Sk + 1) = (21)

= SkSk-iSk-2Sk-3||ek-3L <

< SkSk-1 Sk-2 (Sk-3 II ek-31| + s)+ s (SkSk-1 + Sk + 1) =

= SkSk-1 Sk-2Sk-31| ek-3 II + +S (SkSk-1 Sk-2 + SkSk-1 + Sk + 1) < ... < SkSk-1 ■ ... ■ Sk-r IIek-rII + +S(SkSk-1 ■... ■ Sk-r+1 + ... + Sk + 1). Подставим в оценку (21) r = k- 1, мы получим:

||ek+1 IL < SkSk-1 ■... ■ S111e 11| + s(SkSk-1 ■... ■ S2 + ... + Sk + 1) < Co ||eo|| + O(т).

k

Из второго начального условия (12) из которого следует: ||e1| < ||e0|| и C0 = п sp.

p=1

В силу нашего предположения A > ||А ||B, что равносильно оценке на шаг т:

1

т <( 2г (2 - У) ) в-^. (22)

-||11 г (3 - в) J . (22)

Условие (22) начинает работать в случае, когда число обусловленности матрицы М: д (М) ^ 1, что возникает при достаточно больших значениях Х. В случае, когда число обусловленности д (М) мало, то достаточно выбрать шаг т < 1. Поэтому объединяя

эти два случая мы приходим к оценке на шаг т (17). Теорема доказана.

□

Замечание 5. Заметим, что в случае, когда коэффициентом трения Я постоянный, то условие (17) переходит в условие [27]:

I г

т < ^0 = -о ^ (^ )в — Г).

Замечание 6. Заметим, что в случае в = 2 и у = 1 соотношение (17) переходит

в известное соотношение для классического осциллятора с трением Я и внешним воздействием f [28, С. 205]:

2

т < Т0 = шт ^1,

Замечание 7. В работе [29] при доказательстве сходимости явной конечно-разностной

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

Рассмотрим теперь вопросы устойчивости явной конечно разностной схемы (13). Пусть Гк и Хк два различных решения матричного уравнения (18) с начальными условиями Т и Х0. Тогда справедлива теорема.

Теорема 2. Явная конечно-разностная схема (13) условно устойчива, если выполнено условие (17) и справедлива оценка |Гк — Хк| < С|Г0 — Х0| для любых к, где C > 0 константа, независящая от шага т.

Доказательство. Введем следующее обозначение: ек+1 = Гк+1 — Хк+1. Тогда уравнение (18) может быть записано в форме: е+1 = Мек + . Следуя результатам Теоремы 1, используя следующую оценку:

||М+аЩ|< (3 + А+У = 'к,

мы приходим к оценке:

||ек+1|| < ||М + а^|| ||ек|| < (з + Н^кН =

= 'к||ек| < 'к'к—1 ||ек—11| < 'к'к—1'к—2|ек—2| < ...

< 'к'к—1 ■... ■ 'к—г ||ек—г||.

к

Обозначим г = к — 1, мы получим 11ек+11 < С0 ||е1| < С0 ||е0||, С0 = п 'Р. Последнее

р=1

соотношение следует из второго начального условия (4).

Поэтому, если в вектор начальных условий Х0 добавить возмущение, то это не приведет к значительному увеличению ошибки численного решения. Однако для больших значений Як, возникает большое число обусловленности д ^ 1 и необходимо уменьшать шаг дискретизации т, согласно (24), для малых значений Як достаточно выбрать т < 1. Тогда система (18) устойчива, если выполнено условие (17). Теорема доказана. □

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

Рассмотрим работу явной конечно-разностной схемы (13) на конкретных простейших примерах. Мы с помощью метода двойного пересчета (правила Рунге), покажем, что вычислительная точность метода при уменьшении шага дискретизации т стремиться к единице, что соответствует результатам Леммы 2. Напомним, что правило Рунге для оценки ошибки £ в общем случае можно записать так:

max |x(ti) - x|

£ = -Ч-;-, (23)

ъ 2 p -1 ' v ;

где max |x(ti) — хг-| - максимальная по абсолютному значению ошибка между точным i

x(ti) и численным решением xi, а p - порядок аппроксимации численного метода. Формула (23) справедлива в случае, когда известно точное решение x(ti) задачи Коши. В случае, когда точное решение задачи Коши не известно, например, в силу ее нелинейности, то можно использовать модификацию правила Рунге, метод двойного пересчета:

max |xi — X2i|

£ = —-, (24)

ъ 2 p — 1 v ;

где X2i — численное решение, полученное на шаге т/2.

В нашем случае порядок p в формуле (24), согласно Лемме 2 равен единице. Чтобы получить оценку для вычислительного порядка точности а можно использовать следующую формулу [11]:

а = НЙ • (25)

l0g2(т)

Рассмотрим следующий тестовый пример.

Пример 1. Рассмотрим задачу Коши (3) и (4) со следующими значениями управляющих параметров: непрерывные функции f (x, t) и Я (x, t) выберем f (x (t), t) = —x(t), Я (x (t), t) = 2Aocos (t2), Яо = 0.1,N = 1000,T = 20,т = 0.02,в = 1.8,у = 0.8,x(0) = 0.2,x(0) = 0.1. Мы приходим к следующей задаче Коши:

¿0А(П) + 0.2cos (x2(t)) d00t8x(n)+ x(t) = 0,t g [0,20], (26)

х (0) = 0.2, х (0) = 0.1, (27)

Уравнение (26) описывает нелинейные свободные колебания с памятью. Для задачи Коши (26) и (27) в общем виде конечно-разностная схема (13) может быть записана в виде:

(2А - 1)хк - (А - 0.2ес)(х2) В) х*-А + 0.2ес8 (х2) В

xk+1 = -„ , n ^ ,, _ ( ^ „-- (28)

A k—1

A + 0.2 cos (x2) B j=1

£ ay (xk—j+1 — 2xk—j + xk—j—1) —

0.2Bcos (x2) ( )

Lbj(xk—j+1 — xk—j—0 .

A + 0.2 cos (x2) B j=1

Рис. 2. Осциллограмма для Примера 1

С помощью схемы (28), реализованной в компьютерной программе на языке Maple, построим расчетные кривые — осциллограммы (рис. 2).

Из рис. 2 видно, что колебания имеют затухающий характер. Это обусловлено, тем что колебательный процесс происходит нелинейным трением и с памятью. Поэтому в соотношении для полной механической энергии появляется дополнительное слагаемое, ответственное за диссипацию энергии. Также заметим, что в уравнении (26) отсутствует внешнее воздействие на систему, которое бы подкачивала энергию в колебательную систему.

Исследуем явную конечно-разностную схему (28). Заметим, что для Примера 1 выполнены условия Теоремы 2, т.е. A > ||А||B, а норма матрицы ||M|| = 2.992336728 в матричном уравнении (15), поэтому шаг т можно выбрать меньше единицы.

Выберем в Примере 1 Ао = 4000, остальные значения параметров оставим без изменения и построим согласно расчетной схеме (28) осциллограмму (рис. 3).

0.2015 -0.2010 0.2005 02000

2 4 6 В 10 12 14 16 IS

Рис. 3. Осциллограмма для Примера 1

Здесь мы видим отсутствие сходимости и устойчивости численного решения. Действительно, так как Ао ^ 1 и мы приходим к условию А < ||А ||В и норма матрицы ||М|| = 1.195924608 из уравнения (15). Поэтому необходимо уменьшить шаг дискретизации т согласно условию Теоремы 2, соотношение (17), т.е. т < То = 0.001, что приводит к большому числу вычислений.

Построим фазовую траекторию для Примера 1 по точкам (х (г),х (г)) (рис. 4).

Видно, что фазовая траектория имеет форму закручивающейся спирали в силу затухания колебаний, точка покоя является устойчивым фокусом. Другой подход к построению фазовой траектории был предложен в работах [11],[ 16] (рис. 5).

Рис. 5. Фазовая траектория для Примера 1

Фазовая траектория на рис. 5 построена по обобщенным координатам !х (п),х (г) Такой подход может рассматриваться, как обобщение классического понятия фазовой траектории. Сравнивая рис. 4 и рис. 5 мы можем, сделать вывод о том, что введение обобщенных координат приводит к повороту по часовой стрелке фазовых траекторий при сохранении их формы. По-видимому угол поворота зависит от параметра в.

Проведем оценку вычислительной точности а правилу Рунге (25). Результаты приведены в табл. 1.

Таблица 1

Вычислительная точность а

ß = 1.8, y = 0.8

N Ъ а

10 0.00592 1.712

20 0.00317 1.55

40 0.00166 1.45

80 0.00086 1.39

160 0.00044 1.33

320 0.00022 1.29

640 0.00011 1.26

1280 0.000057 1.24

2560 0.000421 0.91

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

Проверим устойчивость явной конечно-разностной схемы по начальным данным для Примера 1. Абсолютную по модулю ошибку £ мы можем определить по формуле:

£ = тах |х - ^ |, (29)

I

где хе — численное решение, полученное по формуле (28) с учетом добавления малой величины е в первое начальное условие (4). Порядок а найдем по формуле (25). Результаты моделирования приведены в табл. 2.

Таблица 2

Устойчивость по начальным данным

ß = 1.8, у = 0.8

N Ъ а

10 10-5 4.99

20 10-5 3.84

40 10-5 3.12

80 10-5 2.62

160 10-5 2.26

320 10-5 1.99

640 10-5 1.78

1280 10-5 1.6

2560 10-5 1.44

Из табл. 2 мы видим, что максимальная абсолютная ошибка £ в точности совпадает с возмущением е первого начального условия (4), а порядок а стремится к единице. Поэтому условия Теоремы 2 выполняются и явная конечно-разностная схема (28) устойчива.

Рассмотрим другой пример.

Пример 2. Рассмотрим задачу Коши (3) и (4) со следующими значениями управляющих параметров: непрерывные функции f (x, t) и X (x, t) выберем f (x (t), t) = 8 sin (фt) — x(t), X (x(t),t) = 2Xocos (t2), Xo = 0.1,N = 1000,T = 20,т = 0.02,в = 1.8, Y = 0.8,8 = 0.1, ф = 25, x (0) = 0.2, x (0) = 0.1.

В этом примере на колебательную систему действует внешнее гармоническое воздействие с амплитудой 8 и частотой ф. Явная конечно-разностная схема для этого примера имеет вид:

xk+1 —

(2A - 1) xk - (A - 0.2cos (x2) B) jc¿-i _

A

A + 0.2cos (x2) B

k-1

(30)

A + 0.2 cos (x2) B j—1

£ о,- (xk- j+1 - 2xk-j + xk- j-1) -

0.2Bcos (x|

k1

£ (xk-j+1 -xk-j-0 +

5 sin (ф t)

A + 0.2cos (x2) B'

А + 0.2 ^ (х2) В^=1

Так же как и в предыдущем примере условия Теоремы 2 выполняются А > ||А ||В, норма основной матрицы в матричном уравнении (15) ||М|| = 2.992336728 и поэтому схема (30) сходится. Построим осциллограмму для этого примера.

Рис. 6. Осциллограмма для Примера 2

Заметим, что осциллограмма на рис. 6 похожа на осциллограмму, приведенную на рис. 2 для Примера 1. В силу того, что амплитуда внешнего гармонического воздействия мала 8 = 0.1, то мы имеем затухающие колебания.

Выберем в Примере 2. А = 4000, что приводит к нарушению условия в Теоремы 2, т.е. А < ||А||В, ||М|| = 1.195924569 (рис. 7).

Рис. 7. Осциллограмма для Примера 2 при А0 = 4000

В этом случае о сходимости и тем более об устойчивости явной конечно-разностной схемы (30) говорить не приходиться. Построим фазовую траекторию (рис. 8).

Рис. 8. Фазовая траектория для Примера 2. при А0 = 0.1

Видим, что фазовая траектория имеет форму закручивающейся спирали, однако за счет внешнего гармонического воздействия даже при малой амплитуде колебаний 8 фазовая траектория сама совершает колебания с частотой ф. Аналогичную фазовую траекторию с поворотом на некоторый угол по часовой стрелке мы можем получить, построив ее по обобщенным координатам (рис. 9).

Рис. 9. Фазовая траектория в обобщенных координатах для Примера 2 при А0 = 0.1

Рис. 10. Осциллограмма для Примера 2 при А0 = 0.1 и 8 = 50

Рассмотрим Пример 2 в случае, когда амплитуда внешнего воздействия 8 = 50. В этом случае осциллограмма имеет вид, приведенный на рис. 10.

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

Рис. 11. Фазовая траектория для Примера 2 при Я0 = 0.1 и 8 = 50

Рис. 12. Фазовая траектория в обобщенных координатах для Примера 2 при А0 = 0.1 и 8 = 50

Заключение

В работе с помощью численного анализа был исследован класс фрактальных осцилляторов, который характеризуется задачей Коши (3), (4). Была построена нелокальная конечно-разностная схема, изучены вопросы ее устойчивости и сходимости, рассмотрены примеры применения схемы и произведена оценка расчетной точности численного алгоритма. Показано, что явная конечно-разностная схема имеет ограничение на шаг сетки т и первый порядок аппроксимации.

Предложенная схема позволяет найти численное решение задачи Коши (3), (4) и в первом приближении дает хорошие результаты. Для получения более точных результатов необходимо строить другие конечно-разностные схемы более высокото порядка точности или использовать другие численные алгоритмы.

References

[1] Boltzmann L., "Zur theorie der elastischen nachwirkung", Annalen der Physik, 241:11 (1878), 430-432.

[2] Вронский А. П., "Явление последействия в твердом теле", АН СССР. Прикладная математика и механика, 5:1 (1941), 31-56. [Vronskij A. P., "Javlenie posledejstvija v tverdom tele", AN SSSR. Prikladnaja matematika i mehanika, 5:1 (1941), 31-56].

[3] Герасимов А. Н., "Обобщение линейных законов деформации и их приложение к задачам внутреннего трения", АН СССР. Прикладная математика и механика, 12 (1948), 529-539. [Gerasimov A.N., "Obobshhenie linejnyh zakonov deformacii i ih pri-lozhenie k zadacham vnutrennego trenija", AN SSSR. Prikladnaja matematika i mehanika, 12 (1948), 529-539].

[4] Слонимский Г. Л., "О законе деформации высокоэластичных полимерных тел", ДАН СССР, 140 (1961), 343-346. [Slonimskij G. L., "O zakone deformacii vysokojelastichnyh polimernyh tel", DAN SSSR, 140 (1961), 343-346].

[5] Uchaikin V. V., Fractional derivatives for physicists and engineers. V. I: Background and theory, Springer, Berlin, 2013, 373 pp.

[6] Magin R. L., "Fractional calculus models of complex dynamics in biological tissues", Computers & Mathematics with Applications, 59:5 (2010), 1586-1593.

[7] Carvalho A. R. M., Pinto C.M.A., "Non-integer order analysis of the impact of diabetes and resistant strains in a model for TB infection", Communications in Nonlinear Science and Numerical Simulation, 61 (2018), 104-126.

[8] Tarasova V.V., Tarasov V.E., "Economic Interpretation of Fractional Derivatives", Progr. Fract. Differ. Appl, 3:1 (2017), 1-6.

[9] Petras I., Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation, Springer, Beijing and Springer-Verlag Berlin Heidelberg, 2011, 218 pp.

[10] Паровик Р. И., "Дробное исчисление в теории колебательных систем", Современные наукоемкие технологии, 2017, №1, 61-68. [Parovik R. I., "Drobnoe ischislenie v teorii kolebatel'nyh sistem", Sovremennye naukoemkie tehnologii, 2017, №1, 61-68].

[11] Паровик Р. И., Математическое моделирование нелинейных эредитарных осцилляторов, КамГУ им. Витуса Беринга, Петропавловск-Камчатский, 2017, 134 с. [Parovik R. I., Matematicheskoe modelirovanie nelinejnyh jereditarnyh oscilljatorov, KamGU im. Vitusa Beringa, Petropavlovsk-Kamchatskij, 2017, 134 pp.]

[12] Volterra V., "Sur les equations integro-differentielles et leurs applications", Acta Mathe-matica, 35:1 (1912), 295-356.

[13] Oldham K. B., Spanier J., The fractional calculus. Theory and applications of differentiation and integration to arbitrary order, Academic Press, London, 1974, 240 pp.

[14] Нахушев А. М., Дробное исчисление и его применение, Физматлит, М., 2003, 272 с. [Nahushev A.M., Drobnoe ischislenie i ego primenenie, Fizmatlit, M., 2003, 272 pp.]

[15] Kilbas A. A., Srivastava H.M., Trujillo. J. J., Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam, 2006, 523 с.

[16] Meilanov R. P., Yanpolov M.S., "Features of the Phase Trajectory of a Fractal Oscillator", Technical Physics Letters, 28(1) (2002), 30-32.

[17] Паровик Р. И., "Математическое моделирование фрактального осциллятора Ван дер Поля", Доклады Адыгской (Черкесской) Международной академии наук, 17:2 (2015), 57-62. [Parovik R. I., "Matematicheskoe modelirovanie fraktal'nogo oscilljatora Van der Polja", Doklady Adygskoj (Cherkesskoj) Mezhdunarodnoj akademii nauk, 17:2 (2015), 57-62].

[18] Lee J.G., Kim S.W., Bae Y. Ch., "Analysis of Nonlinear Behavior in Fractional Van der Pol Equation with Periodic External Force and Fractional Differential Equation", Journal oftheKIECS, 11:2 (2016), 191-196.

[19] Паровик Р. И., "Об исследовании устойчивости эредитарного осциллятора Ван дер Поля", Фундаментальные исследования, 2016, №3-2, 283-287. [Parovik R. I., "Ob issledovanii ustojchivosti jereditarnogo oscilljatora Van der Polja", Fundamental'nye issle-dovanija, 2016, №3-2, 283-287].

[20] Бутенков С. А., "Математические модели процессов на фрактальных структурах с заданными свойствами на основе методов грануляции", Известия Южного федерального университета. Технические науки, 121:8 (2011), 199-209. [Butenkov S. A., "Matematicheskie modeli processov na fraktal'nyh strukturah s zadannymi svojstvami na osnove metodov granuljacii", Izvestija Juzhnogo federal'nogo universiteta. Tehnicheskie nauki, 121:8 (2011), 199-209].

[21] Caputo M., Elasticita e dissipazione, Zanichelli, Bologna, 1969, 150 pp.

[22] Паровик Р. И., "Существование и единственность задачи Коши для широкого класса эредитарный осцилляторов", Международный научно-исследовательский журнал, 3:10(64) (2017), 112-115. [Parovik R. I., "Sushhestvovanie i edinstvennost' zadachi Koshi dl-ja shirokogo klassa jereditarnyj oscilljatorov", Mezhdunarodnyj nauchno-issledovatel'skij zhurnal, 3:10(64) (2017), 112-115].

[23] Lorenzo C. F., Hartley T. T., Initialization, Conceptualization, and Application in the Generalized Fractional Calculus, Lewis Research Center, NASA, Cleveland, Ohio, USA, 1998, 16 pp.

[24] Ramirez L. E. S., Coimbra C. F. M., "On the selection and meaning of variable order operators for dynamic modeling", International Journal of Differential Equations, 2010 (2010), 846107.

[25] Зайцев В. В., Карлов Ар. В., Нураев Д. Б., "Численный анализ автоколебаний активного фрактального осциллятора", Физика волновых процессов и радиотехнические системы, 16:2 (2013), 45-48. [Zajcev V.V., Karlov Ar.V., Nuraev D. B., "Chislennyj analiz avtokolebanij aktivnogo fraktal'nogo oscilljatora", Fizika volnovyh processov i ra-diotehnicheskie sistemy, 16:2 (2013), 45-48].

[26] Карлов А. В., Регулярные и хаотические колебания в дробных и дискретных осцилляторах, Дис.... канд. физ.-мат. наук: 01.04. 03: Кандидатская диссертация, Самара, 2016, 176 с. [Karlov A. V., Reguljarnye i haoticheskie kolebanija v drobnyh i diskretnyh oscilljatorah, Dis.... kand. fiz.-mat. nauk: 01.04. 03: Kandidatskaja dissertacija, Samara, 2016, 176 pp.]

[27] Parovik R. I., "About one dynamic system, characterizing free oscillations taking into account the variable heredity", International conference on mathematical modelling in applied sciences, 2017, 297-298.

[28] Волков Е.А., Численные методы, Наука, М., 1987, 248 с. [Volkov E.A., Chislennye metody, Nauka, M., 1987, 248 pp.]

[29] Xu Y., Erturk V. S., "A finite difference technique for solving variable-order fractional integro-differential equations", Bulletin of the Iranian Mathematical Society, 40:3 (2014), 699-712.

References (GOST)

[1] Boltzmann L. Zur theorie der elastischen nachwirkung // Annalen der Physik. 1878. Vol. 241, no. 11. P. 430-432.

[2] Вронский А. П. Явление последействия в твердом теле // АН СССР. Прикладная математика и механика. 1941. Т. 5, № 1. С. 31-56.

[3] Герасимов А. Н. Обобщение линейных законов деформации и их приложение к задачам внутреннего трения // АН СССР. Прикладная математика и механика. 1948. Т. 12. С. 529-539.

[4] Слонимский Г. Л. О законе деформации высокоэластичных полимерных тел // ДАН СССР. 1961. Т. 140. С. 343-346.

[5] Uchaikin V. V. Fractional derivatives for physicists and engineers. Volume I. Background and theory. Berlin: Springer, 2013. 373 p.

[6] Magin R. L. Fractional calculus models of complex dynamics in biological tissues // Computers & Mathematics with Applications. 2010. Т. 59. №. 5. pp. 1586-1593.

[7] Carvalho A. R. M., Pinto C. M. A. Non-integer order analysis of the impact of diabetes and resistant strains in a model for TB infection // Communications in Nonlinear Science and Numerical Simulation. 2018. vol. 61. pp. 104-126.

[8] Tarasova V.V., Tarasov V.E. Economic Interpretation of Fractional Derivatives // Progr. Fract. Differ. Appl. 2017. vol.3, No. 1. pp. 1-6.

[9] Petras I. Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation. Beijing and Springer-Verlag Berlin Heidelberg: Springer, 2011. 218 p.

[10] Паровик Р.И. Дробное исчисление в теории колебательных систем // Современные наукоемкие технологии. 2017. № 1. С. 61-68.

[11] Паровик Р.И. Математическое моделирование нелинейных эредитарных осцилляторов. Петропавловск-Камчатский: КамГУ им. Витуса Беринга, 2017. 134 с.

[12] Volterra V. Sur les equations integro-differentielles et leurs applications // Acta Mathemat-ica. 1912. Vol. 35, no. 1. pp. 295-356.

[13] Oldham K. B., Spanier J. The fractional calculus. Theory and applications of differentiation and integration to arbitrary order. London: Academic Press, 1974. 240 p.

[14] Нахушев А. М. Дробное исчисление и его применение. М.: Физматлит, 2003. 272 с.

[15] Kilbas A. A., Srivastava H. M., Trujillo. J. J. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier, 2006. 523 p.

[16] Meilanov R.P., Yanpolov M.S. Features of the Phase Trajectory of a Fractal Oscillator // Technical Physics Letters. 2002. vol. 28(1). pp. 30-32.

[17] Паровик Р.И. Математическое моделирование фрактального осциллятора Ван дер Поля // Доклады Адыгской (Черкесской) Международной академии наук. 2015. Т. 17, № 2. С. 57-62.

[18] Lee J.G., Kim S.W., Bae Y.Ch. Analysis of Nonlinear Behavior in Fractional Van der Pol Equation with Periodic External Force and Fractional Differential Equation // Journal of the KIECS. 2016. Vol. 11, no. 2. P. 191-196.

[19] Паровик Р.И. Об исследовании устойчивости эредитарного осциллятора Ван дер Поля // Фундаментальные исследования. 2016. № 3-2. С. 283-287.

[20] Бутенков С.А. Математические модели процессов на фрактальных структурах с заданными свойствами на основе методов грануляции // Известия Южного федерального университета. Технические науки. 2011. Т. 121, № 8. С. 199-209.

[21] Caputo M. Elasticita e dissipazione. Bologna: Zanichelli, 1969. 150 p.

[22] Паровик Р.И. Существование и единственность задачи Коши для широкого класса эредитарный осцилляторов // Международный научно-исследовательский журнал. 2017. Т. 3, № 10(64). С. 112-115.

[23] Lorenzo C. F., Hartley T. T. Initialization, Conceptualization, and Application in the Generalized Fractional Calculus. Cleveland, Ohio, USA: Lewis Research Center, NASA, 1998. 16 p.

[24] Ramirez L.E.S., Coimbra C.F.M. On the selection and meaning of variable order operators for dynamic modeling // International Journal of Differential Equations. 2010. Vol. 2010. P. 846107.

[25] Зайцев В.В., Карлов Ар. В., Нураев Д.Б. Численный анализ автоколебаний активного фрактального осциллятора // Физика волновых процессов и радиотехнические системы. 2013. Т. 16, № 2. С. 45-48.

[26] Карлов А. В. Регулярные и хаотические колебания в дробных и дискретных осцилляторах: Дис.... канд. физ.-мат. наук: 01.04. 03: Кандидатская диссертация / Самара. 2016. 176 с.

[27] Parovik R. I. About one dynamic system, characterizing free oscillations taking into account the variable heredity // International conference on mathematical modelling in applied sciences. Saint Petersburg: Peter the Great Saint Petersburg Polytechnic University, 2017. P. 297-298.

[28] Волков Е.А. Численные методы. М.: Наука, 1987. 248 с.

[29] Xu Y., Erturk V.S. A finite difference technique for solving variable-order fractional integro-differential equations // Bulletin of the Iranian Mathematical Society. 2014. Vol. 40, no. 3. P. 699-712.

Для цитирования: Паровик Р. И. Численный анализ задачи Коши для широкого класса фрактальных осцилляторов // Вестник КРАУНЦ. Физ.-мат. науки. 2018. № 1(21). C. 93-116. DOI: 10.18454/2079-6641-2018-21-1-93-116

For citation: Parovik R. I. Numerical analysis of the Сauchy problem for a wide class fractal oscillators, Vestnik KRAUNC. Fiz.-mat. nauki. 2018, 21: 1, 93-116. DOI: 10.18454/2079-66412018-21-1-93-116

Поступила в редакцию / Original article submitted: 20.12.2017

Vestnik KRAUNC. Fiz.-Mat. Nauki. 2018. no.1(21). pp. 93-116. ISSN 2079-6641

DOI: 10.18454/2079-6641-2018-21-1-93-116

MSC 34A08, 34K28, 37N30

NUMERICAL ANALYSIS OF THE CAUCHY PROBLEM FOR A WIDE CLASS FRACTAL

OSCILLATORS 1

R. I. Parovik12

1 Vitus Bering Kamchatka State University, 683032, Kamchatskiy kray, 4, Pogranich-naya Str., Russia

Petropavlovsk-Kamchatsky, Russia

2 Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS , 684034, Kamchatskiy kray, 7, Mirnaya Str., Paratunka, Russia

E-mail: romanparovik@gmail.com, parovik@ikir.ru

The Cauchy problem for a wide class of fractal oscillators is considered in the paper and its numerical investigation is carried out using the theory of finite-difference schemes. Fractal oscillators characterize oscillatory processes with power memory or, in general, with heredity, and are described by means of integro-differential equations with difference kernels - memory functions. By choosing memory functions as power functions, integro-differential equations are reduced to equations with derivatives of fractional orders. In this paper, using the approximation of the fractional derivatives of Gerasimov-Kaputo, a non-local explicit finite-difference scheme was developed, its stability and convergence are justified, estimates of the numerical accuracy of computational accuracy are presented. Examples of the work of the proposed explicit-finite scheme are given. It is shown that the order of computational accuracy tends to unity as the number of grid nodes increases and coincides with the order of approximation of the explicit finite difference scheme.

Key words: Cauchy problem, fractal oscillators, hereditary, Gerasimov-Caputo operator, numerical scheme, stability, convergence, Runge rule.

© Parovik R.I., 2018

1 This work was supported by the grant of the President of the Russian Federation No. MK-1152.2018.1 and on the topic of the research of Vitus Bering Kamchatka State University "Application of fractional calculus in the theory of oscillatory processes"No.AAAA-A17-117031050058-9.