CC BY
32
DOI: 10.24143/2072-9502-2019-4-70-80 УДК 504.054
ЧИСЛЕННОЕ РЕШЕНИЕ НЕСТАЦИОНАРНОГО ДРОБНОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ В ЗАДАЧАХ МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТОКСИЧНЫХ ВЕЩЕСТВ В ПОДЗЕМНЫХ ВОДАХ
А. А. Афанасьева, Т. Н. Швецова-Шиловская, Д. Е. Иванов, Д. И. Назаренко, Е. В. Казарезова
Государственный научно-исследовательский институт органической химии и технологии,
Москва, Российская Федерация
В настоящее время во многих областях науки для моделирования различных процессов широко применяется теория дробного исчисления. Дифференциальные уравнения с производными дробного порядка используются при моделировании миграции токсичных веществ в пористых неоднородных средах и позволяют более корректно описать поведение токсичных веществ на больших расстояниях от источника. Аналитическое решение дифференциальных уравнений с производными дробного порядка зачастую оказывается очень сложным или даже невозможным. Предложен численный метод решения дифференциальных уравнений в частных производных дробного порядка по времени для описания миграции токсичных веществ в подземных водах. Для численного решения нестационарного дробного дифференциального уравнения разработана неявная разностная схема, являющаяся аналогом известной неявной разностной схемы Кранка - Николсона. Система разностных уравнений представлена в матричном виде. Решение задачи сводится к многократному решению трехдиагональной системы линейных алгебраических уравнений методом прогонки. Представлены результаты оценки распространения токсичного вещества в подземных водах на основе численного метода для модельных примеров. Выполнено сравнение концентраций вещества, полученных на основе аналитического и численного решения нестационарного одномерного дробного дифференциального уравнения. Результаты, полученные с помощью предлагаемого метода и на основании известного аналитического решения дробного дифференциального уравнения, достаточно хорошо согласуются между собой. Относительная ошибка составляет в среднем 9 %. В отличие от известного аналитического решения разработанный численный метод может использоваться при моделировании распространения токсичных веществ в подземных водах с учетом их биодеградации.
Ключевые слова: дробная производная, производная по времени, дробное дифференциальное уравнение, разностная схема, подземные воды.
Для цитирования: Афанасьева А. А., Швецова-Шиловская Т. Н., Иванов Д. Е., Назаренко Д. И., Казарезова Е. В. Численное решение нестационарного дробного дифференциального уравнения в задачах моделирования распространения токсичных веществ в подземных водах // Вестник Астраханского государственного технического университета. Серия: Управление, вычислительная техника и информатика. 2019. № 4. С. 70-80. DOI: 10.24143/2072-9502-2019-4-70-80.
Введение
В настоящее время многих исследователей интересует использование математических моделей с дробными производными в силу точности моделей при описании объектов в различных областях науки. Причинами возросшего интереса к дробным производным и дробным дифференциальным уравнениям являются следующие: 1) определение дробной производной нелокально и может быть использовано для математического моделирования сред с памятью; 2) дробные производные и дифференциальные уравнения дробных порядков точнее описывают рассматриваемые явления во многих задачах физики, биологии, геологии; 3) дробными дифференциальными уравнениями можно описывать немарковские процессы. Повышенный интерес к дробному исчислению, наблюдаемый в последнее время, позволяет говорить о нем как о новой области науки [1-10].
Особенный интерес к дробным производным проявляют гидрогеологи в связи с проблемой обоснования безопасности хранения опасных химических и радиоактивных отходов в геологических формациях [1, 2]. В соответствии с имеющимися экспериментальными данными изменение концентрации вещества на больших расстояниях от источника не может быть описано законами классической диффузии, поскольку подчиняется степенному закону убывания, а не экспоненциальному, как это считалось ранее. Данное явление связывают с тем, что область распространения подземных вод представляет собой неоднородную гетерогенную среду, которая включает различные по литологическому составу (по водопроницаемости) породы [1, 2].
Как указано в [1, 2], в таких случаях имеет место так называемая аномальная диффузия, которую подразделяют на субдиффузию и супердиффузию. Субдиффузия протекает намного медленнее классической диффузии, поскольку диффузант захватывается ловушками (дефектами, адсорб-ционно- или химически активными центрами), уходит в боковые, тупиковые пути и на некоторое время или навсегда выводится из миграционного процесса. Наоборот, супердиффузия осуществляется со скоростями, существенно превышающими классическую диффузию. Этот режим наблюдается, если в системе есть облегченные пути (например, трещины) или присутствуют процессы случайной или направленной адвекции (увлечение диффузанта потоками флюидов). Эти явления невозможно корректно описать классическим уравнением турбулентной диффузии.
С учетом этого обстоятельства для более точного описания данных процессов необходимо использовать уравнения диффузии в частных производных дробных порядков [1, 2, 10]. Чтобы учесть в математической модели структуру среды или изменения внешнего поля со временем, нужно заменить обычный оператор производной оператором производной вариационного порядка.
Дифференциальные уравнения в частных производных дробных порядков подразделяют на два вида: дифференциальные уравнения с дробной производной по пространству и дробной производной по времени. Чаще всего получить решение для таких уравнений в явном виде не удается. В настоящее время известны аналитические решения дробных дифференциальных уравнений для ограниченного круга задач. Например, в работе [11] приведено аналитическое решение дробного дифференциального уравнения, которое описывает распространение вещества в подземных водах без учета его возможных химических превращений. Поэтому для решения дробных дифференциальных уравнений в частных производных очень эффективными являются приближенные методы [3].
Несмотря на развитый математический аппарат дробного интегрирования и дифференцирования, вопросы численного решения уравнений с дробными производными практически не нашли отражения в существующей литературе [1, 2, 12-15].
В данной работе представлены результаты исследований по разработке численного метода и разностной схемы для решения одномерного дифференциального уравнения с дробной производной по времени.
Наибольшее распространение в практике гидрогеологических расчетов получил метод конечных разностей. Однако данный метод обладает некоторыми недостатками [4, 5]. Избавиться от недостатков данного метода можно, используя различные подходы к построению разностной схемы [7].
В данной статье рассматривается уравнение параболического типа, описывающее процесс диффузии, и один из подходов к построению разностных схем для такого типа уравнений -схема Кранка - Николсона, называемая также схемой с весом 1/2. Устойчивость и сходимость решения конечно-разностных уравнений на основе схемы Кранка - Николсона доказана во многих работах [1, 2]. Применение данной схемы позволяет без существенных вычислительных затрат повысить (по сравнению с другими подходами) порядок аппроксимации разностной схемы [13-19].
Таким образом, наиболее целесообразной является разработка неявной разностной схемы для решения дробного дифференциального уравнения на основе известной схемы Кранка - Николсона.
Неявная разностная схема для численного решения нестационарного одномерного дробного дифференциального уравнения
Рассмотрим общее одномерное уравнение для моделирования миграции токсичных веществ в подземных водах с начальными и граничными условиями вида
дса (х, t) „д с dc .
= D—- -v— - лс,
дх2
dtа
дх
0 < х < L, t > 0,
(1)
начальные условия: с(х, 0) = с0, 0 < х < L; граничные условия: с (0, 0 = 0, с (Д {) = 0, t > 0; где с - концентрация токсичного вещества, мг/м3; с0- начальная концентрация токсичного вещества, мг/м3; D - коэффициент дисперсии в направлении оси х, м2/с; X - константа скорости превращения токсичного вещества, 1/с; V - скорость движения воды в порах почвы в направлении оси х, м/с; I - время, с; L - расстояние, на котором концентрация принимает значение с0, м; а - порядок производной.
В 1967 г. М. Капуто ввел новое и полезное определение дробной производной, которое сегодня широко известно. Приведем это определение для дробной производной по времени порядка а [20, 21]:
дас( х, t) д at
1
г дс( х, t) dt, Г(1 - а )| (Г-^
дс (х, t)
0 < а < 1,
(2)
а = 1,
где Г (.) - Гамма-функция.
Разработаем неявную разностную схему, являющуюся аналогом схемы Кранка - Николсона, для численного решения дробного дифференциального уравнения. Схема Кранка - Николсона может быть записана в виде
дс дс д2с.
— = Г (с, х, t, —, —г); д1 дх дх
1
дс д2с,
ck+1 - ck , ,,, д с
с-CJ- =±. [Flk+1(c, х, tЦ-) + Fk (с, х, t, —
At 2 дх дх2 дх дх
дс д2с
(3)
В качестве сетки возьмем совокупность точек (х1, tk) , координаты которых х1 = I ■ h; I = 0, 1, к, М; 0 < I < М ; М ■ И = L; ^ = к ■ т; 0 < к < N ; N ■ т = Т. На данной сетке будем строить разностную аппроксимацию дифференциального уравнения (1), используя шеститочечный
шаблон: Ш {(хг, tк ), (x/+l, tk+l ),(х1 -1 , tk+1
Следовательно, схема Кранка - Николсона для модели (1) и шеститочечного шаблона имеет вид
дс дх 2
1 ff с( хмД+1) - с(х1 -1tk+1) 1+f с(х1+1tk) - с( )
W
2(h)
2(h)
+ O(h);
//
д2с 1 f( с(х,+1tk+1) - 2с(х, tk+1) + с(х1 -1tk+1) 1 + f с(х1+1tk) - 2с(х1 h) + с(х1 -1tk) Ц
дх2
w
(h)2
(h)2
+ O (h2);
//
с = 1(с( хltk+1) + с( х, tk)) •
(4)
(5)
(6)
Приближенное решение задачи (1) получаем, заменяя значение исходной непрерывной функции на функцию, определяемую в узлах сетки (к - номер узла по времени и I - номер узла
по пространству): = с(х{, tк) и С означает численную аппроксимацию точного решения с(хг, tk ) (или обозначает численное приближение к точному решению).
Аналог неявной разностной схемы Кранка - Николсона для численного решения дробного дифференциального уравнения (1) с учетом выражения для дробной производной Капуто по времени (2) можно сформулировать следующим образом:
да1 с(х, ,tfe+1)
Г(2 - аг )
с(х,,tk+1)-с(х,,tk)+У[с(х,,tk+1-j)-с(х,,tk-j)][j+1)1-а,+ -С/)1-а1+ ] к (7)
^^L V l
j=1
k
а
k
T
Выполним последовательно преобразования правой части уравнения (7):
т-а'+1 Г к 1
——к— \е (х,, tk+l) - с (х,, tk) + 1 [с (x, ) - с (х,, )][(з + 1)1-а<+' - (;)1-а"1 ]1 = 1(2 - аг ) I з=1 ]
Т - а 1+1 т - а к „ к
(8)
с+1 - ск) + ,; к+и I(ск+1-з - ск-з)[(з+1)1-а' - (,)1-а- ].
Г(2-ак+1) ' ' 1(2-ак+1)
г ) ,=1
Запишем уравнение (1) с учетом (4)-(6), (8):
(с-"1 -ск) 2(ск+1-3 -ск-3Ц*1 =
Г(2 - с^Г' '' Г(2 - аГ) ^ ' (9)
= D - 2ск+1 + + с+1 - 2ск + ск-1)(<% - ск+11 + с+1 - £)-к 1(сгк+1 + ск),
7 7,1 1 к+1 1 к+1
где +1 = (з + 1)1-аг -(,)1-аг ; з = 1, 2, ..., N; I = 1, 2, ...,М -1.
Проведем упрощающие преобразования. Примем, что коэффициент перед выражением (с\+1 - ск) в левой части равен единице, и введем следующие обозначения:
Ьк+1 = 1 Г(2 - оО = DГ(2 - ак+1)хак+' ,
2к2 х-аГ 2к
§к1 =
к+1_ V Г(2 - ак+')таг , (10)
4к
ак+1 = -к Г(2 - ак+1)так+' г 2 . После преобразований (8)-(10) система разностных уравнений для первоначальной краевой задачи с дробной производной (1) имеет вид
ьк+1 [(ск:; - 2ск+1+ск--;)+(ск+1 - 2ск+ск-1)]+£+1 кск? - с-1)+
к - ^. (11)
г-1Л 1 иг 1 ^г ^г ^г 1 IV-г с'
3=1
к км, ,]к+1/ к+1 , к\ к+1 к , X ' ( к+1-, к-/\ г,к+1 + (сг+1 - с1 -1)] + Лг (сг + сг ) = с1 - С + 1 (сг 3 - сг3 Н'
Отметим, что при к = 0 слагаемое с суммой исчезает, что упрощает реализацию метода при проведении вычислений.
Выполним преобразование уравнения (11):
с-+1 - Ьгк+1(ск++11 - 2ск+1 + ск_+1) - +1(ск++11 - О - Л,к+1С1к+1 =
к
-¡к+\/к л к . к\-|..к+1/к к \ , лк+\ к . к X""1 / к+1- / к- / \ г,к+1
= Ьг (с-+1 - 2сг + сг-1)] + &г (сг+1 - сг-1) + Л сг + с, -I(С- с- ,
/=1
или
<+/( - ьк+1 + ягк+1) + ск+1(1 + 2 ьк+1 - л к+1) + сгк++11( - ьк+1 - ^гк+1) =
сгк-1( ьк+1 - gk+1) + ск (1 - 2 ьк+1 + л к+1) + сгк+1( ьк+1 + gk+1) -I (ск+1-з - сгк - з )юг/+1
з=1
Коэффициенты в узлах сетки вычисляются по следующим формулам:
к+1 г.к+1 . „к+1 к 7 к+1 „к+1 к+1 1 . 1 . ог.к+1 лк+1
сг-1 : пг =-ьг + gг ; сг-1 : - пг = ь - gг ; сг :1 + даг = 1 + 2ь, - Лг ;
+1 +1 +1 +1 +1 +1 +1 сг : 1 - даг = 1 - 2ьг + ; сг+1 : Р1 =-ь, - gl ; сг+1 : - Рг = ь + gl •
(12)
А:
-а
т
Систему разностных уравнений (12) с учетом преобразований можно представить в следующем виде:
п1е1+11 + ск+1(1 + т,) + с£( - Рг) =
k / \ , k /л \ , k Ч 1 / k+1- j k - /\ l ,k+1
= cl-l(-ni) + ci (1 - ml) + cl+1 Pl (cl - cl )ю
или
щс^1 + cf+1(1 + ml) + С/( - P,) =
j=1
k+1.
k-1
ck-1(-n) + ck(1 - m! - +1) + cl.p, + c>k'k+1 - X ck-j(q j
^-'k+1 ,k+1 ,ю j - ю j+1
)
j=1
Первоначальная краевая задача с дробными производными (2) может быть представлена в матричном виде:
АС1 = ВС0 для к = 0; (13)
k-1
ACk+1 = BCk + X (юj - юj+1)Ck-1 + юС0 для k > 1
j=1
(14)
где Ск = (С1к, Ск, к , СМ-1)Т, к = 0,1,2, ... ; начальные условия и граничные условия:
с? = с0; I = 1, 2,м; с0к = 0; скм = 0.
Матрицы А и 5 имеют следующий вид:
f(1+m) - P1
А =
0
(1 + m2) - P2
(1 +m) - p,
lM -
1 (1 + mM-1 )y
'(1 - ю1 - m1) n
B =
P1 K
(1 - Ю1 - m2) P2
(M-1 )x(M-1); ... 0
(1 - ю- m)
- n
M-1
(1 - ю- mM-1); 1
(M-1 )x(M-1).
n
2
0
0
n
0
0
0
n
0
Решение задачи (13), (14) в конечном счете сводится к многократному решению системы линейных алгебраических уравнений (СЛАУ) с трехдиагональной матрицей. Для решения таких СЛАУ известны специальные методы, например метод прогонки и метод циклической редукции [9]. Наиболее простым и эффективным методом является метод прогонки, который привлекателен своей экономичностью и простотой реализации в последовательном режиме [8]. Именно поэтому метод прогонки использовался в данном случае для решения системы линейных алгебраических уравнений.
Разработанный численный метод решения дробного дифференциального уравнения реализован в программном комплексе [22] для решения задачи распространения токсичных веществ в условиях аномальной диффузии. Алгоритм реализации численного метода представлен на рис. 1.
Пользователь принимает решение о возвращении на один из предыдущих этапов ввода исходных данных
Решение СЛАУ
k-i
ac+1 = bc + -öji) c-1+щс°
j=i
Рис. 1. Алгоритм реализации разработанного численного метода
Данный алгоритм предусматривает проверку корректности исходных и расчетных данных на всех этапах вычислений.
Новый численный метод реализован на языке программирования Microsoft Visual Basic 6.5, в качестве системы управления базой данных использовался Microsoft Access.
Оценка распространения токсичного вещества в подземных водах на основе численного метода
В данном разделе приведены результаты численного эксперимента, иллюстрирующие вычислительные свойства нового метода. Выполнено сравнение концентраций вещества, полученных на основе решения нестационарного одномерного дробного дифференциального уравнения разработанным методом и на основании подхода, представленного в работе [11], для двух значений порядка производной по времени - а.
Результаты приведены для модельного примера, в котором химические превращения токсичного вещества не рассматривались. Параметры расчетной сетки: h = 0,10; т = 0,1. Число расчетных точек по осям - 100.
На рис. 2 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при а = 0,2.
0.4 0.6
Расстояние
_Численное
решение
.. . Аналитическое решение
Рис. 2. Зависимость концентрации токсичного вещества от расстояния (а = 0,2)
На рис. 3 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при а = 0,6.
0.4 0.6
Расстояние
Численное решение
Аналитическое решение
Рис. 3. Зависимость концентрации токсичного вещества от расстояния (а = 0,6)
Из рис. 2 и 3 видно, что графики зависимостей концентрации токсичного вещества от расстояния изменяются симбатно. Результаты, полученные с помощью предлагаемого метода и известного точного решения дробного дифференциального уравнения, достаточно хорошо согласуются между собой. Относительная ошибка составляет в среднем 9 %.
В отличие от известного аналитического решения разработанный численный метод может использоваться при моделировании распространения токсичных веществ в подземных водах с учетом их биодеградации.
Заключение
В результате проведенных исследований разработан новый численный метод решения дифференциального уравнения с производной дробного порядка по времени, а также неявная разностная схема, являющаяся аналогом схемы Кранка - Николсона.
Система разностных уравнений представлена в матричном виде. Решение задачи сводится к многократному решению трехдиагональной системы линейных алгебраических уравнений методом прогонки. Представлены результаты оценки распространения токсичных веществ в подземных водах на основе численного метода при решении модельных задач.
Полученные результаты позволяют рекомендовать разработанный метод численного решения нестационарного одномерного дробного дифференциального уравнения, основанный на использовании линейных интегро-дифференциальных уравнений достаточно общей структуры, для использования на практике при моделировании процессов тепло- и массопереноса в сильно неоднородных средах.
СПИСОК ЛИТЕРА ТУРЫ
1. Головизнин В. М., Киселев В. П., Короткин И. А., Юрков Ю. И. Некоторые особенности вычислительных алгоритмов для уравнений дробной диффузии. М.: Изд-во Ин-та проблем безопасного развития атомной энергетики РАН, 2002. 57 с.
2. Головизнин В. М., Киселев В. П., Короткин И. А. Численные методы решения уравнений дробной диффузии в одномерном случае. М.: Изд-во Ин-та проблем безопасного развития атомной энергетики РАН, 2002. 35 с.
3. Меркулова Н. Н., Михайлов М. Д. Методы приближенных вычислений: учеб. пособие. Томск: Изд-во ТГУ, 2011. Ч. III. 184 с.
4. Самарский А. А., Гулин А. В. Численные методы. М.: Наука, 1989. 432 с.
5. Вержбицкий В. М. Численные методы. М.: Высш. шк., 2001. 382 с.
6. Байков В. А., Жибер А. В. Уравнения математической физики. М.; Ижевск: Изд-во Ин-та компьютерных исследований, 2003. 252 с.
7. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1987. 130 с.
8. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1977. 735 с.
9. Кустикова В. Д. Дифференциальные уравнения в частных производных. Н. Новгород: Изд-во НГУ им. Н. И. Лобачевского, 2011. 54 с.
10. Хенди Ахмед Санд Абделазиз. Численные алгоритмы решения дробных дифференциальных уравнений с запаздыванием: дис. ... канд. физ.-мат. наук. Екатеринбург, 2017. 124 с.
11. Shukla H. S., Tamsir M., Srivastava V. K., Kumar J. Approximate Analytical Solution of Time-fractional order Cauchy-Reaction Diffusion equation // Computer Modeling in Engineering & Sciences. 2014. V. 103. N. 1. Р. 1-17.
12. Самко С. Г., Килбас А. А., Маричев О. И. Интегралы и производные дробного порядка и некоторые их приложения. Минск: Наука и техника, 1987. 687 с.
13. Штумпф C. А., Бахтин М. А. Математические методы компьютерных технологий в научных исследованиях: учеб. пособие. СПб.: Изд-во НИУ ИТМО, 2012. 148 с.
14. Fadugba S. E., Nwozo C. R. Crank-Nicolson Finite Difference Method for the Valuation of Options // The Pacific Journal of Science and Technology. 2013. V. 14. N. 2. Р. 136-145.
15. Wani S. S., Thakar S. H. Crank-Nicolson Type Method for Burgers Equation // International Journal of Applied Physics and Mathematics. 2013. V. 3. N. 5. Р. 324-328.
16. Kilbas A. A., Trujillo J. J. Differential equations of fractional order: methods, results and problems // Applicable Analysis. 2001. V. 78. N. 1. Р. 153-192.
17. Atangana1 А., Kilicman А. On the Generalized Mass Transport Equation to the Concept of Variable Fractional Derivative // Mathematical Problems in Engineering. 2014. Р. 1-9.
18. Ерохин С. В. Математическое моделирование напряженно-деформированного состояния вязко-упругих тел с использованием методов дробного исчисления: дис. ... канд. техн. наук. М., 2016. 130 с.
19. Бэгли Р. Л., Торвик П. Дж. Дифференциальное исчисление, основанное на производных дробного порядка - новый подход к расчету конструкций с вязкоупругим демпфированием // Аэрокосмическая техника. 1984. Т. 2. № 2. С. 84-91.
20. Caputo M. Linear model of dissipation whose Q is almost frequency independent // Geophis. J. R. Astr. Soc. 1967. V. 13. P. 529-539.
21. Caputo M., Mcinardi F. Linear models of dissipation in an elastic solids // Riv. Nuovo Cimento. 1971. V. 1. Р. 161-196.
22. Свидетельство о государственной регистрации программы для ЭВМ № 2019614483. Программное обеспечение для моделирования распространения загрязняющих веществ в подземных водах / Афанасьева А. А., Назаренко Д. И., Швецова-Шиловская Т. Н., Казарезова Е. В.; зарег. 26.03.2019; опубл. 05.04.2019.
Статья поступила в редакцию 28.08.2019
ИНФОРМАЦИЯ ОБ АВТОРАХ
Афанасьева Александра Алексеевна - Россия, 111024, Москва; Государственный научно-исследовательский институт органической химии и технологии; канд. техн. наук; начальник отдела; [email protected].
Швецова-Шиловская Татьяна Николаевна - Россия, 111024, Москва; Государственный научно-исследовательский институт органической химии и технологии; д-р техн. наук, профессор; начальник отделения; [email protected].
Иванов Дмитрий Евгеньевич - Россия, 111024, Москва; Государственный научно-исследовательский институт органической химии и технологии; канд. техн. наук; ведущий научный сотрудник; [email protected].
Назаренко Денис Игоревич - Россия, 111024, Москва; Государственный научно-исследовательский институт органической химии и технологии; канд. техн. наук; начальник отдела; [email protected].
Казарезова Елена Викторовна - Россия, 111024, Москва; Государственный научно-исследовательский институт органической химии и технологии; младший научный сотрудник; [email protected].
NUMERICAL CALCULATION OF NONSTATIONARY FRACTIONAL DIFFERENTIAL EQUATION IN PROBLEMS OF MODELING TOXIC SUBSTANCES DISTRIBUTION IN GROUND WATERS
A. A. Afanasyeva, T. N. Shvetsova-Shilovskaya, D. E. Ivanov, D. I. Nazarenko, E. V. Kazarezova
State Research Institute оf Organic Chemistry and Technology, Moscow, Russian Federation
Abstract. At present, the theory of fractional calculus is widely used in many fields of science for modeling various processes. Differential equations with fractional derivatives are used to model the migration of pollutants in porous inhomogeneous media and allow a more correct description of the behavior of pollutants at large distances from the source. The analytical solution of differential equations with fractional order derivatives is often very complicated or even impossible. There has been proposed a numerical method for solving fractional differential equations in partial derivatives with respect to time to describe the migration of pollutants in groundwater. An implicit difference scheme is developed for the numerical solution of a non-stationary fractional differential equation, which is an analogue of the well-known implicit Crank-Nicholson difference scheme. The system of difference equations is presented in matrix form. The solution of the problem is reduced to the
multiple solution of a tridiagonal system of linear algebraic equations by the tridiagonal matrix algorithm. The results of evaluating the spread of pollutant in groundwater based on the numerical method for model examples are presented. The concentrations of the substance obtained on the basis of the analytical and numerical solutions of the unsteady one-dimensional fractional differential equation are compared. The results obtained using the proposed method and on the basis of the well-known analytical solution of the fractional differential equation are in fairly good agreement with each other. The relative error is on average 9%. In contrast to the well-known analytical solution, the developed numerical method can be used to model the spread of pollutants in groundwater, taking into account their biodegradation.
Key words: fractional derivative, time derivative, fractional differential equation, difference scheme, groundwater.
For citation: Afanasyeva A. A., Shvetsova-Shilovskaya T. N., Ivanov D. E., Nazarenko D. I., Kazarezova E. V. Numerical calculation of nonstationary fractional differential equation in problems of modelling toxic substances distribution in ground waters. Vestnik of Astrakhan State Technical University. Series: Management, Computer Science and Informatics. 2019;4:70-80. (In Russ.) DOI: 10.24143/2072-9502-2019-4-70-80.
REFERENCES
1. Goloviznin V. M., Kiselev V. P., Korotkin I. A., Iurkov Iu. I. Nekotorye osobennosti vychislitel'nykh algoritmov dlia uravnenii drobnoi diffuzii [Characteristics of computational algorithms for fractional diffusion equations]. Moscow, Izd-vo In-ta problem bezopasnogo razvitiia atomnoi energetiki RAN, 2002. 57 p.
2. Goloviznin V. M., Kiselev V. P., Korotkin I. A. Chislennye metody resheniia uravnenii drobnoi diffuzii v odnomernom sluchae [Numerical methods for solving equations of fractional diffusion in one-dimensional case]. Moscow, Izd-vo In-ta problem bezopasnogo razvitiia atomnoi energetiki RAN, 2002. 35 p.
3. Merkulova N. N., Mikhailov M. D. Metody priblizhennykh vychislenil: uchebnoe posobie [Approximate computation methods: tutorial]. Tomsk, Izd-vo TGU, 2011. Part III. 184 p.
4. Samarskii A. A., Gulin A. V. Chislennye metody [Numerical procedures]. Moscow, Nauka Publ., 1989. 432 p.
5. Verzhbitskii V. M. Chislennye metody [Numerical procedures]. Moscow, Vysshaia shkola Publ., 2001. 382 p.
6. Baikov V. A., Zhiber A. V. Uravneniia matematicheskoi fiziki [Equations of mathematical physics]. Moscow, Izhevsk, Izd-vo In-ta komp'iuternykh issledovanii, 2003. 252 p.
7. Samarskii A. A., Nikolaev E. S. Metody resheniia setochnykh uravnenii [Methods for solving grid equations]. Moscow, Nauka Publ., 1987. 130 p.
8. Tikhonov A. N., Samarskii A. A. Uravneniia matematicheskoi fiziki [Equations of mathematical physics]. Moscow, Nauka Publ., 1977. 735 p.
9. Kustikova V. D. Differentsial'nye uravneniia v chastnykh proizvodnykh [Partial differential equations]. Nizhnii Novgorod, Izd-vo NGU im. N. I. Lobachevskogo, 2011. 54 p.
10. Khendi Akhmed Sand Abdelaziz. Chislennye algoritmy resheniia drobnykh differentsial'nykh uravnenii s zapazdyvaniem. Dissertatsiia ... kand. fiz.-mat. nauk [Numerical algorithms for solving fractional differential equations with delay. Diss. ... Cand.Phys.-Math.Sci.]. Ekaterinburg, 2017. 124 p.
11. Shukla H. S., Tamsir M., Srivastava V. K., Kumar J. Approximate Analytical Solution of Time-fractional order Cauchy-Reaction Diffusion equation. Computer Modeling in Engineering & Sciences, 2014, vol. 103, no. 1, pp. 1-17.
12. Samko S. G., Kilbas A. A., Marichev O. I. Integraly i proizvodnye drobnogo poriadka i nekotorye ikh prilozheniia [Fractional integrals, derivatives and their applications]. Minsk, Nauka i tekhnika Publ., 1987. 687 p.
13. Shtumpf C. A., Bakhtin M. A. Matematicheskie metody komp'iuternykh tekhnologii v nauchnykh issle-dovaniiakh: uchebnoe posobie [Mathematical methods of computer technology in research: teaching guide]. Saint-Petersburg, Izd-vo NIU ITMO, 2012. 148 p.
14. Fadugba S. E., Nwozo C. R. Crank-Nicolson Finite Difference Method for the Valuation of Options. The Pacific Journal of Science and Technology, 2013, vol. 14, no. 2, pp. 136-145.
15. Wani S. S., Thakar S. H. Crank-Nicolson Type Method for Burgers Equation. International Journal of Applied Physics and Mathematics, 2013, vol. 3, no. 5, pp. 324-328.
16. Kilbas A. A., Trujillo J. J. Differential equations of fractional order: methods, results and problems. Applicable Analysis, 2001, vol. 78, no. 1, pp. 153-192.
17. Atangana1 A., Kilicman A. On the Generalized Mass Transport Equation to the Concept of Variable Fractional Derivative. Mathematical Problems in Engineering, 2014, pp. 1-9.
18. Erokhin S. V. Matematicheskoe modelirovanie napriazhenno-deformirovannogo sostoianiia viazkoupru-gikh tel s ispol'zovaniem metodov drobnogo ischisleniia. Dissertatsiia ... kand. tekhn. nauk [Mathematical modeling of stress-strain state of viscoelastic bodies using fractional calculus methods. Diss. ... Cand.Tech.Sci.]. Moscow, 2016. 130 p.
19. Begli R. L., Torvik P. Dzh. Differentsial'noe ischislenie, osnovannoe na proizvodnykh drobnogo poriadka -novyi podkhod k raschetu konstruktsii s viazkouprugim dempfirovaniem [Differential calculus based on fractional derivatives as a new approach to calculation of structures with viscoelastic damping]. Aerokosmicheskaia tekhnika, 1984, vol. 2, no. 2, pp. 84-91.
20. Caputo M. Linear model of dissipation whose Q is almost frequency independent. Geophis. J. R. Astr. Soc, 1967, vol. 13, pp. 529-539.
21. Caputo M., Mcinardi F. Linear models of dissipation in an elastic solids. Riv. Nuovo Cimento, 1971, vol. 1, pp. 161-196.
22. Afanas'eva A. A., Nazarenko D. I., Shvetsova-Shilovskaia T. N., Kazarezova E. V. Programmnoe obespechenie dlia modelirovaniia rasprostraneniia zagriazniaiushchikh veshchestv v podzemnykh vodakh [Software for modeling distribution of pollutants in groundwater]. Svidetel'stvo o gosudarstvennoi registratsii programmy dlia EVM № 2019614483 ot 26.03.2019. opubl. 05.04.2019.
The article submitted to the editors 28.08.2019
INFORMA TION ABOUT THE A UTHORS
Afanasyeva Alexandra Alekseyevna - Russia, 111024, Moscow; State Research Institute of Organic Chemistry and Technology; Candidate of Technical Sciences; Head of the Department; [email protected].
Shvetsova-Shilovskaya Tatyana Nikolayevna - Russia, 111024, Moscow; State Research Institute of Organic Chemistry and Technology; Doctor of Technical Sciences; Professor; Head of the Department; [email protected].
Ivanov Dmitriy Evgenevich - Russia, 111024, Moscow; State Research Institute of Organic Chemistry and Technology; Candidate of Technical Sciences; Leading Researcher; [email protected].
Nazarenko Denis Igorevich - Russia, 111024, Moscow; State Research Institute of Organic Chemistry and Technology; Candidate of Technical Sciences; Head of the Department; [email protected].
Kazarezova Elena Victorovna - Russia, 111024, Moscow; State Research Institute of Organic Chemistry and Technology; Candidate of Technical Sciences; Junior Researcher; [email protected].
