ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОГО КЛАССА НЕЛИНЕЙНЫХ РЕАКЦИОННО-ДИФФУЗИОННЫХ УРАВНЕНИЙ ПАРАБОЛИЧЕСКОГО ТИПА С ЗАПАЗДЫВАНИЕМ

Математика.Прикладная математика.

Механика

УДК 519.63:519.688:004.94

А.В. Павельчук, Н.Л. Габрелян, А.Г. Масловская

ЧИСЛЕННОЕ РЕШЕНИЕ ОДНОГО КЛАССА НЕЛИНЕЙНЫХ РЕАКЦИОННО-ДИФФУЗИОННЫХ УРАВНЕНИЙ ПАРАБОЛИЧЕСКОГО ТИПА С ЗАПАЗДЫВАНИЕМ

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

NUMERICAL SOLUTION FOR ONE CLASS OF THE NONLINEAR PARABOLIC REACTION-DIFFUSION EQUATIONS WITH DELAY

The finite-difference computational scheme of one class of initial-boundary value problem for the one-dimensional nonlinear parabolic equation describing reaction-diffusion system with delay is presented in the paper. The algorithm is designed with use of the Crank-Nicolson numerical scheme. The program implementation of the problem is illustrated by a test example. The practical estimates of errors as -well as optimum values of computational parameters -were also specified.

Key -words: reaction-diffusion equation, delay parabolic partial differential equation, computational scheme, finite-difference method, Crank -Nicolson method.

Введение

Отличительной особенностью многих диффузионных процессов, протекающих в сложных физических и технических системах, является зависимость скорости изменения характерных величин не только от состояния объекта в текущий момент времени, но и от развития процесса в предшествующий период. Если характеристики исследуемой модели системы зависят от того, какие значения они принимали в определенный момент прошлого, то такие модели системы классифицируют как модели с запаздывающей обратной связью [1-3].

Наиболее известные прикладные задачи, для решения которых требуется учет эффекта запаздывания, описывают сложные реакционно-диффузионные явления при изучении демографических процессов, химических реакций, радиационных повреждений твердых тел, тепломассопере-носа, диффузии носителей заряда в полупроводниках и диэлектриках и др. [2-3]. В качестве математических операторов в таких моделях выступают дифференциальные операторы, приводящие математическую постановку к начально-граничным задачам для функционально-дифференциальных уравнений с частными производными параболического и гиперболического типов. При интерпретации результатов моделирования запаздывание, или наследственность, как правило, ассоциируется с конечной скоростью распространения возмущений или обусловлена инерционной природой самой системы, которая отвечает за формирование отклика на внешнее воздействие с некоторым временным лагом [4]. Так, к моделям с запаздыванием различного вида приводят инженерные приложения [5], задачи теории теплопроводности [6] и расчет характеристик зарядки диэлектриков [7].

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

Постановка задачи и конструирование вычислительной схемы

Рассмотрим начально-граничную задачу для одномерного нелинейного реакционно-диффузионного уравнения параболического типа в следующей обобщенной постановке:

= + / (" (*, г), и (х, I -1*)), (х, ()е(а, Ь )х(0, Т], (1)

и (а, () = у^) ,и (Ь, ^ = , (е(0 ,Т ], (2)

и (х, ^ = V (х, ^, х е[а, 6], t е[_г1\\о] > (3)

где и (х, t) - искомая функция; а - положительно определенный параметр; t - время запаздывания;

а, Ь - параметры, определяющие размер расчетной области; Т - параметр, определяющий временной диапазон.

Численное решение будем строить с использованием метода конечных разностей. Введем в рассмотрение конечно-разностную сетку

Пкт = * О, ,

Ь (Л Т

где 0.к {х;|х1. = /Л, 0 < / < ; { = кх,-ц < к <М}• ; к = ^ - шаг по координате; т = — шаг по времени; д =--индекс, соответствующий временному лагу.

В сеточной области О.1п непрерывной задаче (1)-(3) поставим в соответствие конечно-разностную, для этого вместо функций и (х,, и, (х,, ихх (х,^, /[и(х,{),и(х,t - {)) , V(х,, <р(7) ,

(// (/) рассмотрим их дискретные аналоги. Вычислительный алгоритм будем конструировать на основе явно-неявной схемы Кранка - Николсона, которая обладает абсолютной устойчивостью и имеет второй порядок аппроксимации по координате и второй - по времени [10]. Аппроксимация задачи (1)-(3) будет иметь вид:

ик+1 - и"

а ( к+1 Г. к+1 , к+1\\ , а ( к о к , к \\

(и,+1 - 2 и, + и, _1 ) + — (+1 - 2 щ + и, _1) +

т 2к2 у 1+1 & 7 2к

+/(ик, ик+1~?), 1 < / < N -1,1 < к < М -1,

ик = ук, икы =щк, 1 <к <М,

ик = ук, 0 </ <N , -ц < к < 0. (6)

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

«м + 2к 1-1

ui--- и,.,, = и,.

к+1 „ к+1

& + ^^ "2и& + ^] +1 (и&&&Ь

Для каждого нового временного слоя 4+1 и всех внутренних узлов расчетной области получим систему (N-1) линейных алгебраических уравнений с^-1 неизвестными. В первом и последнем уравнениях требуется учет граничных условий I рода согласно выражению (5). Для старта вычислительного процесса используются значения функции из предыстории процесса согласно (6). Общая структурная схема нумерации узлов расчетной области и направление выполняемых прогонок показаны на рис. 1.

Рис. 1. Схема нумерации узлов и направление прогонок.

После решения системы методом прогонки на каждом временном слое дополним вектор решения значениями функции на границах: ик0 = (рк при х = а и икы =щк при х = Ь . Следует заметить, что алгоритмической особенностью решения уравнения с запаздыванием является использование для расчета значений искомой функции на временном слое к +1 значений данной функции для слоя к +1 - q (как схематически показано на рис. 1). Это приводит к необходимости хранения значений

функции ик на q предыдущих временных слоях в динамическом массиве - матрице «окна запаздывания», которая перезаписывается на каждом временном слое.

Решение тестовой задачи

С целью демонстрации результатов программного решения класса задач, определенного в постановке (1)-(3), рассмотрим тестовый пример, для которого известно аналитическое решение. Сформулируем математическую постановку задачи в виде:

= 20 + и2 (х, г - 0.1) + ^ (х, г), (х, г) е (0,1) х (0,1],

и(0,х) = 0, и) = х, Г е(0,1], и (X,г1 ) = х2Х, X е [0,1], X е [—0.1,0],

(9) (10)

^(х,X) = х2 - 4Х - х4 (х - 0.1)2 - функция источника.

Точное решение начально-граничной задачи (8)-(10) может быть записано в форме: и (х, х ) = х2х. (11)

На рис. 2 представлены графики аналитического и приближенного решений задачи (8)-(10).

—#— приближенное решение 0 точное решение

Рис. 2. Сравнение графиков точного и приближенного решений для задачи (8)-(10), в момент времени ¿=1.

Для практической оценки отклонения решения и[к, полученного с помощью численного метода, от точного, вычисленного в узлах сеточной области и (xi, Хк ) , используем два вида норм - норму-максимум и евклидову норму:

и - и = тах\\и

\\\\и - и =

Е \\и(х>&ь)-&

I 0<1<N,0<к<М

Проведем исследование зависимости практической оценки погрешности решения при различных значениях характерных параметров вычислительного процесса - шагов по координате Л и времени т. На рис. 3 приведены графики зависимости ошибок как функций значений шага по координате (рассмотрен частный случай - шаг по времени выбирался равным шагу по координате).

Результаты показывают, что увеличение числа разбиений интервала расчетной области (0, 1) до значения Л=100 приводит к уменьшению абсолютной погрешности до значений 0.00062 (норма-максимум) и 0.00402 (евклидова норма). Дальнейшее же измельчение шага не приводит к улучшению точности получаемого решения - на графиках присутствует четко выраженный минимум. Поэтому в данных условиях задачи выбранные значения шага Л и числа разбиений Л^ представляются оптимальными.

°0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

Рис. 3. Графическая визуализация практической оценки погрешности численного решения.

Заключение

Таким образом, в работе сконструирована вычислительная схема решения первой краевой задачи для одномерного нелинейного реакционно-диффузионного уравнения параболического типа с запаздыванием, включенным в реакционное слагаемое. С целью экономии вычислительных затрат алгоритм построен на основе эффективной схемы Кранка - Николсона и процедуры пересчета матрицы «окна», содержащей значения искомой функции в период запаздывания. Программная реализация алгоритма для решения рассмотренного класса задач продемонстрирована на тестовом примере. Практическая оценка погрешности найденного решения позволила установить оптимальное значение числа разбиений расчетной области. Выполненная работа послужит основой для алгоритмизации и программной реализации обобщенной физико-математической модели процесса электронно-стимулированной зарядки полярных диэлектриков [11-12], модифицированной к случаю присутствия запаздывания по времени в реакционно-диффузионном уравнении модели.