CC BY
180
Вычислительная математика
УДК 517.98
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ КРАЕВОЙ ЗАДАЧИ ДЛЯ ДВУМЕРНОГО УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ С ПРОИЗВОДНЫМИ ДРОБНОГО ПОРЯДКА
В. Д. Бейбалаев1, М. Р. Шабанова2
1 Дагестанский государственный университет,
367025, Махачкала, ул. М. Гаджиева, 43 а.
2 Институт проблем геотермии Дагестанского НЦ РАН,
367030, Махачкала, пр-т Шамиля, 39 а.
E-mails: kaspij_03@mail.ru, meil-sha@yandex.ru
Получено решение краевой задачи для двумерного уравнения теплопроводности с производными дробного порядка по времени и по пространственным переменным методом сеток. Разработаны явные и неявные разностные схемы. Доказаны критерии устойчивости этих разностных схем. Показано, что порядок аппроксимации по времени равен единице, а по пространственным переменным равен двум. Предложен метод решения с помощью дробных шагов. Доказано, что модуль перехода, соответствующий двум полушагам, аппроксимирует модуль перехода для. данного уравнения.
Ключевые слова: численные методы, устойчивость, аппроксимация дробных производных, дифференциальные уравнения дробного порядка.
Отсутствие адекватных количественных моделей процессов тепломассо-переноса в многокомпонентных системах со сложной пространственно-временной структурой, для которых характерны эффекты памяти, пространственные нелокальности и самоорганизации, привело к тенденции пересмотра основных положений теории тепломассопереноса [1]. Значительные успехи при рассмотрении кинетических явлений в системах с фрактальной структурой связаны с использованием формализма интегро-дифференцирования дробного порядка. Повышенный интерес к дифференциальным уравнениям дробного порядка обусловлен их физической интерпретацией, что и вызвало их широкое применение в механике, физике, биофизике, экономике — практически во всех областях естествознания. Как указано в [2], несмотря на долгую историю развития математического аппарата дробного дифференцирования, аналитические методы решения уравнений с производными дробного порядка оказываются малоэффективными, а теория численных методов их решения носит фрагментарный характер и далека от завершения. Таким образом, разработка численных методов решения краевых задач для дифференциальных уравнений с частными производными дробного порядка является актуальной.
Ветлугин Джабраилович Бейбалаев (к.ф.-м.н.), ст. преподаватель, каф. прикладной математики. Мумина Руслановна Шабанова, научный сотрудник, лаб. математического моделирования и мониторинга геотермальных объектов.
Пусть С — прямоугольник {0 < х < а, 0 < у < Ъ} с границей Г. В области И = О х (О, Т] рассмотрим эволюционную (нестационарную) краевую задачу для двумерного уравнения теплопроводности с производными дробного порядка по времени и по пространственным переменным.
Задача. Найти решение и(Ь,х,у) € С2(Д уравнения
В&и = СхГ)$хи + СуГ)$уи + /(1,х,у), (¿,ж,у)€Д (1)
удовлетворяющее начальному и граничным условиям:
и(О, х, у) = ф(х, у), (х, у) € О,
и(г,х,у)\г =<р(х,у), (х, у) € С,
где и = и(1;,х,у) —искомая, /(¿,ж,у) —заданная функции; Сх, Су —коэффициенты теплопроводности; 0 < а ^ 1, 1 < /3 ^ 2; П^и(1,х,у) = 0^хи{Ь,х,у) = (Л)21^, 0^хи(Ь,х,у) = (щ)21оу 13 —частные дробные производные Римана—Лиувилля;
ТП , ч 1 /"* и{з,х,у) ,
1п+и{Ь, х, у) = —— / --------—¿8
0+1 У) Г (а) Л (¿-в)1"“
— частные дробные интегралы Римана—Лиувилля по соответствующей переменной [3].
Задачу (1), (2) будем решать численным методом. Численным методам решения краевых задач для уравнений теплопроводности с производными дробного порядка посвящено много работ. В работе [4] разработаны разностные схемы решения краевой задачи для уравнения теплопроводности с производной дробного порядка по времени. Разработке численных методов решения краевых задач для дифференциальных уравнений с частными производными дробного порядка посвящены также работы [2,5-19].
Пусть /(ж) — некоторая функция с областью определения И. Тогда имеет место конечно-разностное определение летниковской производной порядка /3 € М в точке х € Д/) [20]:
- х~тгк)■
к=О
где 1 < /3 < 2; до = 1, дк = (~1)кр(р — 1) •••(/? — к + 1 )/к\.
Согласно теореме Летникова [20], если / € С2(О), то производная Грюн-вальда—Летникова совпадает с производной Римана—Лиувилля.
Для аппроксимации дробных производных Римана—Лиувилля по переменным х, у при 1 < Р < 2 на отрезках [0, а], [0,6] воспользуемся формулой Грювальда—Летникова со смещением [8, 9]:
[х/Щ
= Д+о ЙР ^ ^ (Х ~ ^ ~ ^
к=О
где к = х/М. Как установлено в [8, 9], формула (3) обеспечивает более точную аппроксимацию, чем стандартная формула Грюнвальда—Летникова.
Воспользовавшись формулой (3), для производных дробного порядка Ри-мана—Лиувилля по пространственным переменным в случае 1 < /3 < 2 получим
п+1
-^Охи(^’Х’У)\Хп ~ Тд '^^1(1]и{^^хп-]+1^у) ^
(4)
т+1
^>0уи(^^х,у)\ут ~ -Гр qju(t,x,ym-j+l),
3=0
Где Xn—jJгl ~ Хп 1)Л-) Уm—j-\-1 ~ 2/т, 0' 1)/?..
Пользуясь достаточным признаком существования дробной производной Римана—Лиувилля [3] при 0 < а < 1, на отрезке получим
, 1 ( и(и,Х,у) и'(х,8)й8 \
о* ’1^ Г(1 - а) V (^+1 - **)“ + 1к (и+1 -8)«)'
где
гоо
Г(а) = / ха~1е~хг1х. Jo
Представляя производную и\т,х,у) на отрезке [^,^+1] в виде конечной разности
/йи\ и(гк+1, ж, у) - и(гк,х, у)
\с1т) к Т ’
разностную аппроксимацию производной дробного порядка а на отрезке [Ьк, ¿^+1] можно записать в виде
рау(+ т 7л| ~ ^ и(^к,х,у)
~и(4+1,ж,у) -и{гк,х,у) Г^+1 (18
т Л* (4+1 - «)с
_ ц(4+1,а;,у) - аи(и,х,у)
Г(1 - а)(1 - а)т“ ’ 1 ;
Для нахождения решения задачи (1), (2) в области _0={(ж, у, ¿) : 0 ^ ж ^ а,
О ^ у ^ Ь, 0 ^ ^ Т} введём сетку
'г^Нх,Ну,т — (Жго) 2Лти ^к) ■ %п — П'Нху Ут — ТпНу, — кт,
п = 0,1,..., И, кх = 1/Л^; т = 0,1,..., М, Ну = 1/М;
к = 0,1,... К, т =
с шагом ¡гх по ж, Ну по у и г по Обозначим ~ и(1к,хп,ут), икп_^+1 та
~ и(1^к1 Xn—j-^-l1 Ут) 1 Un_j_~ u{tk,Xn,ym—j-lt-\)l /п,т ~ /(£&) Ут)- ВОСПОЛЬ-зовавшись равенствами (4) и (5), для уравнения (1) запишем явную и неявную разностные схемы:
(/с+1) (к) ^ / п+1
ип,т ~ (Шп,т _ Сх ( (к) _ {к) ул !к) \
1 ггг, И ип,т ' / ; 4.3 ап— 7 + 1,гтг }'
Г(2 — <у)та ^\п+1’т
ґ-ч / т+1 ч
+ -ф ~ Рип)п + ^2 ^піг-і+і) + ^
У j=2
(/с+1) (/с) / п+1
Цга,ш -дЦп,т Оя/ (А-Ц) _ о (к+1) , V«-7,^+1) Ї +
Г(2 — а)га /г/3 \ га+1>т ^ап,т
/^ї / т+1 \
, Ч/ / (М-1) _ о (к+1) , V о-7У(;г+1) I + (7
,/З I га,т,+ 1 Уап,т ' / у Щип,т-^+1 ) ' <1п,т• Vі
Ґіу \ _о /
Ч ' І=2
Теорема 1. Явная разностная схема (6) устойчива, когда
т»(С^ + 9у) а + 1
(2 + /3)Г(2 — а) ’
где 0 < а: ^ 1, 1 < /3 ^ 2.
Доказательство. Решение -и^1 можно представить в виде ик+1 = = 5'и;г + Г(2 — а)та/к, где 5 —оператор перехода с одного временного слоя на другой. Условием устойчивости по начальным данным разностной схемы (6) является ||5|| ^ 1. Действительно, ||м;г+1|| = 1154^11 ^ Ц^Ц • Цг^Ц или ||м;г+1|| ^
^ \\ик\\ ^ ^ ||«°||-
Следовательно, начальные возмущения затухают. Условие Цб'Ц ^ 1 есть условие того, что спектр оператора лежит внутри круга единичного радиуса на комплексной плоскости. Это означает, что тах|А| < 1, где Л — собственное число оператора перехода в.
Для нахождения собственных значений оператора перехода решение и$т представим в виде возмущения:
ип)п = ехр(шхпкх) ехр(шутку), (8)
где г — мнимая единица, Л — собственные числа оператора перехода [21]. В результате подстановки равенства (8) в (6) определяются собственные значения оператора перехода:
Л<а-Г<2-^г°(48ш^+/і-2)-
_ Г(2-а)С,т°/48Іп2^,+/;_2ч (9)
Ну
Тогда из неравенства (9) получим, что в случае
все собственные числа оператора перехода по абсолютной величине не превосходят 1. Условие (10) равносильно условию теоремы. □
Теорема 2. Неявная разностная схема (7) абсолютно устойчива.
Доказательство. Решение ик+1 можно представить в виде Аик+1 = = ик + Г(2 — а)та/к, тогда 5 = А~1 — оператор перехода с одного временного слоя на другой. Условием устойчивости по начальным данным разностной схемы (7), как было отмечено в доказательстве теоремы 1, является шах |А| < 1, где Л — собственное число оператора перехода 5.
Для доказательства теоремы ещё раз решение и^}т представим в виде возмущения (8). В результате подстановки равенства (8) в (7) для модуля перехода получим выражение
\ ^ (л I Г(2 “ а)С^а (л • 2 Uxhx \
А ^ а I Н----------р--------------------------(4 sm —-1- /3 - 2 J -
V hx
Г(2 - а)Сутс
Из неравенства (11) следует, что неявная разностная схема (7) абсолютно устойчива, так как шах|Л| < 1 при любых шагах т, Нх, Ну. □
Разлагая функции ик^, и^1, икик^, ик^\ в ряд Тейлора
и подставляя полученные соотношения в разностную схему (6), получим
^(ік,хг,у3) = Сх^(ік,хг,у3)+
д^и
+ Cy-^^(tkXi,yj) + fkj + а(т2 ") + b(hx) + c(hy),
где а, b, с — некоторые константы. Следовательно, разностная схема (6) аппроксимирует уравнение (1) с порядком 2 —о; по времени и вторым порядком по координатам х, у.
В случае квадратной сетки, т. е. когда hx = hy = h, Сх = Су = С, для разностной схемы (6) имеем
/ п+1
и{п^] = (а- 27/3)ugL + l{Unll,m + <l3un-j+l,m+
V 3=2
m+1 •.
+ unln+1 + ^2 Cl3Un}n-j+1 ) + Г(2 “ (12)
o —9 /
3=2
uo% = ¥>(0, Уш), u['ff = <p(a, ym),
unl = V(%n,0), u^M = <p(xn, b), (13)
где 7 = Г(2 — а)таС ¡hP.
Un)m — Ф{хт, Уп)і
Явная разностная схема (12) устойчива, когда
(.a + l)h,P
т €
2 • (2 + /3)Г(2 — а)С ’
где 0 < о: ^ 1, 1 < /3 ^ 2.
В случае двух и более пространственных переменных применение неявных схем вызывает трудности, связанные с решением систем линейных уравнений.
Снова рассмотрим квадратную сетку, но теперь для неявной разностной схемы (7). В этом случае (7) примет вид
(к-\-1) (к) ^ / п-\-1
ип,т ~ ОШп'Гп С ( (к+1) _ о.Ак+1) , 0:+1) ,
ап-\-1,т Уап,т ' / у Чз 7 + 1,т '
гтг+1
Г(2 — а)та h13 ч 2
J=2
Решение можно найти с помощью дробных шагов, то есть переход с fc-того слоя на (к + 1)-й разобьём на два полушага:
(к-\-1/2) (к) / п-\-1 \
Un,m ~ aUn,m С_ ( (fc+1/2) _ о (fc+1/2) (fc+i/2) \ , .
Т(2—а)та hi3 \ ra+l>m Z_^43an-j+l,mJ ’ v±J/
(fc+1) (fc) ^ / m+1 4
•Un,m ~ «%,ш W (fc+1) _ о (fc+1) , V- 1 + f(fc) Пб)
— (х}та \ ra>m+1 n,m ' / y i13an,m-j+1 у ^ Jn,m- \±и)
Разностные схемы (15) и (16) являются неявными и содержат неизвестные значения на (к + 1/2)-м шаге в схеме (15) и на (к + 1)-м — в схеме (16). При этом в схеме (15) неизвестные значения функции вычисляются прогонкой по слоям у = const, в схеме (16)—прогонкой по слоям х = const.
Сеточные уравнения (15) и (16) не аппроксимируют уравнение (1) на дробных шагах, но результирующая схема для перехода с fc-того слоя на (к + 1)-й является аппроксимирующей.
Покажем, что модуль перехода, соответствующий двум полушагам, действительно аппроксимирует модуль перехода для уравнения (1). С помощью метода гармоник Фурье [21] можно получить выражения для модулей перехода разностных схем (15) и (16):
( 4Г(2 — а)таС соф о/с^Ллл-1
Л Ч1+ /,/ -г8Ш (-г)) ’
и ( 4Г(2 — а)таС . 2(U)ih\\-i
л <Ч1+ 5Ш (—)) •
Для совокупности переходов (15) и (16) модуль перехода Л = . Таким
образом, схема (15), (16) абсолютно устойчива.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Нигматуллин Р. Р. Дробный интеграл и его физическая интерпретация // ТМФ, 1992. — Т. 90, №3. — С. 354-368; англ. пер.: Nigmatullin R. R. Fractional integral and its physical interpretation// Theoret. and Math. Phys., 1992. — Vol. 90, No. 3. —P. 242-251.
2. Головизнин В. М., Киселев В. П., Короткий И. А., Юрко Ю. И. Прямые задачи неклассического переноса радионуклидов в геологических формациях // Изв. РАН. Энергетика, 2004. - №4. - С. 121-130.
3. Самко С.Г., Килбас А. А., Маричев О. И. Интегралы и производные дробного порядка и некоторые их приложения. — Минск: Наука и техника, 1987. — 498 с.
4. Кольцова Э. М., Василенко В. А., Тарасов В. В. Численные методы решения уравнений переноса во фрактальных средах // Журн. физ. химии, 2000. — Т. 74, № 5. — С. 954-956.
5. Головизнин В. М., Киселев В. П., Короткий И. А., Юрко Ю. И. Некоторые особенности вычислительных алгоритмов для уравнения дробной диффузии: Препринт № IBRAE-2002-01. - М.: ИБРАЭ РАН, 2002. - 57 с.
6. Головизнин В. М., Киселев В. П., Короткий И. А. Численные методы решения уравнения дробной диффузии в одномерном случае: Препринт № IBRAE-2002-10. — М.: ИБРАЭ РАН, 2002. - 35 с.
7. Головизнин В. М., Короткий И. А. Методы численных решений некоторых одномерных уравнений с дробными производными// Дифференц. уравнения, 2006. — Т. 42, №7. — С. 907-913; англ. пер.: Goloviznin V.M., Korotkin I. A. Numerical methods for some one-dimensional equations with fractional derivatives 11 Differ. Equ., 2006. — Vol. 42, No. 7. — P. 967-973.
8. Meerschaert М. М., Tadjeran C. Finite difference approximations for two-sided space-frac-tional partial differential equations// Appl. Numer. Math., 2006. — Vol. 56, No. 1. — P. 80-90.
9. Meerschaert М. М., Tadjeran C. Finite difference approximations for fractional advection-dispersion flow equations// J. Com,put. Appl. Math., 2004. — Vol. 172, No. 1. — P. 65-77.
10. Tadjeran C., Meerschaert М. М., Scheffler H.-P. A second-order accurate numerical approximation for the fractional diffusion equation// J. Comput. Phys., 2006. — Vol. 213, No. 1. — P. 205-213.
11. Lynch V.E., Carreras B.A., del-Сastill-Negrete D., Ferreira-Mejias K.M., Hicks H. R. Numerical methods for the solution of partial differential equations of fractional order //
J. Comput. Phys., 2003. — Vol. 192, No. 2. — P. 406-421.
12. Liu Q., Liu F., Turner I., Anh V. Approximation of the Levy-Feller advection-dispersion process by random walk and finite difference method // J. Comput. Phys., 2007. — Vol. 222, No. 1. — P. 57-70.
13. Таукенова Ф. И., Шхануков—Лафишев М. X. Разностные методы решения краевых задач для дифференциальных уравнений дробного порядка // Ж. вычисл. матем. и ма-тем. физ., 2006. — Т. 46, № 10. — С. 1871-1881; англ. пер.: Taukenova F. Shkhanukov-Lafishev М. Kh. Difference methods for solving boundary value problems for fractional-order differential equations// Comput. Math. Math. Phys., 2006. — Vol. 46, No. 10. — P. 1785-1795.
14. Лафишева M.M., Шхануков—Лафишев М. X. Локально-одномерная разностная схема для уравнения диффузии дробного порядка// Ж. вычисл. матем. и матем. физ., 2008. — Т. 48, № 10. — С. 1878-1887; англ. пер.: Lafisheva М.М., Shkhanukov-Lafishev M.Kh. A locally one-dimensional difference scheme for a fractional-order diffusion equation // Comput. Math. Math. Phys., 2008. — Vol. 48, No. 10. — P. 1875-1884.
15. Алиханов А. А. Разностные методы решения краевых задач для волнового уравнения с дробной производной по времени// Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2008. - №2(17). - С. 13-20.
16. Бейбалаев В. Д. Математическая модель теплопереноса в средах с фрактальной структурой// Матем. моделирование, 2009. — Т. 21, №5. — С. 55-62.
17. Бейбалаев В. Д. Численный метод решения задачи переноса с двусторонней производной дробного порядка // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2009. — Т. 1(18). - С. 267-270.
18. Назаралиев М. А., Бейбалаев В. Д. Численные методы решения краевой задачи для уравнения теплопереноса с производной дробного порядка// Вестн. ДГУ, 2008. — №6. - С. 46-53..
19. Бейбалаев В. Д. Численный метод решения математической модели теплопереноса
в средах с фрактальной структурой // Фундаментальные исследования, 2007. — № 12. — С. 249-251.
20. Нахушев А. М. Элементы дробного исчисления и их применение. — М.: Физматлит, 2003. - 272 с.
21. Самарский А. А. Теория разностных схем: 2-е изд. — М.: Наука, 1983. — 616 с.
Поступила в редакцию 08/11/2010; в окончательном варианте — 28/1У/2010.
MSC: 65N12, 34A08, 26A33, 45K05
NUMERICAL METHOD OF VALUE BOUNDARY PROBLEM DECISION FOR 2D EQUATION OF HEAT CONDUCTIVITY WITH FRACTIONAL DERIVATIVES
V.D. Beybalaev1, M. V. Shabanova2
1 Daghestan State University,
43 a, M. Gadzhiyev St., Makhachkala, 367025, Russia.
2 Institute of Geothermy Problems, Dagestan Research Center of RAS,
39 a, pr. Shamilya, Makhachkala, 367030, Russia.
E-mails: kaspij_03@mail.ru, meil-sha8yandex.ru
In this work a solution is obtained for the boundary problem for two-dimensional thermal conductivity equation with derivatives of fractional order on time and space variables by grid method. Explicit and implicit difference schemes are developed,. Stability criteria of these difference schemes are proven. It is shown that approximation order by time equal but by space variables it equal two. A solution method is suggested using fractional steps. It is proved that the transition module, corresponding to two half-steps, approximates the transition module for given equation.
Key words: numerical methods, stability, approximation of fractional derivatives, fractional differential equations.
Original article submitted 08/11/2010; revision submitted 28/IV/2010.
Vetlugin D. Beybalaev (Ph. D. (Phys. & Math.)), Senior Lecturer, Dept, of Applied Mathematics. Mumina R. Shabanova, Research Assistant, Lab. of Mathematical Modeling & Monitoring of Geothermal Objects.
