Неявная схема расщепления для некоторого класса дифференциальных уравнений параболического типа с возмущением Текст научной статьи по специальности «Математика»

УДК 519.6:51-73:004.942

А.В. Павельчук, А.Г. Масловская

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

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

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

IMPLICIT SPLITTING SCHEME FOR SOME CLASS OF PARABOLIC DIFFERENTIAL EQUATIONS WITH PERTURBATION

There was proposed acomputingscheme for implementation of the model described by parabolic differential equationwith perturbation subject to cylindrical symmetry of multidimensional problem. The algorithm was constructed using the implicitsplitting finite-difference scheme as well as thesweep method applied to solve the system of the finite-difference equation. The scheme permits mathematical modelling of complex processes like «reaction-diffusion» processes.

Key words: splitting scheme, parabolic differential equationwith perturbation, computational scheme, implicit finite difference scheme, sweep method.

Введение

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

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

ных соотношениях сеточных характеристик в задачах с граничными условиями такие схемы становятся условно устойчивыми [3-7].

Использование явных схем для аппроксимации параболических уравнений приводит к жесткому ограничению на шаг по времени для выполнения требования устойчивости и, следовательно, к очень большим временным затратам при вычислении схемы [5]. Применение неявных схем снимает ограничение, однако трудоемкость решения систем линейных алгебраических уравнений, возникающих при этом, может сделать невыгодным их использование для некоторого класса задач. Из экономичных конечно-разностных схем наибольшее распространение получили: метод переменных направлений, метод дробных шагов Яненко, схема центрально-симметричного метода [7].

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

Постановка задачи. Конечно-разностная аппроксимация задачи. Устойчивость схемы

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

в прямоугольнике: 0 < г < R, 0 < г < Z, 0 < t < Т , где с1, с2, с3, с4 - константы. Для замыкания постановки задачи (1) дополним ее начальным и (г, г,0) = 0 и смешанными краевыми условиями:

ди

-= 0 при г = 0 и 2 = 0,

Введем следующую пространственно-временную сетку с шагами И1, Н2, т соответственно по переменным г, г, t:

( А = (г. = /И1, / = 0,1,= ]к2, ] = 0, J, ^ = кт, к = 0,1,2,...} и на этой сетке будем аппроксимировать дифференциальную задачу (1)-(2) на верхнем временном слое ^+1 = (к + 1)т :

(1)

дп

и = 0 при г = Я и 2 = 2 .

(2)

+с

4

2<+1 + «*+1

+ -

Учитывая порядок аппроксимации разностных операторов, использованных при составлении разностной схемы, видим, что она имеет первый порядок аппроксимации по времени и второй - по каждой из координат 0(т, й12, й2).

Исследуем устойчивость неявной разностной схемы (3), аппроксимирующей дифференциальное уравнение (1) с помощью спектрального метода, описанного, например, в [5]. Представляем решение в виде гармоники:

и!.. = Аке'ше3, а е [0,2я],/ е [0,2п],

i, j

Xk+1eimelPj — Ake'melPj

= c -с2(ХУшeia])Xk+1eimeia] —

^k+\e'ai e'PJ — eia(i-i) eiPJ X^e™ e'^J — Xk+leiaiei^( J'1)

-c3 (-+-) +

h h2

Xk+ieia(.i+1)eiPJ — 2Xk+1eiaieiPj + xk+1eia(i—1) eiPJ

+c'(-ц-+

1 Xk+1eiaielPj — Xk+1eia(l—1)elPj

+--+

2 h

Xk+1eiaieip( J+1) — 2Xk+1eiaieiPj + Xk+1eiaieiP( J—1) +-h-

Для анализа схемы исключаем компонент Cj, который не влияет на устойчивость схемы.

Кроме того, метод замороженных коэффициентов позволяет проанализировать устойчивость схемы для уравнений с нелинейными коэффициентами [10]. Нелинейность определяется как с2и = (c2u)u, c2 = const, где выражение в скобках играет роль константы C в зависимости от r, z, t. После некоторых преобразований получим:

1

X = -

