Translate this page:
Please select your language to translate the article


You can just close the window to don't translate
Library
Your profile

Back to contents

Cybernetics and programming
Reference:

On the influence of the method of approximation of an unknown function on the stability of numerical methods for solving the anomalous diffusion equation

Litvinov Vladimir Andreevich

PhD in Physics and Mathematics

Docent, the department of Informatics and Special Technology, Barnaul Law Institute of the Ministry of Internal Affairs of Russia

656038, Russia, Altaiskii krai, g. Barnaul, ul. Chkalova, 49

lva201011@yandex.ru
Other publications by this author
 

 

DOI:

10.25136/2644-5522.2019.2.29201

Received:

11-03-2019


Published:

27-05-2019


Abstract: The subject of the research is numerical algorithms for solving fractional partial differential equations. The object of the study is the stability of several algorithms for the numerical solution of the anomalous diffusion equation. Algorithms based on the difference representation of the fractional Riemann-Liuville derivative and the Caputo derivative for various orders of accuracy are considered. A comparison is made of the results of numerical calculations using the analyzed algorithms for a model problem with the exact solution of the anomalous diffusion equation for various orders of the fractional derivative with respect to the spatial coordinate. The results of the work were obtained on the basis of the analysis of the constructed difference schemes for stability, the conducted numerical experiments and a comparative analysis of the data obtained. The main conclusions of the study are the advantage of using the approximation of the fractional Caputo derivative compared to using the difference scheme for the fractional Riemann-Liouville derivative in the numerical solution of the anomalous diffusion equation. The paper also indicates the importance of choosing the method of difference approximation of the second derivative, which is a derivative of the Caputo.


Keywords:

numerical method, fractional derivative, stability, algorithm, difference scheme, diffusion, approximation, derivative Of Caputo, Riemann-Liouville derivative, differential equation


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

Процессы диффузионного переноса в пористых средах, жидких полимерах, стеклах часто отличаются от описываемых классическим уравнением диффузии. Такие процессы принято называть процессами аномальной диффузии. Аномальная диффузия, в отличие от классической, обладает пространственной нелокальностью и эффектами памяти. Математически это приводит к замене в уравнении диффузии обычных производных, производными дробного порядка [1]. Фактически это превращает классическое уравнение диффузии в частных производных в интегро-дифференциальные уравнения:

В научной литературе описано множество алгоритмов аппроксимации дробных производных для проведения численных расчетов. К классическим изданиям из данной серии можно отнести монографию китайских авторов Ч. Ли и Ф.Зенг [2]. При этом появляются новые работы, связанные с численными методами решения аномального уравнения диффузии, например, [3]. Развиваются и аналитические методы, например, [4–6]. При этом вопрос о целесообразности применения того или иного алгоритма для решения конкретной задачи представляется открытым. В настоящей работе проведено исследование влияния выбранного способа аппроксимации пространственной дробной производной в аномальном уравнении диффузии на устойчивость решения.

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

Введем не обязательно равномерную сетку на оси X: x1<x2…<xN. На узлах данной сетки представим искомую функцию в виде полинома первой степени. В этом случае интеграл, входящий в выражение (2), можно представить в виде суммы

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

где коэффициенты Dmk выражаются через Pmk и Qmk. Сама же вторая производная по координате при равномерной сетке заменяется конечной разностью:

В отличие от классического уравнения диффузии, устойчивые разностные схемы аппроксимации которого приводят к системам линейных уравнений с трехдиагональными матрицами, использование выражения (4) приводит к системе линейных уравнений с «почти треугольной» матрицей. Предполагается, что выражение (4) записывается на «верхнем» слое по времени при аппроксимации производной по времени.

В случае равномерной сетки по координате формула (3) соответствует первому порядку точности по h=xkxk-1. Для повышения точности аппроксимации неизвестной функции по координате можно воспользоваться интерполяционным полиномом Лагранжа более высокой степени n:

В результате получится выражение аналогичное (4), но приводящее к несколько большему объему вычислений.

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

В выражении (6) слева записана дробная производная Римана-Лиувилля, а справа — Капуто, которые для функций, обращающихся в ноль на бесконечности совпадают [7]. В работе [2] рассмотрены различные схемы аппроксимации дробной производной Капуто.

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

В этом случае для аппроксимации производной на отрезке [xk-1, xk] использованы значения функции в точках xk-1, xk и xk+1. Можно получить аналогичную формулу, когда в качестве базовых значений функции для вычисления производной используются точки xk-2, xk-1 и xk.. В этом случае получается система линейных уравнений с треугольной матрицей (за исключением одного элемента). В том случае, когда в уравнении диффузии используется производная Римана-Лиувилля, в работе [2] предлагается преобразовывать её к производной Капуто и далее применять разностные схемы к ней. Алгоритмы, основанные на использовании формул (3) и (4) не рассматривались.

В настоящей работе были проанализированы четыре алгоритма. Первый реализован на основе формулы (3). Второй является его аналогом, но для аппроксимации подынтегральной функции использована квадратичная интерполяция. Третий алгоритм основан на использовании формулы (7), а четвертый на основе её аналога, но с использованием второго набора базисных точек для аппроксимации производной.

Тестирование алгоритмов проводилось для расчета решения, соответствующего гауссовому начальному значению:

Точное решение уравнения (2) можно получить при помощи преобразования Фурье по координате. Используя свойства преобразования Фурье для дробных производных [1,7] можно записать уравнение для трансформанты Фурье:

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

В дальнейшем результаты вычислений по формуле (8) будем называть точным решением. Проведем численные расчеты с использованием описанных выше алгоритмов и сравним результаты с точным решением.

Для получения оценки устойчивости алгоритмов можно воспользоваться «принципом максимума» [8, c.360]. Представим исследуемую разностную схему в виде:

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

Для устойчивости по правой части необходимо дополнительно выполнение условия:

В случае уравнения диффузии с дробными производными наборы коэффициентов a и b будут разными для различных узловых точек m. Анализ, проведенный на основе формул (3) и (7) показывает, что существуют С>0 и g>0 при которых условия (10) и (11) выполняется для всех узловых точек первых трех рассматриваемых алгоритмов. В то же время не представилось возможным получить для произвольного µ оценку констант µ и С. Дальнейший анализ устойчивости алгоритмов проводился экспериментально путем проведения пробных расчетов.

Отметим сразу, что при использовании четвертого алгоритма не удалось получить устойчивых решений, хотя бы грубо напоминающих точное.

На рис. 1 приведены результаты расчетов для m=1,25 и m=1,95 по первому и второму алгоритмам в сравнении с точным решением. Приведенные результаты демонстрируют, что оба алгоритма позволяют достичь хорошей точности для заданных при расчетах значениях параметров. В тоже время отметим некоторые особенности алгоритмов. Для получения примерно одинаковой точности при расчете по первому алгоритму потребовалось задать шаг по координате в 5 раз меньше чем при расчете по второму алгоритму, где используется квадратичная аппроксимация неизвестной функции.

Рис. 1. Расчеты по алгоритму 1 и 2. Линии 2, 5 — точное решение для m=1,25 и m=1,95, соответственно. Пунктирная линия и х — первый алгоритм; о — второй алгоритм.

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

Из данного соотношения вытекает необходимость уменьшения шага по координате h при уменьшении шага по времени. В свою очередь величина шага по времени t определяет порядок точности аппроксимации производной по времени.

Хотя, как указывалось ранее, вычислительные затраты на реализацию одного шага по второму алгоритму несколько больше, чем по первому, более высокий порядок аппроксимации второго алгоритма компенсирует эти затраты. При этом для m ³ 1,3 при всех значениях hи t, позволяющих получить хорошую точность решения, второй алгоритм устойчив. Проблемы возникают при m близких к единице. При этом в отличие от первого алгоритма условие устойчивости второго алгоритма приближенно выглядит как

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

На рис. 2 приведены результаты расчетов для алгоритма, соответствующего аппроксимации дробной производной Капуто. Представленные результаты соответствуют a=1; h=0.4; t=0,02; t=15. Проводится сравнение с точными значениями для m=1.1 и m=1.9. Из данных рисунка видно, что точность численного расчета падает с уменьшением порядка дробной производной. Это связано с увеличением роли движения частиц без рассеяния, которое описывается первой производной по координате.

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

Рис. 2. Результаты расчетов по третьему алгоритму. Линии — точное решение; пунктир — численное решение для m=1.1, о — для m=1.9.

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

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

В заключение заметим, что существуют и приближенные методы решения уравнений с дробными производными. Одним из таковых является метод вариационного интерполирования, применение которого к уравнениям с дробными производными продемонстрировано в работах [5,6].

References
1. Uchaikin V.V. Metod drobnykh proizvodnykh. – Ul'yanovsk: Izd-vo «Artishok», 2008. – 512 s.
2. Changpin Li, Fanhai Zeng Numerical Methods for Fractional Calculus.– Taylor & Francis Group.– 2015. – 300 p.
3. Lukashchuk S.Yu. Dvukhsetochnye parallel'nye algoritmy dlya resheniya drobno-differentsial'nykh uravnenii anomal'noi diffuzii Vestnik Yuzhno-Ural'skogo gosudarstvennogo universiteta. Seriya: Vychislitel'naya matematika i informatika, 2012, №47(306) – S.83–98.
4. Ivashchenko D.S. Chislennye metody resheniya pryamykh i obratnykh zadach dlya uravneniya diffuzii drobnogo poryadka po vremeni: diss. kand. fiz.mat. nauk. – Tomsk: Institut matematiki SO RAN, 2008. – 187 s.
5. Litvinov V.A. Variatsionnoe interpolirovaniya reshenii drobnykh differentsial'nykh uravnenii Izvestiya vysshikh uchebnykh zavedenii. Fizika, 2018, t. 261, №8(728). – S.39–44.
6. Uchaikin V.V., Litvinov V.A. Nelineinaya teoriya vozmushchenii na osnove variatsionnogo printsipa: model'nye primery. Izvestiya vysshikh uchebnykh zavedenii. Prikladnaya nelineinaya dinamika, 2018, t. 26, № 6. – S. 82–98.
7. Samko S.G., Kilbas A.A., Marichev O.I. Integraly i proizvodnye drob-nogo poryadka i nekotorye ikh prilozheniya. – Minsk: Nauka i tekhnika, 1987. – 688 s.
8. Kalitkin N.N. Chislennye metody: ucheb. posobie. –2-e izd., ispravlennoe. – SPb.: BKhV-Peterburg, 2011. – 592 s.