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

-----------------------------------------------^

ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ

УДК 519.6

В. В. Благовидов, А. И. Лобанов

МНОГОЗНАЧНЫЙ МЕТОД РЕШЕНИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ С ДРОБНЫМИ ПРОИЗВОДНЫМИ

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

A multivalue method for decision of differential equations including fractional derivatives is offered. A research of approximation and stability of the method is carried out.

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

Key words: differential equations, fractal derivity, multivalue method, stability and convergence of the difference scheme.

Введение

Краевые задачи для дифференциальных уравнений с дробными производными возникают при описании процессов переноса в средах с фрактальной размерностью (см. [1—3]). В таких уравнениях порядок дробной производной связан с фрактальной размерностью среды (см. [2]). Способы вычисления фрактальной размерности описаны в работе [4].

Метод для решения дифференциальных уравнений д^ = f (и, £), а є (0,1) (да ° —

регуляризованная дробная производная порядка а) предложен в [5] и применен для решения первой начально-краевой задачи для уравнения диффузии дробного порядка в многомерных областях. Порядок аппроксимации O(hа). Разностные аппроксимации оператора дробного дифференцирования выводятся в работе [6].

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

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

Цель работы — создание многозначного метода для численного решения задачи Коши для уравнения дробного порядка вида

Щх^) = /(и> х) / (!)

и(0) = V , (2)

