CC BY
27
УДК 517.3
Численное решение обыкновенного дифференциального уравнения в случае дробной разности матричного оператора *
Ильин И.А.1, 2, Нощенко Д.С.2, Пережогин А.С.1, 2
1 Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, Камчатский край, c. Паратунка, ул. Мирная, 7
2 Камчатский государственный университет имени Витуса Беринга, 683032, г. Петропавловск-Камчатский, ул. Пограничная, 4
E-mail: drew72156@yandex.ru
Рассмотрено численное моделирование решения обыкновенного дифференциального уравнения в случае обобщения разностного оператора. Численные расчеты приведены для уравнения первого порядка. Установлена связь дробной степени матричного оператора СЛАУ с операторами Грюнвальда - Летникова.
Ключевые слова: численное моделирование, дробный разностный оператор, формула Грюнвальда - Летникова
(с) Ильин И.А., Нощенко Д.С., Пережогин А.С., 2012
MSC 34A08
Computation of differential equation with fractional matrix
Ilyin I.A.1, 2, Noshchenko D.S.2, Perezhogin A.S.1, 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, 683032, Petropavlovsk-Kamchatsky, Pogranichnaya st., 4, Russia
E-mail: drew72156@yandex.ru
Fractional power computation of matrix is represented. Numerical calculations for ordinary differential equation and its fraction analog is done. Results are compared with Grunvald -Letnikov operator.
Key words: computation, fractional differential calculus, Grunvald - Letnikov operator
(c) Ilyin I.A., Noshchenko D.S., Perezhogin A.S., 2012
*Работа выполнена при поддержке Министерства образования и науки РФ в рамках программы стратегического развития ФГБОУ ВПО «Камчатский государственный университет имени Витуса Беринга» на 2012-2016 гг.
Введение
Математическое моделирование физических процессов, протекающих в реальных системах с фрактальными свойствами, выполняется с помощью аппарата дробного дифференцирования [4]. В таких задачах применяются операторы с дробным порядком а, которые в случае целого порядка оператора совпадают с классическими дифференциальными операторами.
В настоящей работе предложен подход к построению модели с дробным оператором через замену классической производной на разностное соотношение и дальнейшее вычисление дробной степени матрицы в системе линейных алгебраических уравнений. Численное моделирование выполнено для обыкновенного дифференциального уравнения первого порядка с начальным условием.
Постановка задачи
Рассмотрим задачу Коши для обыкновенного дифференциального уравнения:
О? + у(г ) = 0, у(0) = 1. (1)
Для этого уравнения запишем разностную схему:
Іу Уі+1 - Уі . „ 1
— «----------, і = 0..п, у0 = 1.
аг т
Таким образом, численную схему для уравнения (1) можно записать в матричном виде:
(і о о о .о Уо (і о о о. .о Уо о
-і і о о .о Уі і о о о. .о Уі о
і о -і і о .о У2 о і о о. .о У2 о
т о о -і і . о Уз + о о і о. . о Уз о
.о. о о о .. .і. W .о. о о о. .. .о. W .о.
(2)
Введем обозначение для матрицы в (2), которая соответствует оператору дифференцирования системы (1):
A=
( і о -і і
о
о
(3)
\о О .. 1/
Вычислим дробную степень а для матрицы А. Положим а = 2. Используем метод Ньютона решения нелинейных уравнений. Вычислительная схема матрицы В = \/Л:
Вп+1 = В+Л'В-1, Во = Л. (4)
2
Используя итерационную схему (4) и систему (2), можно получить численную схему для разностного аналога дробной степени исходного дифференциального оператора уравнения (1).
Реализуем итерационный процесс (4) с начальным приближением Во в виде (3). Алгоритм вычисления запрограммирован в системе компьютерной алгебры Maxima [2]. Приведем пример для матрицы размерности 4x4 (5).
B=
/ 1 0 0 0
-.5 1 0 0
-.125 -.5 1 0
- .0625 -.125 -.5 1
(5)
После этого мы получаем систему (2) с основным оператором представленным матрицей В.
/ 1 0 0 0
1 -.5 1 0 0
-.125 -.5 1 0
- .0625 -.125 -.5 1
У1 +
У2
(т* 0 0 ^ /У0^
0 0 0 У1
0 0 0 У2
0 0 0 W
= 0.
(6)
Решение системы с дробной степенью матрицы А имеет вид, отличный от классического решения системы (2). На рис. 1 приведены решения аналогичных систем с размерностью матриц 100x100, шаг г = 0.05.
Рис. 1. Решение систем линейных уравнений: 1 - классическое решение системы 2-го порядка (2); 2 - линия решение системы порядка 3/2 (6)
Решение системы с матрицей В имеет особенность в начальной точке. При этом асимптотическое поведение решения отличается более длительным режимом затухания в отличии от классического решения.
Сравнение с формулой Грюнвальда - Летникова
Обобщение разностного оператора на случай дробного порядка выполняется с помощью формальной замены целого порядка на дробный. В общем случае, разностный оператор п-порядка с шагом т записывается в виде:
8n f (x) = У Г(п + 1)______________
f () L k!r(n - k + 1)
f (x - kT).
(7)
Полагая порядок оператора п дробным, обозначим через V дробный порядок оператора и получим обобщение [5]:
8v f (x)= У r(v + 1)
f (x) k=0 k!r(v - k + 1)
f (x - kT).
(В)
Корень квадратный из разностного оператора первого порядка (З):
( 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \
-0.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.125 -0.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0625 -0.125 -0.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0
-.039063 -0.0625 -0.125 -0.5 1.0 0.0 0.0 0.0 0.0 0.0
-.027344 -.039063 -0.0625 -0.125 -0.5 1.0 0.0 0.0 0.0 0.0
-.020508 -.027344 -.039063 -0.0625 -0.125 -0.5 1.0 0.0 0.0 0.0
-.016113 -.020508 -.027344 -.039063 -0.0625 -0.125 -0.5 1.0 0.0 0.0
-.013092 -.016113 -.020508 -.027344 -.039063 -0.0625 -0.125 -0.5 1.0 0.0
^—0.01091 -.013092 -.016113 -.020508 -.027344 -.039063 -0.0625 -0.125 -0.5 LO.y
1
Hv f (x) = lim — У v
r(v + 1)
т^0 Tv k=0 k!r(v - k + 1)
f (x - kт)
Проведем вычисления в пакете Maxima для коэффициентов входящих в формулу (В) с показателем дробности v = 2
nu : 1/2;
eq : (-1^k*gamma(nu+1)/gamma(k+1)/gamma(nu-k+1);
makelist( subst( j, k, eq ),j,0,10),eval,float;
В результате вычисления для коэффициентов дробного разностного оператора получаем значения, которые совпадают со значениями матрицы B:
[1, -.5, -.125, -.0625, -.039063]. (9)
Получаем согласование вычисления квадратного корня из матрицы A c формулой (В), которая получена формальной процедурой обобщения разностного оператора на случай дробного порядка.
Возьмем определение оператора Грюнвальда - Летникова Hv и установим связь с вычисленным корнем квадратным из матрицы A:
(10)
Если не вычислять предел при т ^ 0, то Н2 совпадает с выражением корня квадратного из матрицы А. При этом ряд, входящий в формулу Грюнвальда - Летникова, обрывается при задании начального условия. Полученный результат говорит о том, что вычисление значения в текущей точке зависит от всех предыдущих значений.
Матричный оператор степени больше чем единица
Умножим матрицу (3) на (5) в случае размерности 10x10, и в результате матричный оператор будет иметь вид:
AB =
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-1.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.375 -1.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0625 0.375 -1.5 1.0 0.0 0.0 0.0 0.0 0.0 0.0
.02344 0.0625 0.375 -1.5 1.0 0.0 0.0 0.0 0.0 0.0
.01172 .02344 0.0625 0.375 -1.5 1.0 0.0 0.0 0.0 0.0
.006836 .01172 .02344 0.0625 0.375 -1.5 1.0 0.0 0.0 0.0
.004395 .006836 .01172 .02344 0.0625 0.375 -1.5 1.0 0.0 0.0
.003021 .004395 .006836 .01172 .02344 0.0625 0.375 -1.5 1.0 0.0
^.002182 .003021 .004395 .006836 .01172 .02344 0.0625 0.375 -1.5 1.0.
Вычислим коэффициенты входящие в определение оператора Грюнвальда - Лет-никова при а = |. Получим значения коэффициентов:
[1.0, -1.5,0.375,0.0625, .02344, .01172, .006836, .004395, .003021, .002182, .001637].
Коэффициенты в точности совпадают с результатом, который получен в произведении А ■ В.
Рассмотрим задачу Коши для обыкновенного дифференциального уравнения второго порядка:
Л2у
+ y(t )= 0, y(0) = 1, y'(0) = 0.
(11)
В качестве аппроксимации оператора второго порядка выберем вторую разность:
Л 2у Уі+1 - 2у,- + у-1 .
*2 « --------т2-----------------, ' = °"“, (12)
тогда систему линейных уравнений для уравнения (11) можно записать в виде:
1 0 0 0 /У>\ - __L т2 0 0 0 /У>\
-2 1 0 0 У1 + 1 т2 0 0 0 У1
1 -2 1 0 У2 0 1 0 0 У2
0 1 -2 1 У3 V 0 0 1 0 У3
= 0. (13)
В качестве вычислительного эксперимента зададим значение а = | и система
линейных алгебраических уравнений примет вид:
3
т 2
1 0 0 0
-1.5 1 0 0
.375 -1.5 1 0
^.0625 .375 -1.5 1
/W\
У1
У2
+
(-\ 0 0 0\
т 2
1 0 0 0
2т 3/2 0
1 0 0
/V\
У1
У2
= 0.
Построим графики решений для системы второго порядка и - (рис. 2).
1
2
1
Рис. 2. Решение систем линейных уравнений: 1 - решение системы 2-го порядка (13); 2 - решение системы порядка 3/2 (14)
Как и в случае системы для системы с порядком дробности а = - имеется особенность в начальной точке. Однако, асимптотический характер решения совпадает с тем, который хорошо известен для дробного осциллятора [3]. Подобная динамика решения хорошо согласуется с задачами Коши для обыкновенных уравнений с дробной производной Римана - Лиувилля. Чтобы устранить особенность в начальной точке, можно использовать методику предложенную в работе [1]. В целом, представленные результаты могут быть использованы при численном моделировании динамики систем с памятью. В частности, в теории вязко-упругости можно выполнить моделирование релаксационных процессов для сложных природных сред и композитных материалов.
Заключение
Полученные численные решения имеют принципиальное отличие от экспоненциального решения за счет более медленного спадания решения на асимптотике. Характер решения указывает на особенности динамики таких систем, которые проявляются в наличии эффектов памяти. Связь, установленная между разностным оператором дробного порядка Грюнвальда - Летникова и дробной степенью матрицы, указывает на прямую связь двух подходов.
Библиографический список
1. Ильин И.А., Нощенко Д.С., Пережогин А.С. О дробной степени разностного оператора для обыкновенного дифференциального уравнения // Математическое моделирование фрактальных процессов, родственные проблемы анализа и информатики: материалы II Междунар. конф. молод. ученых. Нальчик; Эльбрус, 2012. С. 107-109.
2. Электронный ресурс: URL: maxima.sourceforge.net/
3. Нахушев А.М. Дробное исчисление и его применение. М.: Физматлит, 2003. 272 с.
4. Тарасов В.Е. Модели теоретической физики с интегро-дифференцированием дробного порядка. М., Ижевск: РХД, 2011. 568 с.
5. Учайкин В.В. Метод дробных производных. Ульяновск: Артишок, 2008. 512 с.
Поступила в редакцию / Original article submitted: 03.12.2012