где Cj = Т(2С3? — C4); C2 = C3!; C3 =—CL- C4 =—CL- E = 1 — e^ ; E2 = 1 — e^ ; Pj = 4sin2 ^; 2h1

P2 = sin2 P . 2 2

1+Ct + CE + c2 E2 + C3 p + C4 P2

C3 — C4) . C = C3T • C = — C4T

2h 2 h2 ' 3 hj2 3 4 h22

Собственные числа оператора перехода удовлетворяют необходимому условию устойчивости: || < 1 при любых значениях т, к1,к2. Следовательно, неявная разностная схема (3) абсолютно устойчива.

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

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

Для построения неявной схемы будем использовать метод расщепления. Метод дробных шагов позволяет представить неявную разностную схему (3) в виде двух подсхем, каждая из которых может быть решена с помощью метода прогонки. Разобьем пополам интервал т между точками ^ и

г

^+1 на разностной сетке и обозначим полученную промежуточную точку как ^+12. Запишем на первом полушаге интервала х неявную разностную схему, которая будет учитывать только производную второго порядка по координате г . Получим первую подсхему:

к+1/2 к к+1/2 _ к+12 И , И ■ , , , И, ■ И, 1 ■

-^ = с, - с2ик и *+12 - с, -У--^ +

х ',3 ',3 к

(4)

ик+12 - 2ик+12 + ик+12 с ик+12 - ик+12 4 7 +с и.+1,- 2и ',- + и.-1,- + с± и',з и'-1,-

4 к 2 к ' Запишем на первом полушаге интервала X неявную разностную схему, которая будет учитывать только производную второго порядка по координате у . Получим вторую подсхему:

и к+1 - и к+12 и к+1 - и к+1 и к++1 - 2ик+1 + и к+\

- -сз-;-+ с4-Г2-. (5)

х 3 к2 4

Складывая подсхемы (4) и (5), получаем соотношения, отличающееся от неявной разностной схемы (3) тем, что вторая производная по координате г аппроксимируется в нем не на (к +1) -м шаге

по времени, а на шаге (к +1/2) . Таким образом, дифференциальное уравнение (1) может быть аппроксимировано с помощью последовательного разрешения двух подсхем (4) и (5), являющихся схемой расщепления.

Характеристика первой подсхемы

Приведем (4) к виду, удобному для использования метода прогонки:

ик+12 / и 2 И'+1,3 +

V

, схх 2схх с.х к

1 + + -2---+ сххик 3

к к 2к1

ик+12+

+

С4Х С3Х С4Х

2к к к

к+1/2 к . И'-1,з - И',з + с1х.

Введем обозначения:

с4т схх 2с.х с.х к с.х схх с.х к к а. -—V; Ь. -1 + — + —4---— + с2хик.; с , ---3---4—; ек. - и . + ах .

' к к к 2к 2 ',3 2^ к И ',3 ',3 1

Рекуррентное пограничное соотношение для первой подсхемы (4) имеет вид:

