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

DOI: 10.18454/2079-6641-2015-11-2-88-95

УДК 517.925.42

конечно-разностные схемы для фрактального осциллятора с переменными дробными порядками

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

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

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

E-mail: [email protected]

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

Ключевые слова: конечно-разностная схема, сходимость, устойчивость

(с) Паровик Р.И., 2015

MSC 37С70

finite-difference scheme for fractal oscillator with a variable fractional

order

R.I. Parovik1, 2

1 Institute of Cosmophysical Researches and Radio Wave Propagation Far-Eastern Branch, Russian Academy of Sciences, 684034, Kamchatskiy Kray, Paratunka, Mirnaya st., 7, Russia

2 Vitus Bering Kamchatka State University, 683031, Petropavlovsk-Kamchatsky, Pogranichnaya st., 4, Russia

E-mail: [email protected]

The paper deals with the explicit finite difference schemes for the fractional oscillator. The questions of approximation, stability and convergence of these schemes.

Key wards: finite-difference scheme, convergence, stability

(c) Parovik R.I., 2015

Введение

Численные методы являются важным инструментом в математическом моделировании различных процессов. Они дают возможность получить приближенное решения любого модельного уравнения, а численный алгоритм можно легко реализовать на ЭВМ. Например, как правило, не всегда удается получить аналитическое решение для линєйного обыкновенного дифференциального уравнения с непостоянными коэффициентами и поэтому невозможно для него получить семейство расчетных кривых.

В этой работе мы рассмотрим один из численных методов - конечно-разностные схемы для решения фрактального или эредитарного осциллятора, рассмотренного в работах автора [1, 2]. Линейные эредитарные осцилляторы можно представить в виде осцилляционного уравнения с операторами дробного дифференцирования Герасимова-Капуто порядков 1 < в < 2 и 0 < у < 1 с коэффициентом трения X и внешним воздействием f (t):

гж d0tx(т) =

d0etx (т) + d0Ytx (т)+ A (t) x (t) = f (t), x (0) = xo,x (0) = yo,

1 t An

(1)

t x(n)(T) dT

j ----n+1, n < a < n + 1 — оператор Герасимова-Капуто;

0tV4'y r (n - a) 0 (x - т)с

A (t) — некоторая известная функция, которая определяет вид линєйного эредитарного осциллятора.

В работах [3, 4] были предложены некоторые конечно-разностные схемы для решения линейных и нелинейных фрактальных осцилляционных уравнений. Можно отмєтить, что построение таких разностных схем было основано на аппроксимации дробной производной Римана-Лиувилля, дискретной производной Грюнвальда-Летникова [5].

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

d0et(t}x (т) + d0rt(t}x (т)+ A (t) x (t) = f (t), x (0) = x0,x (0) = y0,

(2)

Дифференциальные уравнения с производными переменного дробного порядка встречаются в различных приложениях [6, 7, 8, 9].

Для определенности будем считать в уравнении (2) A(t) = юв(t).

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

Фрактальный осциллятор с переменными дробными порядками дается следующим уравнением с начальными условиями:

d0et (t)x (т) + d0Yt(t)x (т) + а)в (t)x (t) = f (t), x (0) = x0, x (0) = y0, (3)

Уравнение (3) характеризует фрактальный осциллятор, частный случай которого рассмотрен в работе [10]. Производные дробного переменного порядка в уравнении (3) можно аппроксимировать следующим образом [11]:

(t)x (П)

т Ц/

j-1 г