где ^хи(г) = ГТ—) Г и(' Г(х) — гамма-фУнкния-г ( а) о (х — г)

2. Построение многозначного метода

Для решения задачи Коши (1 — 2) составим трехкомпонентный вектор Нордсика [8, с. 379], используя в качестве компонентов вектора члены разложения функции и в ряд Тейлора:

( ,2 ЛТ

U„ =

где h — шаг по независимой переменной t.

Вестник Российского государственного университета им. И. Канта. 2010. Вып. 10. С. 59 — 68.

Чтобы связать значения вектора Нордсика на и-м слое с его знач1 запишем формулу Тейлора:

У2 У3

ип = ии-1 + Уии-1 + — ип-1 + — еп , У 6 [О, к].

2 6

на предыдущем слое,

(3)

Первый компонент вектора Нордсика выразим через компоненты вектора на предыдущем слое. Для выражения оставшихся компонент вычислим производные интерполяционного полинома при у = к :

2 у3

кип = кип-1 + у ип-1 +^- е„ , у2 •• у2 •• у3

у Мп = у ип-1 +у еп .

(4)

(5)

Систему (3—5) можно записать в матричном виде:

' 1 1 1 ^ ' 1/3 '

ип = п 1 2 ип-1 + 1

V П п 1V V 1 V

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

Для решения задачи (1 — 2) необходимо построить разностный оператор, аппроксимирующий &ахи(1). Для этого разобьем отрезок [0, х] на N отрезков длины к = х / N (в случае переменного шага к надо лишь ввести следующие изменения: х0 = 0, кп = хп - хп-1), тогда хп = ки, п = О, N. Воспользуемся формулой [1, с. 46]

ш)=£ “(Г>;к:<х-ь>‘-а +~4 к* )(х -, г-*.

к=0 Г( к +1 -а) Г(и -а) ь

В рассматриваемом случае (Ь = 0) выразим оператор дробной производной оах и(Ь) двумя способами, приняв и равным нулю в первом случае и единице — во втором:

ип X

ПЛ п

^«(і) =

Г(1 -а) Г(1 -а) П

1-а •П хп ,

| и(і)(х - і) а йі,

■ | и(і)(х - і)1 а йі.

(6)

(6')

Г(1 -а) Г(2 -а) Г(2 -а ) 0

Приближенно вычислим интеграл в выражении (6). Выразим интеграл в правой части через сумму интегралов по элементарным отрезкам:

хп п кк

| и(і)(х - і)-айі = Х | и(і)(хп - і)-а йі.

П к=1кк-к

Заменим и(і) на значения компонентов вектора Нордсика в узлах

кк п кк

к

1-а

X | и(і)(Хп - і)-а йі «X І 2

к=1(к-1)к 2

ч1-а^ п-1

ик-1 + ик ( - і)-а аі =

к=1(к-1)к

(

2(1 -а)

ипп

1-а

1 -

п -1

+ ип + Х ик ( - к +1)1 а-(п - к -1)1 а) к к=1 I

Тогда разностный аналог (6) имеет виц

naxu(t);

Un x

Пл n

h

1-a

Г(1 -а) 2Г(2 -а)

c<|un (n1 “- (n -1)1 a) + Un + n£^Uk ((n - 2 +1)1 а- (n - k -1)1 а).

'в

(7)

Аналогично уравнения для (6') получим

, 1-а 7„2-а f n-1

Unxn + Цпх,

+

h2

Un +ZUk((n-k+1)2 “-(n-k-1)2 “U. (7')

Г(1 -а) Г(2-а) 2Г(3-а) [ к=1

В выражении (7') отсутствует к2-аи0 (и2-а - (и -1)2-а)(2Г(3 -а)) 1, поскольку и0 = 0. Утверждение. Равенства (7) и (7') аппроксимируют оператор дробного дифференцирования с порядком 0(к2) и 0(к) соответственно.

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

кк

J = 1 j I - U(0] (х„ - t)-a dt.

к=1 кк-к

Оценим погрешность на элементарном отрезке. Для этого рассмотрим разложения щ и и(Ь) в ряд Тейлора в окрестности ик-1:

ик=и(хк-1+к)=ик-1+ик-1(хк- хк-1)+щк-1(хк х-1) +0(к3),

(Ь - хк-1 )2 з

и(Ь) = ик-1 + ик-1(Ь - хк-1) + ик-1-7-^ + 0(к ) ,

таккак t - xti , /г. Тогда получим

kh

Jk = Uk-1 j

kh-h

xk xk-1 2 (t Xk-1)

(xn -1) а dt +

kh

-Uk-1 j

kh-h

(xk - xk-1 )2 (t - xk-1)

2

( -1)-а dt + O(h5-а) = J1 + J2 + O(h5-а).

1

Введем новую переменную интегрирования | = t - xk-1 и оценим Jk :

! kh

J1 = Uk-1- j (4)4)(x« -xk-1 -C“ d^ = ~Uk-1

ah3

■Ah-h

4(x« - xk-0

1+a

+O(|3) = O(^).

Для модуля J2 имеем оценку

k-11 t„3-a

(n - a +1)1 a-(n - k)

\1-a

Наконец, оценим модуль J: |м (t )|

h3 a^ (n-2 +1)1 a-(n-2)

\1-a

k=1L

1-a

= O(h2)+max I U(t)|h2 = O(h2).

4 fe[n,xM]

Такая же оценка аппроксимации справедлива и приближенного вычисления интеграла, входящего в формулу (6'). Но, учитывая опущенное слагаемое h2~aUn (n2-a -(n-1)2-a)(2r(3-a)) 1 (из-за того что Un = П),

X

4

к 2-%(И2-а- (п - 1)2-а )п,2-ал 2

2-а >

□ к2-а(п2-а - (п - 1)2-а) =

□ х^-а (2 - а)к = 0(к).

' 1 - ±Л і а - 2

1 -

х

V Лп V

2Г(3 -а)

_ „2-а /„ г„\2-а _ -.2-

= Хп (хп к) = Хп

Таким образом, равенства (7') аппроксимирует оператор дробного дифференцирования с порядком О (к). Утверждение доказано. □

Введем обозначения:

Б1 = X ик (( - к+1)1-а-(и - к - 1)1-а), Б? = X2 ик ((и - к+1)2-а-(п - к - 1)2-а) .

к=1 ’ к=1 ' 7

С учетом систем (4, 5) равенства (7, 7') примут вид

спхи>

Г(1 -а) 2Ц2-а)

(8)

<^ «їх) + х{и ( -(и-1)1-а) +(1+21-х)ип-1 +кип-1 +| еп+БИ |;

- {(1 + 21 “ )ип-1 + кеп + БП } . (8')

Оп“пи(і )>

ип хп-. + ип хпх-а + к1-а

Г(1 -а) Г(2 -а) 2Г(2 -а)

Рассмотрим проекцию функции правой части (1) на разностную сетку /т = /(ит, хи). Так как она зависит от неизвестного значения ии, то для решения уравнения линеаризуем в окрестности

ии-1 /и = /и-1+(/иТи-1и -«и-1)+о((«и -«и-1)2)=£.1 + (/иии -«и-1)+о(к2). (9)

Получаем явное выражение для вычисления еи и уточнения и0 :

(, 1-

Аи

(- (п - 1)1-а) к3

Л

2Г(2 -а)

1-а

х

Г(2 -а)

4Г(2 -а)

і 3-а

к3

2Г(3 -а)

ип 1 = ( Ф,

еп V

(1п)

где с учетом выражения (9)

фп « /п-1 + (/и )«-1(ип - ип-1) -

хп а (1 + 21 а )ип-1 + кип-1 + БП .1-

Vп « /пП-1 + (&)ИП-1 (ип - ип-1) -

Г(1 -а) 2Г(2 -а)

ипхп а (1 + 22 а)ип-1 + ;„2-

Г(1 -а)

2Г(3 -а)

с точностью до погрешности аппроксимации интегральных сумм.

Для того чтобы система (10) имела единственное решение, необходимо и достаточно, чтобы детерминант матрицы Аи был отличен от 0:

к

4-2 а

-{(1 -а)и1 а+ (и -1)1 а} п.

det А„ = --

4Г(2 -а)Г(3 -а)1

Это заведомо верно при а < 1.

Найдем еи, воспользовавшись правилом Крамера:

(( 1-а / _ \1-а \

1

det Ап

-det

(а-(и - 1)1-а)

2и

1-а

ф.

^п V

1-а

2Г(2 -а)

-(апип - апип-1 - кЪпип-1 + 2к2Спип-1 + Яп ) .

<

е.„ =

п

Здесь введены следующие обозначения:

Г(3 -а) [и1- + (и - 1)1-а]

О-п —

ка (/и )И-1, Ьп —

,(і + 21-а)(2-а)

2 [(1 -а)И а+(п -1)1 а]

2 [(п1-а - (и - 1)1-а )(1 + 21-а) - 2(2 - а)п1-а ]

(1 -а)П а + (и-1)1 а

(1 -а)иа а +(и -1)1 а

Оп —

к аГ(3 -а)

(1 -а)пг а +(и -1)1 а

(и1-а+ (и - 1)1-а )/пП-1 -

1-ач ^ (и1 -а+ (и - 1)1-а )ка

Г(1 -а)

и0 +

+ и----(«_1)— 2 -2(2-а)и1-ак 1-а5и

2Г(3-а) п V ’ п

Вместо выражения (9) также можно использовать итерационный метод Ньютона.

Используем еп для вычисления значения вектора Нордсика на и-м слое. Вместо вектора

(1/3, 1, 1)т будем использовать Ь — (?1,12,-3)Т. Значения параметров 11,12,-3 находим из требований устойчивости метода. Так как вектор Нордсика на следующем шаге представим в виде

Г1 1 1 1 Г -11 к2 '

и — 0 1 2 и„-1 + -2 Г апип - апип-1 - кЬпип-1 + ~Спип-1 + 2

V 0 0 1У V -3 у V У

то для линеаризованного уравнения легко имеем явное выражение

Г1 1 - -1 Ьп 1 -

1 - -1 ап 1 -

и„ — 1 - -1 Ьп 0 1 - -2 Ьп +-2 ап 1 - -1 ап 2 + -2 Сп +

1 --А 0 --3Ьп + 7-^-3ап V 1 - -1 ап 1 + -3 Сп +

Л

1 - -1 Сп

1 - 11 ап

1 - -Л

1 -1

ип-1 +

V-3 у

1 - -1 ап

1ип У

3. Устойчивость

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

Г 1 1 --А 1 -

1 - -1 ап 1 -

6и — 0 1 - -1 Ьи 1 - -2 Ь„ +—^ -2 ап 1 - -1 ап 2 + -2 Сп +

0 V 1 --А --3 Ь„ +—^ -3 ап 1 - -1 ап 1 + 3 п +

1 - -1Сп 1 - 11 ап 1 - -1Сп

-2 ап

-3 ап

1 - 11 ап У

и-1 •

(11)

Рассмотрим спектр оператора послойного перехода (11). Собственные числа матрицы есть

— 1 А п — Рп ~^[Рп - где

А1 — 1, ^2 3 — , где

Рп — 2 - Ьи-2 + 7"^ ап12 + Сп-3 + ^

1 а„ - 1 а„ -

г —

' п

г 1 - Ь1 іг

1 - ЬЛ + —Г а„/2

1 - ап11 у

1 + С„?1

1 + Сп-3 +

Л

ап-3

1 - ап-1 У

+ а„/э

1 - апк

2 + Сп -2 + ---------ап-2

1 - ап-1 У

Собственные значения матрицы (11) по модулю не превосходят 1 при -1 — 0,5, -2 — 0,3, -3 — 0,1 для значений ка(/и)пп-1 є[-0,5;0,4] и при -1 —-0,1, -2 —-0,3, -3 — 0,1 для ка(/и)пп-1 є[-0,2;0,5]. Значения коэффициентов находились численно с использованием МаЬЬаЪ И2009а.

4. Тестирование метода

Работоспособность метода проверялась при решении тестовых задач. Так, известно [7, с. 47] что

х '---4 / \а 2а-1 2а-1

(х-аТ , .-г,

О-

соэл/ х - а Г +-----------

= г_£Л_ + 2 2 (х-а) 4 /2а_1 (\lx-a|, геО , Иеа ^0.

2

т"» _ _.-а ,,ч _ _ СОБуїх + 0,1

В качестве теста выбиралась задача иахи(т) — 0,1 +-----------------------------------------------------------. с начальным данными

ф + 0,1

0 11+а 2а-1

и(0) — —-+ 2а-1/2 • 0,1 4 /а-1/2(л/0Д).

Г(1 + а) а 1/2

На рисунках 1—4 представлены результаты численных расчетов.

Х

Рис. 1. Точное (сплошная линия) и приближенное (пунктирная линия) решения тестовой задачи при

а — 0,3, к — 1

Х

Рис. 2. Абсолютная погрешность решения при а — 0,3, к — 1

100 200 300 400

Х

Рис. 3. Точное (сплошная линия) и приближенное (пунктирная линия) решения тестовой задачи при а — 0,9, к — 1

Рис. 4. Абсолютная погрешность решения при а = 0,9, к = 1

5. Численная оценка реально достигнутого порядка аппроксимации

Реальный порядок аппроксимации метода находился численно. Для различных значений а задача считалась при разных значениях шага. Затем вычислялась евклидова норма ошибки егг(к) на последних 10 % узлов сетки, чтобы оценить реально достигнутый: порядок аппроксимации погрешности метода после завершения разгонного участка.

На рисунке 5 приведен достигнутый порядок аппроксимации в зависимости от порядка производной а.

о С. О О о о о о о о о о о С1 . С1 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

ПОРЯДОК ДРОБНОЙ ПРОИЗВОДНОЙ

Х

< 2.5

1.5

Рис. 5. Зависимость порядка аппроксимации метода от порядка дробной производной

------------------------------------------------------------------Ш------------------------

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

Список литературы

1. Чукбар К. В. Стохастический перенос и дробные производные // Журн. эксперим. и теор. физики. 1995. Т. 108, вып. 5, № 11. С. 1875 — 1884.

2. Кобелев В. Л., Кобелев Я. Л., Романов Е. П. Недебаевская релаксация и диффузия в фрактальном пространстве // Докл. РАН. 1998. Т. 361, № 6. С. 755 — 758.

3. Кочубей А. Н. Диффузия дробного порядка // Дифференциальные уравнения. 1990. Т. 26, № 4. С. 660—670.

4. Mandelbrojt B. B. The fractal geometry of nature. San-Francisco, 1983.

5. Таукенова Ф. И., Шхануков-Лафишев М. Х. Разностные методы решения краевых задач для дифференциальных уравнений дробного порядка / / Журн. вычислительной математики и математической физики. 2006. Т. 46, № 10. С. 1871 — 1881.

6. Oldham K. B., Spanier J. The fractional calculus (theory and applications of differentiation and integration to arbitrary order). N.Y.; London, 1974.

7. Самко С. Г., Килбас А. А., Маричев О. И. Интегралы и производные дробного порядка и некоторые их приложения. Минск, 1987.

8. Хайрер Э., Нерсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. М., 1990.

9. Нахушев А. М. Уравнения математической биологии. М., 2003.

10. Нахушев А. М. Дробное исчисление и его применение. М., 2003.

Об авторах

Валерий Валерьевич Благовидов — студ., МФТИ, e-mail: val_blag@rambler.ru

Алексей Иванович Лобанов — д-р физ.-мат. наук, проф., МФТИ, e-mail: Alexey@crec. mipt.ru

Authors

Valeriy Blagovidov — student, MIPT, Moscow, e-mail: val_blag@rambler.ru

Prosessor Aleksey Lobanov — MIPT, Moscow, e-mail: Alexey@crec.mipt.ru