В настоящее время многих исследователей интересует использование математических моделей с дробными производными в силу точности моделей при описании объектов в различных областях науки. Причинами возросшего интереса к дробным производным и дробным дифференциальным уравнениям являются следующие: 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]. Таким образом, наиболее целесообразной является разработка неявной разностной схемы для решения дробного дифференциального уравнения на основе известной схемы Кранка - Николсона. Неявная разностная схема для численного решения нестационарного одномерного дробного дифференциального уравнения Рассмотрим общее одномерное уравнение для моделирования миграции токсичных веществ в подземных водах с начальными и граничными условиями вида (1) начальные условия: c(x, 0) = c0, 0 < x < L; граничные условия: где c - концентрация токсичного вещества, мг/м3; c0- начальная концентрация токсичного вещества, мг/м3; D - коэффициент дисперсии в направлении оси х, м2/с; λ - константа скорости превращения токсичного вещества, 1/с; ν - скорость движения воды в порах почвы в направлении оси x, м/с; t - время, с; L - расстояние, на котором концентрация принимает значение с0, м; α - порядок производной. В 1967 г. М. Капуто ввел новое и полезное определение дробной производной, которое сегодня широко известно. Приведем это определение для дробной производной по времени порядка α [20, 21]: (2) где Г (.) - Гамма-функция. Разработаем неявную разностную схему, являющуюся аналогом схемы Кранка - Николсона, для численного решения дробного дифференциального уравнения. Схема Кранка - Николсона может быть записана в виде (3) В качестве сетки возьмем совокупность точек , координаты которых ; ; ; ; ; ; . На данной сетке будем строить разностную аппроксимацию дифференциального уравнения (1), используя шеститочечный шаблон: Ш . Следовательно, схема Кранка - Николсона для модели (1) и шеститочечного шаблона имеет вид (4) (5) (6) Приближенное решение задачи (1) получаем, заменяя значение исходной непрерывной функции на функцию, определяемую в узлах сетки (k - номер узла по времени и l - номер узла по пространству): и означает численную аппроксимацию точного решения (или обозначает численное приближение к точному решению). Аналог неявной разностной схемы Кранка - Николсона для численного решения дробного дифференциального уравнения (1) с учетом выражения для дробной производной Капуто по времени (2) можно сформулировать следующим образом: (7) Выполним последовательно преобразования правой части уравнения (7): (8) Запишем уравнение (1) с учетом (4)-(6), (8): (9) где ; . Проведем упрощающие преобразования. Примем, что коэффициент перед выражением в левой части равен единице, и введем следующие обозначения: (10) После преобразований (8)-(10) система разностных уравнений для первоначальной краевой задачи с дробной производной (1) имеет вид (11) Отметим, что при k = 0 слагаемое с суммой исчезает, что упрощает реализацию метода при проведении вычислений. Выполним преобразование уравнения (11): или (12) Коэффициенты в узлах сетки вычисляются по следующим формулам: Систему разностных уравнений (12) с учетом преобразований можно представить в следующем виде: или Первоначальная краевая задача с дробными производными (2) может быть представлена в матричном виде: для k = 0; (13) для , (14) где ; начальные условия и граничные условия: Матрицы А и В имеют следующий вид: ; 1 Решение задачи (13), (14) в конечном счете сводится к многократному решению системы линейных алгебраических уравнений (СЛАУ) с трехдиагональной матрицей. Для решения таких СЛАУ известны специальные методы, например метод прогонки и метод циклической редукции [9]. Наиболее простым и эффективным методом является метод прогонки, который привлекателен своей экономичностью и простотой реализации в последовательном режиме [8]. Именно поэтому метод прогонки использовался в данном случае для решения системы линейных алгебраических уравнений. Разработанный численный метод решения дробного дифференциального уравнения реализован в программном комплексе [22] для решения задачи распространения токсичных веществ в условиях аномальной диффузии. Алгоритм реализации численного метода представлен на рис. 1. Рис. 1. Алгоритм реализации разработанного численного метода Данный алгоритм предусматривает проверку корректности исходных и расчетных данных на всех этапах вычислений. Новый численный метод реализован на языке программирования Microsoft Visual Basic 6.5, в качестве системы управления базой данных использовался Microsoft Access. Оценка распространения токсичного вещества в подземных водах на основе численного метода В данном разделе приведены результаты численного эксперимента, иллюстрирующие вычислительные свойства нового метода. Выполнено сравнение концентраций вещества, полученных на основе решения нестационарного одномерного дробного дифференциального уравнения разработанным методом и на основании подхода, представленного в работе [11], для двух значений порядка производной по времени - α. Результаты приведены для модельного примера, в котором химические превращения токсичного вещества не рассматривались. Параметры расчетной сетки: h = 0,10; τ = 0,1. Число расчетных точек по осям - 100. На рис. 2 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при α = 0,2. Рис. 2. Зависимость концентрации токсичного вещества от расстояния (α = 0,2) На рис. 3 представлены графики зависимости концентрации вещества от расстояния, полученные с помощью разработанного численного метода и на основании аналитического решения, приведенного в работе [11], при α = 0,6. Рис. 3. Зависимость концентрации токсичного вещества от расстояния (α = 0,6) Из рис. 2 и 3 видно, что графики зависимостей концентрации токсичного вещества от расстояния изменяются симбатно. Результаты, полученные с помощью предлагаемого метода и известного точного решения дробного дифференциального уравнения, достаточно хорошо согласуются между собой. Относительная ошибка составляет в среднем 9 %. В отличие от известного аналитического решения разработанный численный метод может использоваться при моделировании распространения токсичных веществ в подземных водах с учетом их биодеградации. Заключение В результате проведенных исследований разработан новый численный метод решения дифференциального уравнения с производной дробного порядка по времени, а также неявная разностная схема, являющаяся аналогом схемы Кранка - Николсона. Система разностных уравнений представлена в матричном виде. Решение задачи сводится к многократному решению трехдиагональной системы линейных алгебраических уравнений методом прогонки. Представлены результаты оценки распространения токсичных веществ в подземных водах на основе численного метода при решении модельных задач. Полученные результаты позволяют рекомендовать разработанный метод численного решения нестационарного одномерного дробного дифференциального уравнения, основанный на использовании линейных интегро-дифференциальных уравнений достаточно общей структуры, для использования на практике при моделировании процессов тепло- и массопереноса в сильно неоднородных средах.