г (3 - в/

£ (k + 1)2-в/ - к2-в/

k=0

x/—k+1

2xj k+x

/-k+x/-k-1

(4)

doY(t )x (П)

г—! С1 г

£ (k + 1)1-Yj - k2-YJ

Г (2 - Yj) k=0

Также справедливо другое представление [11]: j_1 т-0k

к=о г (3 - Pk)

xj-k+1 xj-k

d00t(t )x (n) ~ £

(k + 1)2 0j - k2 0j (xj-k+1 - 2xj-k + xj-k-1) , (5)

j-1

d0r x (n) ~ £

k=0

T Yk

г (2 -Tk )L

(k + 1)1 Yj - k1 Yj (xj-k+1 - xj-).

Подставляя в уравнение (3) аппроксимации (4) и (5), получим две явные конечноразностные схемы. Первая схема:

Х1 = Ao xo + fo, j = 0,

j-1

xj+1 = Ajxj - Bjxj-1 - Bj £ Pk (xj-k+1 - 2xj-k + xj-k-1) -

k= 1

j-1 .

-Cj £ qk (xj—k+1- xj-k) + fj+1j

k=1

2A1 + B\ - rnPj A1 B1

A' = -■ -■ , BJ = j , CJ = .i.„i,

A1 + B{

A1+ B{ A1 + B{

aJ т 0j dJ Ят Yj

A, = ----—, B1 =

1 г (3 - Pj)’"1 г (2 - Yj)’

Pi! = (k + 1)2-0j - k2-0J, qk = (k + 1)1-Yj - k1-YJ, j = 1,..., N - 1.

, а вторая схема:

Х1 = A0 X0 + f0, J = 0,

(6)

j-1

xJ+1 = AJxJ - BxJ-1 - C £

-Pk

k=1 г (3 - Pk)

Pk (xj—k+1 - 2xJ-k + xJ-k-1) -

j-1 AT -Yk

-D £ г (2 _ Yfe) qk (xJ-k+1- x-k)+fj+l, 2A0 + B0 - ®ejD a^ 1

A / =---;--------, B = —----, C = D =

A0 + B0

A0 + B0

A0 + B0

, т 00 Ят 10

A0 = ----Б~г, B0 =

(7)

Г (3 - 00)’ 0 Г (2 - %)’

pk = (k + 1)2-0j - k2-0j, qk = (k + 1)1-Yj - k1-Yj, J = 1,...,N - 1.

Заметим, что согласно второму начальному условию из (3): Х1 = ту0 + Х0. Поэтому мы это условие будем учитывать при численном моделировании.

Аппроксимация, устойчивость и сходимость явных конечно-разностных схем

В этом разделе рассмотрим вопросы аппроксимации, устойчивости и сходимости явных конечно-разностных схем (6) и (7). Сначала отметим свойства коэффициентов, входящих в эти схемы.

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

1) 0 < Bj, Cj, B, C, D, pJk, qj < 1;

2) A,Ao,Bo,Ai,Bi,Aj,,aPj,b{ > 0.

Рассмотрим структуру явной конечно-разностной схемы (6), которую можно переписать в матричном виде следующим образом:

M ■ X = F,

(8)

M-

i 0 0 0

К - 2A°i - B°i Ai + Bi 0 0

Ai Po ові - 2Ai - Bi A1 +b1 0

Ai Pi Ai ^ - 2 p2)- BW\ + A2 (pi - 2) + Bi (qi - 1) Ai + B-

X = [ло ТУ0 + x0^. . , Xn]T , F = [fo, /l,. . ., InY .

В силу Леммы диагональные элементы треугольной матрицы M положительны (определитель больше нуля) и, следовательно она обратима. Поэтому матричное уравнение (8) разрешимо и численное решение для уравнения (3) может быть получено. Аналогичные рассуждения справедливы и для явной конечно-разностной схемы

(7).

Далее перепишем схему (6) в виде:

xj+i = (AAj BBj CCj)xj Djxj-i

(9)

j-i . j-1 .

-Bj £ Pk Yj-k+i - 2xj-k + xj-k-i) - Cj £ qJk {xj-k+i - xj-k) + fj,

k=2 k=2

где AAj = Aj- Bjp{- Cjqi, BBj = pi Bj, CCj = qi Cj, Dj = Bj- 2В.р{ - ji.

Следуя методике работы [12], введем ошибку, как разность: є,+1 = xj+1 -xcj+1, где

x-j+i - точное значение решения x (t) в точке t = tj+i, которую подставим в уравнение (9). Ошибка будет удовлетворят следующему неравенству:

|£j+i |<|AAj - BBj - CCj Це, |.

(10)

Следуя методике работы [12], введем ошибку, как разность: ej+i = xj+i -Xj+1, где xj+i - точное значение решения x(t) в точке t = tj+i, которую подставим в уравнение (9). Ошибка будет удовлетворят следующему неравенству:

|£j+i |<|AAj - BBj - CCj ||£j |.

(11)

Сформулируем следующую теорему.

Теорема 1. Явная конечно-разностная схема (6) устойчива, если выполняется условие:

(12)

A, < < A, (3 - pi) + Bj (2 - qj

Доказательство. Перепишем условие (11) в виде:

2A, + B, — — A, pi — Bi'qi

|£j+i| <

A, + Bj

|£j|.

(13)

Очевидно, что в случае Aj — < 0 с учетом Леммы выполняется условие:

2A, + Bj — a>Pj — Aj p, — Bi'qi

Aj + Bj

< і.

С другой стороны условие

2Aj' + Bj — a)Pj — Aj pj — Bi'qi A, + B,

> -1

может быть выполнено при условии:

rnPj < A^3 — p,) + Bj (^2 — qj) .

Поэтому согласно (12) мы имеем:

2Aj' + Bj — aPj — Aj pj — Biqi

< 1

Aj + Bj

Следовательно ошибка ограничена и с уменьшением шага т выполняется условие: £j+1 ^ 0 при j ^ ~.

□

Аналогично можно сформулировать теорему для явной конечно-разностной (7).

Теорема 2. Явная конечно-разностная схема (7) является устойчивой, если выполняется условие:

Ao < < 3Ao + 2Bo — Aj pj — Bjqj', (14)

т—pi Ат—Yi

где Ai = Г(3 — ві),Bi = г(2 — у,) •

Доказательство этой теоремы проводится аналогично с учетом Леммы.

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

Сходимость явных схем (6) и (7) с учетом их устойчивости по теоремам 1 и 2 следует из теоремы Лакса-Рихтмайера [13], согласно которой, существует сходимость первого порядка численного решения к точному. Действительно, можно показать, что выполняется неравенство |x(j — Xj | < Ст, где С - произвольная константа, которая не зависит от шага т.

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

Для решения вопросов аппроксимации, устойчивости и сходимости явных схем (6) и (7), рассмотрим следующую задачу Коши:

дв (t)x (т) + Xd0Yt(t)x (т) + ав (t >x (t)

212-в (t) 2X t2-Y(t)

г (3 - в (t)) + Г (3 - Y(t))

+ ав (t)t2,

(15)

с однородными начальными условиями x

2t2-в (t)

f (t) = Г (3 - в (t)) ■

(0) = X (0) = 0 и правой частью: 2X t2-Y (t)

Г(3 - Y(t))

+ юв (t )t2

Задача (15) имеет точное решение:

x (t) = t2.

(16)

Рассмотрим численное решение задачи (15), полученное по явным конечно-разностным

схемам (6) и (7) с учетом следующих управляющих параметров: в 0 3t

0.8 - , t є (0, T), T = 1000,N = 2000, т = 0.5,X = 0.02, a = 1, tj =

Результаты моделирования приведены на рис.1 и в табл. 1.

(t ) = 18 - 0T5t, Y (t ) = jT, j = 0, ••• , N - 1.

Рис. і. Расчетные конечно-разностные явные схемы: 1 - схема (6); 2 - схема (7)

Таблица 1

Сходимость явных схем (6) и (7)

t Точное решение Ошибка для (6) Ошибка для (7)

0.05 0.0025 0.0425 0.0425

0.1 0.01 0.084 0.084

0.15 0.0225 0.119 0.12

0.2 0.04 0.15 0.152

0.25 0.0625 0.175 0.18

0.3 0.09 0.19 0.2

Из табл. 1 следует, что численные решения, полученные по явным схемам (6) и (7) достаточно хорошо аппроксимируют точное решение (16). Так же можно отметить, что первая конечно-разностная явная схема (6) несколько лучше аппроксимирует точное решение, чем вторая явная схема (7), что показано на рис. 1. Поэтому отмєтим, что расчеты по формулам (6) и (7), можно считать схожими.

Можно так же отметить, что при уменьшении шага дискретизации т абсолютные ошибки для схем (6) и (7) стремятся к нулю. Действительно результаты моделирования приведены в табл. (2).

Таблица 2

Сходимость явных схем (6) и (7)

т Точное решение Ошибка для (6) Ошибка для (7)

0.05 0.0025 0.0425 0.0425

0.025 0.000625 0.0218 0.0218

0.0125 0.00015 0.011 0.011

0.00625 0.00007 0.0074 0.0074

Заметим, что в этом примере аппроксимация имеет первый порядок. Действительно по определению аппроксимация р-го порядка характеризуется неравенством |£j| = |x (tj) — Xj| < Стр, где С - произвольная константа, которая не зависит от шага т. Прологарифмировав это неравенство приходим к соотношению: ln |ej| < lnС+рlnт. Построив график прямой в двойном логарифмическом масштабе, определим порядок аппроксимации р как угол ее наклона (рис. 2).

Рис. 2. Порядок аппроксимации явных конечно-разностных схем: а) - схема (6); б)-схема (7)

Из рис.2 видно, что порядок аппроксимации для первой схемы (6) р = 1.065, а для второй схемы (7) - р = 1.038.

Заключение

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

Результаты моделирования на конкретном примере показали, что расчеты по формулам (6) и (7), можно считать примерно схожими и следовательно их можно использовать в моделировании колебательных систем.

Библиографический список

1. Паровик Р. И. Численный анализ некоторых осцилляционных уравнений с производной дробного порядка // Вєстник КРАУНЦ. Физико-математические науки. 2014. Т. 9, № 2. С. 30-35.

2. Паровик Р. И. Математическое моделирование эредитарного осциллятора // Компьютерные исследования и моделирования. 2015. Т. 7, № 5. С. 1002-1023.

3. Momani S., Odibat Z. Numerical comparison of methods for solving linear differential equations of fractional order // Chaos, Solitons & Fractals. 2007. Vol. 31, no. 5. P. 1248-1255.

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

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. Mainardi F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena // Chaos, Solitons & Fractals. 1996. Vol. 7, no. 9. P. 1461-1477.

7. Coimbra C.F.M. Mechanics with variable-order differential operators // Annalen der Physik. 2003. Vol. 12, no. 11-12. P. 692-703.

8. Sheng H., Sun H.G., Coopmans C. et al. A physical experimental study of variable-order fractional integrator and differentiator // The European Physical Journal Special Topics. 2011. Vol. 193, no. 1. P. 93-104.

9. Shichang Ma, Yufeng Xu, Wei Yue. Numerical solutions of a variable-order fractional financial system // Journal of Applied Mathematics. 2012.

10. Мейланов Р. П. К теории фильтрации в пористых средах с фрактальной структурой // Письма ЖТФ. 1996. Т. 22, № 23. С. 40-43.

11. Sun H. G., Chen W., Wei H., Chen Y. Q. A comparative study of constant-order and variable-order fractional models in characterizing memory property of systems // The European Physical Journal-Special Topics. 2011. Vol. 193, no. 1. P. 185-192.

12. Yufeng Xu, Vedat Suat Erttirk. 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.

13. Самарский A.A. Теория разностных схем. М.: Наука, 1977. 656 с.

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