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

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

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

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

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

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

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

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

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

В случае равномерной сетки по координате формула (3) соответствует первому порядку точности по h=xk– xk-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].