к+1/2 . и к+1/2 , к+12 к (г\

аИ,-+1,/з- + ь И,-,+ 1 + сгИн/з- - еи1. (6)

Прогоночные коэффициенты а,Д определяются из соотношения: ик++12 -аии+ Д следующим образом. Перепишем последнее соотношение в виде:

и -З2 - а-мк+12 + Д-1 . (7)

Подставим полученное соотношение в (6):

к+1/2 , 7 к+1/2 . к+12 . о к

аи,-+1,/з- + ьИ ,-,+' +са-ии+ + с,Р,- £,,з.

Выражаем ик

Ик+12 - аг Ик+12 + - С' Д-1

И ¿, з --1-И '+1, з +Т-.

Ьг + С'а,-1 Ьг + Са-1

Тогда пограничные коэффициенты имеют вид:

а. Д -^^^. (8)

ь + са ьг + са-1

к+12 : 1', 3 :

Выражения (8) позволяют рассчитать значения прогоночных коэффициентов на 1 -м шаге по координате г , если известны их значения на (1 -1) -м шаге по координате г . Следовательно, чтобы определить значения прогоночных коэффициентов на любом шаге по координате г, необходимо знать их значения на первом шаге, а1, Д . Для определения значений коэффициентов на первом шаге и решения на правой границе используются граничные условия по координате г .

Соотношения (7) и (8) включают переменную 3, поэтому необходимо задать внешний цикл по этой переменной: 3 = 2,..., — 1. Следовательно, при решении первой подсхемы (4) метод прогонки будет использован — 2 раза. Результатом решения первой подсхемы (4) схемы расщепления являются значения функции и на шаге по времени (k +1/2) , необходимые для решения второй подсхемы.

Характеристика второй подсхемы

Вторая подсхема (5) схемы расщепления абсолютно устойчива и решается с помощью метода прогонки.

Приведем (5) к виду, удобному для использования метода прогонки:

Сл Т

4 1,k+l

h2

u,,.,, +

( о Л

i+СзТ+2ст

h2 h2 j

wk++1 -

^ с3т + с4т ^

V h2 h2 j

2 i, j+i

2

Введем обозначения:

~ ст г с3т 2с т

а. =—4— ; b . = i+ — + —с. = -

j h2 j h2 h2 j

<+ii = <+i2.

с3т + с4т

V h2 h2 j

~k+i/2 k+i/2 ; sU+ = uu+ .

~ k+i . г k+i . ~ k+i ~k+i/2 /П\

au+ b.u.. + си... = s ' . (9)

j i,j+i j i,j j i,j-i i,j v '

Тогда рекуррентное пограничное соотношение для второй подсхемы (5) имеет вид:

к+1 ~ к+1 I. ■ = а и л

1,3 3 1,3 +1

<+i = j +ii + ß,, (i0)

а, _ ёк+12 — СД11 где а =— --3-, Д =-33 (11)

Ьз + Сзаз -1 Ь3 + Сзаз -1

- прогоночные коэффициенты.

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

Соотношение (10) и (11) включают переменную 1, поэтому необходимо задать внешний цикл по этой переменной 1 = 2,..., Ыг — 1. Следовательно при решении второй подсхемы (5) метод прогонки будет использован Ыг — 2 раза. Результатом второй подсхемы (5) схемы расщепления является значения функции и на (к +1) -м шаге по времени.

Заключение

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

1. Мартинсон, Л.К., Малов, Ю.И. Дифференциальные уравнения математической физики - М.: МГУ, i996. - i85 с.

2. Тихонов, А.Н., Самарский, А.А. Уравнения математической физики - М.: Наука, i977. - 736 с.

3. Формалев, В.Ф., Ревизников, Д.Л. Численные методы. - М.: ФИЗМАТЛИТ, 2004. - 400 с.

4. Самарский, А.А., Гулин, А.В. Устойчивость разностных схем. - М.: Наука, 1973. - 415 с.

5. Самарский, А.А. Теория разностных схем. - М.: Наука, 1989. - 479 с.

6. Марчук, Г.И. Методы расщепления для решения нестационарных задач // Журнал вычислительной математики и математической физики. - 1995. - Т. 35, № 6. - С. 843-849.

7. Holden, H., Karlsen, K.H., Lie, K.-A., Risebro, N.H. Splitting Methods for Partial Differential Equations with Rough Solutions. Analysis and MATLAB programs. - Zurich (Switzerland): European Mathematical Society Publishing House, 2010. - 226 p.

8. Maslovskaya, A., Pavelchuk, A. Simulation of dynamic charging processes in ferroelectrics irradiated with SEM // Ferroelectrics. - 2015. - V. 476. - P. 157-167.

9. Maslovskaya, A., Pavelchuk, A. Simulation of heat conductivity and charging processes in polar dielectrics induced by electron beam exposure // In: IOP Conf. Series: Materials Science and Engineering, 2015. - V. 81. - P. 012119 (doi: 10.1088/1757-899X/81/1/012119).

10. Годунов, С.К., Рябенький, В.С. Разностные схемы (введение в теорию) - М.: Наука, 1973. - 400 с.