Журнал вычислительной математики и математической физики, 2019, T. 59, № 1, стр. 50-62

Аналитико-численный подход для описания периодических по времени движущихся фронтов в сингулярно возмущенных моделях реакция–диффузия–адвекция
В. Т. Волков, Д. В. Лукьяненко, Н. Н. Нефедов

В. Т. Волков^1, *, Д. В. Лукьяненко^1, **, Н. Н. Нефедов^1, ***

¹ МГУ, физ. ф-т
119992 Москва, Ленинские горы, Россия

^* E-mail: volkovvt@mail.ru
^** E-mail: lukyanenko@physics.msu.ru
^*** E-mail: nefedov@physics.ru

Поступила в редакцию 24.08.2018

DOI: 10.1134/S0044466919010150

Аннотация

В настоящей работе представлен аналитико-численный подход к изучению движущихся фронтов в сингулярно возмущенных моделях типа реакция–диффузия–адвекция. Предложен метод генерации динамически адаптированной сетки для эффективного численного решения задач данного класса. Метод основан на априорной информации о движении и свойствах фронта, полученной в результате строгого асимптотического анализа сингулярно возмущенной параболической задачи. В частности, существенными параметрами, которые учитываются при построении сетки, являются оценки местоположения переходного слоя, его ширина и структура. Предлагаемый аналитико-численный подход позволяет значительно сэкономить вычислительные ресурсы, сократить время счета и повысить стабильность работы вычислительного процесса по сравнению с классическими подходами. Рассмотрен пример, демонстрирующий основные идеи и методику применения предлагаемого подхода. Библ. 21. Фиг. 5.

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

1. ВВЕДЕНИЕ

Решения сингулярно возмущенных параболических задач часто содержат характерные узкие пограничные и внутренние слои (стационарные или движущиеся фронты.) Численное исследование таких задач с помощью разностных схем требует использования сеток с очень большим количеством узлов. Это приводит в некоторых случаях к значительным вычислительным затратам, а также к неустойчивости численных алгоритмов. Для преодоления обеих проблем мы предлагаем асимптотико-численный подход к решению задач с движущимися внутренними слоями в нелинейных моделях типа реакция–диффузия–адвекция. Основная идея подхода базируется на следующем: чем меньше параметр $\varepsilon $ в сингулярно возмущенной задаче, тем более грубым и неустойчивым получается его численное решение; в то же время, тем более точную априорную информацию о решении возможно получить с помощью асимптотического анализа. Соответствующее сочетание асимптотического подхода и вычислительных методов позволяет повысить эффективность численных расчетов.

Эта идея была использована в последнее время многими авторами для задач со стационарными пограничными и внутренними слоями, например [1]–[10], где были предложены методы построения специальных сеток, учитывающих особенности решения. В случае движущихся внутренних слоев ряд разностных схем описан в работах [11]–[15], а также в работах [16] и [17], где были рассмотрены примеры периодических задач для сингулярно возмущенных параболических уравнений.

В данной работе представлен эффективный аналитико-численный подход для исследования периодических по времени задач типа реакция–диффузия–адвекция, решения которых содержат движущиеся внутренние слои (движущиеся фронты). Предложен метод построения динамически адаптированной сетки, основанный на априорной информации о решении, полученной в результате строгого асимптотического анализа задачи, развитого в [18]. Указанный подход позволяет cэкономить вычислительные ресурсы и увеличить скорость построения решения с приемлемой точностью, что особенно важно, например, при решении обратных задач для сингулярно возмущенных уравнений. Оптимизация решения обратной задачи для уравнения реакция–диффузия–адвекция с использованием динамически адаптированной сетки реализована авторами в [19].

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

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

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

Статья имеет следующую структуру. В разд. 2 кратко описан асимптотический алгоритм, позволяющий получить априорную информацию, которая используется в дальнейшем для построения динамически адаптированной сетки. В разд. 3 кратко описаны основные идеи построения сетки. Наконец, в разд. 4 представлен численный эксперимент, обсуждаются его результаты и некоторые нюансы построения сетки в рассматриваемом примере.

2. АСИМПТОТИЧЕСКИЙ АНАЛИЗ

Основные идеи предлагаемого подхода продемонстрируем на примере следующей задачи:

(1)

$\begin{gathered} \varepsilon \left( {\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{2}}}} - \frac{{\partial u}}{{\partial t}}} \right) = A(u,x,t)\frac{{\partial u}}{{\partial x}} + B(u,x,t), \\ (x,t) \in D: = \left\{ {x \in ( - 1,1);\quad t \in \mathbb{R}} \right\}, \\ u( - 1,t) = {{u}_{{{\text{left}}}}}(t),\quad u(1,t) = {{u}_{{{\text{right}}}}}(t),\quad t \in \mathbb{R}, \\ u(x,0) = u(x,t + T),\quad x \in [ - 1,1],\quad t \in \mathbb{R}, \\ \end{gathered} $

где $\varepsilon $ – малый параметр $(0 < \varepsilon \ll 1)$, а функции $A(u,x,t)$, $B(u,x,t)$, ${{u}_{{{\text{left}}}}}(t)$ и ${{u}_{{{\text{right}}}}}(t)$ – достаточно гладкие и $T$-периодические по переменной $t$.

Известно (см. [18]), что при определенных условиях сформулированная задача допускает решение вида движущегося фронта: на интервале (–1, 1) существует точка ${{x}_{{tr}}}(t,\varepsilon )$, движущаяся по периодическому во времени закону, в окрестности которой наблюдается узкий внутренний переходный слой. А именно, слева от указанной точки (при $ - 1 < x < {{x}_{{tr}}}(t,\varepsilon )$) решение близко к некоторому уровню ${{\varphi }^{l}}(x,t)$, а справа (при ${{x}_{{tr}}}(t,\varepsilon ) < x < 1$) – к уровню ${{\varphi }^{r}}(x,t) \ne {{\varphi }^{l}}(x,t)$ для всех $t$. Основной целью настоящей работы является развитие метода построения динамически адаптированной сетки для эффективного численного решения задач подобного типа с использованием априорной информации о локализации, ширине и структуре движущегося фронта, полученной путем строгого асимптотического анализа задачи (1), выполненного в работе [18]. Ниже мы кратко опишем основные идеи получения асимптотического приближения решения задачи (1) и используем ряд результатов и формул из [18].

Следуя [18], сформулируем основные условия.

Положив $\varepsilon = 0$ в (1), получаем вырожденное уравнение

(2)

$A(u,x,t)\frac{{du}}{{dx}} + B(u,x,t) = 0,$

где $t$ играет роль параметра. Определим в области $D: = \left\{ {x \in ( - 1,1),\;t \in \mathbb{R}} \right\}$ две функции ${{\varphi }^{l}}(x,t)$ и ${{\varphi }^{r}}(x,t)$ как решения следующих задач:

(3)

$\begin{gathered} {{\varphi }^{l}}(x,t)\,{\text{:}}\;\;A(u,x,t)\frac{{du}}{{dx}} + B(u,x,t) = 0,\quad u( - 1,t) = {{u}_{{{\text{left}}}}}(t), \\ {{\varphi }^{r}}(x,t)\,:\;\;A(u,x,t)\frac{{du}}{{dx}} + B(u,x,t) = 0,\quad u(1,t) = {{u}_{{{\text{right}}}}}(t). \\ \end{gathered} $

Условие 1. Пусть для всех $(x,t) \in \bar {D}$ существуют $T$-периодические по переменной $t$ решения ${{\varphi }^{l}}(x,t)$ и ${{\varphi }^{r}}(x,t)$ задач (3), причем для всех $(x,t) \in \bar {D}$ выполнены неравенства

(4)

$\begin{gathered} {\text{а }})\quad {{\varphi }^{l}}(x,t) < {{\varphi }^{r}}(x,t), \hfill \\ {\text{б }})\quad A({{\varphi }^{l}}(x,t),x) > 0,\quad A({{\varphi }^{r}}(x,t),x) < 0. \hfill \\ \end{gathered} $

Определим следующую функцию в виде

(5)

$I(x,t): = \int\limits_{{{\varphi }^{l}}(x,t)}^{{{\varphi }^{r}}(x,t)} A(u,x,t)du.$

Условие 2. Пусть уравнение

(6)

$I(x,t) = 0$

имеет $T$-периодическое решение ${{x}_{0}}(t)$, причем для всех $t \in \mathbb{R}$ имеем

(7)

$\begin{gathered} {\text{a}})\quad - {\kern 1pt} 1 < {{x}_{0}}(t) < 1, \hfill \\ {\text{б }})\quad \int\limits_{{{\varphi }^{l}}({{x}_{0}}(t),t)}^s A(u,{{x}_{0}}(t),t)du > 0\quad {\text{п р и }}\quad s \in \left( {{{\varphi }^{l}}({{x}_{0}}(t),t),\;{{\varphi }^{r}}({{x}_{0}}(t),t)} \right). \hfill \\ \end{gathered} $

Условие 3. Пусть на решении ${{x}_{0}}(t)$ уравнения (6) выполнено неравенство

(8)

$\frac{{\partial I}}{{\partial x}}\left( {{{x}_{0}}(t),t} \right) < 0\quad {\text{п р и }}\;{\text{в с е х }}\quad t \in \mathbb{R}.$

Положение переходного слоя (фронта) будем описывать функцией ${{x}_{{tr}}}(t,\varepsilon )$, которую представим в виде ряда по параметру $\varepsilon $

(9)

${{x}_{{tr}}}(t,\varepsilon ) = {{x}_{0}}(t) + \varepsilon {{x}_{1}}(t) + {{\varepsilon }^{2}}{{x}_{2}}(t) + \ldots ,$

где ${{x}_{i}}(t)$, $i = 1,2,\; \ldots $, суть $T$-периодические функции, определяемые в процессе построения асимптотического приближения.

В дальнейшем использованы следующие обозначения:

$\begin{gathered} \mathop D\nolimits^l : = \left\{ {x \in [ - 1,{{x}_{{tr}}}(t,\varepsilon )],\;t \in \mathbb{R}} \right\}, \\ \mathop {\bar {D}}\nolimits^r : = \left\{ {x \in [{{x}_{{tr}}}(t,\varepsilon ),1],\;t \in \mathbb{R}} \right\}, \\ \end{gathered} $

$\begin{gathered} \varphi (x,t): = \frac{1}{2}\left( {{{\varphi }^{l}}(x,t) + {{\varphi }^{r}}(x,t)} \right), \\ \xi : = \frac{{x - {{x}_{{tr}}}(t,\varepsilon )}}{\varepsilon }. \\ \end{gathered} $

Асимптотика решения задачи (1) в областях $\mathop {\bar {D}}\nolimits^l $ и $\mathop {\bar {D}}\nolimits^r $ строится в виде

(10)

${{U}^{{l,r}}}(x,t,\varepsilon ) = {{\bar {u}}^{{l,r}}}(x,t,\varepsilon ) + {{Q}^{{l,r}}}(\xi ,t,\varepsilon ) = \sum\limits_{i = 0}^\infty \,{{\varepsilon }^{i}}\left( {{{{\bar {u}}}^{{l,r}}}(x,t) + Q_{i}^{{l,r}}(\xi ,t)} \right).$

Здесь ${{\bar {u}}^{{l,r}}}(x,t,\varepsilon )$ – регулярные функции, описывающие решение вдали от точки перехода ${{x}_{{tr}}}(t,\varepsilon )$ в областях ${{\bar {D}}^{l}}$ и ${{\bar {D}}^{r}}$; функции ${{Q}^{{l,r}}}(\xi ,t,\varepsilon )$ описывают переходный слой (движущийся фронт) и существенны вблизи точки ${{x}_{{tr}}}(t,\varepsilon )$; $\xi = (x - {{x}_{{tr}}}(t,\varepsilon )){\text{/}}\varepsilon $, причем $\xi \leqslant 0$ для функций с индексом $l$ и $\xi \geqslant 0$ для функций с индексом $r$.

Решения ${{U}^{l}}(x,t,\varepsilon )$ и ${{U}^{r}}(x,t,\varepsilon )$ в областях ${{\bar {D}}^{l}}$ и ${{\bar {D}}^{r}}$ гладко сшиваются в точке ${{x}_{{tr}}}(t,\varepsilon )$, и члены ряда (10) определяются путем применения асимптотической процедуры, описанной в [18].

Функции $\bar {u}_{0}^{l}(x,t) = {{\varphi }^{l}}(x,t)$ и $\bar {u}_{0}^{r}(x,t) = {{\varphi }^{r}}(x,t)$ суть решения вырожденного уравнения (2) при $x \in \left[ { - 1,{{x}_{{tr}}}(t,\varepsilon )} \right)$ и $x \in \left( {{{x}_{{tr}}}(t,\varepsilon ),1} \right]$ соответственно. Далее, $\bar {u}_{k}^{{l,r}}(x,t)$, $k \geqslant 1$, находятся из линейных уравнений [18] и при каждом $k$ могут быть получены в явном виде.

Главный член асимптотики внутреннего переходного слоя (движущегося фронта) определяется из задачи

(11)

$\begin{gathered} \frac{{{{\partial }^{2}}\tilde {Q}}}{{\partial {{\xi }^{2}}}} - \left( {A(\tilde {Q},{{x}_{{tr}}}(t,\varepsilon ),t)} \right)\frac{{\partial{ \tilde {Q}}}}{{\partial \xi }} = 0,\quad \xi \in ( - \infty , + \infty ),\quad t \in \mathbb{R}, \\ \tilde {Q}(0,{{x}_{{tr}}}(t,\varepsilon )) = \varphi ({{x}_{{tr}}}(t,\varepsilon ),t), \\ \tilde {Q}( - \infty ,{{x}_{{tr}}}(t,\varepsilon )) = {{\varphi }^{l}}({{x}_{{tr}}}(t,\varepsilon ),t), \\ \tilde {Q}( + \infty ,{{x}_{{tr}}}(t,\varepsilon )) = {{\varphi }^{r}}({{x}_{{tr}}}(t,\varepsilon ),t), \\ \end{gathered} $

где

(12)

$\tilde {Q}(\xi ,{{x}_{{tr}}}(t,\varepsilon )) = \left\{ \begin{gathered} {{\varphi }^{l}}({{x}_{{tr}}}(t,\varepsilon ),t) + Q_{0}^{l}(\xi ,t)\quad {\text{п р и }}\quad \xi \leqslant 0, \hfill \\ {{\varphi }^{r}}({{x}_{{tr}}}(t,\varepsilon ),t) + Q_{0}^{r}(\xi ,t)\quad {\text{п р и }}\quad \xi \geqslant 0. \hfill \\ \end{gathered} \right.$

Требование условия 1 обеспечивает существование единственного решения задачи (11). Кроме того, имеют место экспоненциальные оценки

(13)

$\left| {\tilde {Q}(\xi ,t) - {{\varphi }^{{l,r}}}({{x}_{{tr}}}(t,\varepsilon ),t)} \right| \leqslant C{{e}^{{ - \varkappa |\xi |}}}\quad {\text{п р и }}\quad \xi \to \pm \infty ,\quad t \in \mathbb{R},$

где $C$ и $\varkappa $ – положительные константы.

Функции $Q_{k}^{{l,r}}(\xi ,t)$, $k = 1,2,3,\; \ldots $, определяются из линейных задач, могут быть выписаны в явном виде и удовлетворяют экспоненциальным оценкам

(14)

$\left| {Q_{k}^{{l,r}}(\xi ,t)} \right| \leqslant C{{e}^{{ - \varkappa |\xi |}}}\quad {\text{п р и }}\quad \xi \to \pm \infty ,\quad t \in \mathbb{R},$

где $C$ и $\varkappa $ – положительные константы.

Уравнение для определения нулевого приближения ${{x}_{0}}(t)$ положения движущегося фронта имеет вид

(15)

$\int\limits_{{{\varphi }^{l}}({{x}_{0}}(t),t)}^{{{\varphi }^{r}}({{x}_{0}}(t),t)} A(u,{{x}_{0}}(t),t)du \equiv I({{x}_{0}}(t),t) = 0.$

Требование условия 2 гарантирует существование $T$-периодического решения ${{x}_{0}}(t) \in ( - 1;1)$ задачи (15).

Для ${{x}_{1}}(t)$ и ${{x}_{k}}(t)$, $k = 2,3,\; \ldots $ – следующих членов ряда (9), описывающего положение движущегося фронта, – получаются линейные уравнения

(16)

${{x}_{k}}(t){{I}_{x}}({{x}_{0}}(t),t) = {{\Phi }_{k}}(t),$

где ${{\Phi }_{k}}(t)$ – известные достаточно гладкие $T$-периодические функции. В частности, явное выражение для ${{\Phi }_{1}}(t)$ приведено в [18]. Условие 3 обеспечивает существование $T$-периодических решений ${{x}_{k}}(t)$ задач (16) для всех $k = 1,2,3,\; \ldots $ .

Основной результат асимптотического анализа задачи (1) содержится в следующей теореме из [18].

Теорема 1. Пусть выполнены условия 1–3. Тогда при достаточно малых $\varepsilon $ существует $T$-периодическое по переменной $t$ решение $u(x,t,\varepsilon )$ задачи (1) с внутренним переходным слоем, локализованным в окрестности точки ${{x}_{0}}(t)$, а именно

$\mathop {lim}\limits_{\varepsilon \to 0} u(x,t,\varepsilon ) = \left\{ \begin{gathered} {{\varphi }^{l}}(x,t)\quad п р и \quad x \in [ - 1;{{x}_{0}}(t)),\quad t \in \mathbb{R}, \hfill \\ {{\varphi }^{r}}(x,t)\quad п р и \quad x \in ({{x}_{0}}(t);1],\quad t \in \mathbb{R}, \hfill \\ \end{gathered} \right.$

причем

$\left| {u(x,t,\varepsilon ) - {{U}_{n}}(x,t,\varepsilon )} \right| \leqslant C{{\varepsilon }^{n}},\quad (x,t) \in \bar {D},$

где положительная постоянная C не зависит от ε, а

(17)

${{U}_{n}}(x,t,\varepsilon ) = \left\{ \begin{gathered} U_{n}^{l}(x,t,\varepsilon )\quad п р и \quad - 1 \leqslant x \leqslant \sum\limits_{i = 0}^{n + 1} {{\varepsilon }^{i}}{{x}_{i}}(t),\quad t \in \mathbb{R}, \hfill \\ U_{n}^{r}(x,t,\varepsilon )\quad п р и \quad \sum\limits_{i = 0}^{n + 1} {{\varepsilon }^{i}}{{x}_{i}}(t) \leqslant x \leqslant 1,\quad t \in \mathbb{R}. \hfill \\ \end{gathered} \right.$

В формуле (17) функции $U_{n}^{{l,r}}(x,t,\varepsilon )$ – частичные суммы ряда (10), где $\xi $ заменено на

${{\xi }_{{n + 1}}} = \frac{{x - \sum\limits_{i = 0}^{n + 1} {{\varepsilon }^{i}}{{x}_{i}}(t)}}{\varepsilon }.$

Подводя итог, сформулируем основные результаты асимптотического анализа задачи (1), которые будут использованы в дальнейшем при построении динамически адаптированной сетки для численного решения.

Положение и скорость движения фронта. Два главных члена асимптотики по $\varepsilon $ локализации фронта ${{x}_{0}}(t) + \varepsilon {{x}_{1}}(t)$ определены в (15) и (16) (скорость движения фронта есть производная этих функций).

Оценка ширины переходного слоя может быть получена из (13) и (14). Ввиду экспоненциального характера примыкания функций, описывающих переходный слой, к уровням ${{\varphi }^{{r,l}}}(x,t)$ при $\xi \to \pm \infty $, в качестве оценки ширины слоя можно использовать

(18)

$d = C\left| {\varepsilon ln\varepsilon } \right|.$

Ширина области, где должно быть произведено сгущение сетки, зависит также от параметра $\varkappa $ из (13). При заданных функциях $A(u,x,t)$ и $B(u,x,t)$ и граничных условиях ${{u}_{{{\text{left}}}}}(t)$, ${{u}_{{{\text{right}}}}}(t)$) этот параметр может быть оценен. Для примера, который рассмотрен ниже в разд. 4 статьи, можно получить, что $\varkappa > 1$.

Структура внутреннего переходного слоя определяется видом решения задачи (11) и оценками (13), (14).

В разд. 4 предлагаемый подход построения динамически адаптированной сетки проиллюстрирован на конкретном примере задачи типа (1) с $A(u,x,t) = - u$ и $B(u,x,t) = u \times b(t)$, где $b(t) = 2 + cos(4\pi t)$:

$\begin{gathered} \varepsilon \left( {\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{2}}}} - \frac{{\partial u}}{{\partial t}}} \right) = - u\frac{{\partial u}}{{\partial x}} + ub(t),\quad x \in ( - 1;\;1),\quad t \in \mathbb{R}, \\ u( - 1,t) = {{u}_{{{\text{left}}}}}(t) \equiv - 8 + sin(4\pi t), \\ u(1,t) = {{u}_{{{\text{right}}}}}(t) \equiv 8 - 2sin(4\pi t). \\ \end{gathered} $

Для данного примера из (3) получим

(19)

${{\varphi }^{l}}(x,t) = - 8 + sin(4\pi t) + (x + 1)\left( {2 + cos(4\pi t)} \right),$

(20)

${{\varphi }^{r}}(x,t) = 8 - 2sin(4\pi t) + (x - 1)\left( {2 + cos(4\pi t)} \right).$

Условие 1 выполнено, так как для всех $x \in [ - 1,\;1]$ имеет место

a) ${{\varphi }^{l}}(x,t) - {{\varphi }^{r}}(x,t) = - 16 + 2\left( {2 + cos(4\pi t)} \right) + 3sin(4\pi t) < 0;$

б) $A({{\varphi }^{l}}(x,t),x,t)) = 8 - (x + 1)\left( {2 + cos(4\pi t)} \right) - sin(4\pi t) > 0,$

$A({{\varphi }^{r}}(x,t),x,t) = - 8 - (x - 1)\left( {2 + cos(4\pi t)} \right) + 2sin(4\pi t) < 0.$

Функция $I(x,t)$, определенная в (5), имеет вид

$\begin{gathered} I(x,t) = \int\limits_{{{\varphi }^{l}}(x,t)}^{{{\varphi }^{r}}(x,t)} - {\kern 1pt} udu = \frac{1}{2}\left( {{{\varphi }^{l}}(x,t) - {{\varphi }^{r}}(x)} \right)\left( {{{\varphi }^{l}}(x,t) + {{\varphi }^{r}}(x)} \right) = \\ = \;\frac{1}{2}\left( {2x\left( {2 + cos(4\pi t)} \right) - sin(4\pi t)} \right)\left( { - 16 + 2\left( {2 + cos(4\pi t)} \right) + 3sin(4\pi t)} \right), \\ \end{gathered} $

и решение уравнения (15), определяющего главный член асимптотики положения фронта, есть

(21)

${{x}_{0}}(t) = \frac{{sin(4\pi t)}}{{4 + 2cos(4\pi t)}} \in [ - 1;1]\quad {\text{п р и }}\quad t \in \mathbb{R}.$

Нетрудно проверить, что условия 2 и 3 для ${{x}_{0}}(t)$ выполнены.

В следующем приближении из уравнения (16) с учетом явного выражения для ${{\Phi }_{1}}(t)$ из [18] получаем

(22)

${{x}_{1}}(t) = - \frac{1}{{2b(t)}}\left[ {2x_{0}^{'}(t) + \left( {\bar {u}_{1}^{l}({{x}_{0}}(t),t) + \bar {u}_{1}^{r}({{x}_{0}}(t),t)} \right)} \right],$

где

$\begin{gathered} \bar {u}_{1}^{l}(x,t) = \frac{d}{{dt}}\left[ {\frac{{{{u}_{{{\text{left}}}}}(t)}}{{b(t)}}} \right]ln\left( {1 + \frac{{b(t)}}{{{{u}_{{{\text{left}}}}}(t)}}(x + 1)} \right) + \frac{d}{{dt}}\left[ {ln(b(t))} \right](x + 1) \hfill \\ \hfill \\ \end{gathered} $

$\begin{gathered} \bar {u}_{1}^{r}(x,t) = \frac{d}{{dt}}\left[ {\frac{{{{u}_{{{\text{right}}}}}(t)}}{{b(t)}}} \right]ln\left( {1 + \frac{{b(t)}}{{{{u}_{{right}}}(t)}}(x - 1)} \right) + \frac{d}{{dt}}\left[ {ln(b(t))} \right](x - 1). \hfill \\ \hfill \\ \end{gathered} $

3. ПОСТРОЕНИЕ ДИНАМИЧЕСКИ АДАПТИРОВАННОЙ СЕТКИ

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

Сформулируем основные этапы предлагаемого подхода. Вначале строится равномерная сетка ${{T}_{M}}$ только по переменной $t$ с числом узлов $M + 1$ (что эквивалентно $M$ интервалам сетки): ${{T}_{M}} = \left\{ {{{t}_{m}},\,0 \leqslant m \leqslant M{\text{:}}\;{{t}_{m}} = 0 + \tfrac{{T - 0}}{M}(m - 1)} \right\}$. Также возможно использование квазиравномерных сеток, описанных в [21], без каких-либо изменений в дальнейшем алгоритме. Заметим, что выбор правильного шага по времени $\tau = (T - 0){\text{/}}M$ в случае равномерной сетки очень важен. В процессе построения динамически адаптированной сетки возможно возникновение ситуации, когда область локального сгущения сетки на очередном временном слое не имеет пересечения (или данное пересечение слишком мало) с областью локального сгущения на предыдущем слое по времени. Во избежание указанной проблемы при переходе на новый слой нужно проверять условие на максимальный шаг по времени $\tau $, а именно, фронт не должен покидать текущую область локального сгущения сетки в пределах одного шага по времени:

(23)

$\tau \leqslant \mathop {\max }\limits_{t \in [0,T]} \left| {x_{{tr}}^{'}({{t}_{i}})} \right| \leqslant C\left| {\varepsilon ln\varepsilon } \right|.$

Таким образом, $\tau = O(\left| {\varepsilon ln\varepsilon } \right|)$ или $M = O({{\left| {\varepsilon ln\varepsilon } \right|}^{{ - 1}}})$, что означает, что сетка будет являться $\varepsilon $-зависимой по пространственной переменной.

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

Фиг. 1.

Динамически адаптированная сетка.

Далее, уточним некоторые детали процесса построения динамически адаптированной сетки ${{X}_{N}}{{T}_{m}} \equiv {{X}_{N}}({{t}_{m}})$, $0 \leqslant m \leqslant M$, для слоя $m$ сетки ${{T}_{M}}$.

1. Вначале вводится базовая равномерная сетка $X_{{{{N}_{0}}}}^{{(0)}}$ по переменной $x$ с шагом ${{h}_{0}} = (1 - ( - 1)){\text{/}}{{N}_{0}}$ и ${{N}_{0}} + 1$ узлами (т.е. ${{N}_{0}}$ интервалами): $X_{{{{N}_{0}}}}^{{(0)}} = \left\{ {x_{n}^{{(0)}},\,0 \leqslant n \leqslant {{N}_{0}}{\text{:}}\;x_{n}^{{(0)}} = 0 + n{{h}_{0}}} \right\}$. Число узлов ${{N}_{0}}$ выбирается на основе априорной информации о ширине переходного слоя и может быть получено с использованием оценки (18):

${{N}_{0}} = \left[ {\frac{{1 - ( - 1)}}{{\left| {\varepsilon ln\varepsilon } \right|}}} \right].$

2. Затем строится семейство кусочно-равномерных сеток $\left\{ {X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n,n + 1}}} \right\}$, $0 \leqslant n \leqslant {{N}_{0}} - 2$, где $(n + 1)$-й и $(n + 2)$-й интервалы базовой сетки дополнительно разбиты на ${{N}_{{{\text{int}}}}}$ интервалов. Таким образом, суммарное число интервалов теперь равно $N = {{N}_{0}} - 2 + 2{{N}_{{{\text{int}}}}}$:

$X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n,n + 1}} = \left\{ {\mathop {{{x}_{m}}}\limits_{}^{} ,\;0 \leqslant m \leqslant {{N}_{0}} - 2 + 2{{N}_{{{\text{int}}}}}} \right.,\;0 \leqslant n \leqslant {{N}_{0}}:0 = {{x}_{0}} = x_{0}^{{(0)}} < {{x}_{1}} = x_{1}^{{(0)}} < \ldots < {{x}_{n}} = $

$ = x_{n}^{{(0)}} < {{x}_{{n + 1}}} = x_{n}^{{(0)}} + 1\frac{{x_{{n + 1}}^{{(0)}} - x_{n}^{{(0)}}}}{{{{N}_{{{\text{int}}}}}}} < {{x}_{{n + 2}}} = x_{n}^{{(0)}} + 2\frac{{x_{{n + 1}}^{{(0)}} - x_{n}^{{(0)}}}}{{{{N}_{{{\text{int}}}}}}} < \ldots < {{x}_{{n + {{N}_{{{\text{int}}}}} - 1}}} = $

$ = x_{n}^{{(0)}} + ({{N}_{{{\text{int}}}}} - 1)\frac{{x_{{n + 1}}^{{(0)}} - x_{n}^{{(0)}}}}{{{{N}_{{{\text{int}}}}}}} < < {{x}_{{n + {{N}_{{{\text{int}}}}}}}} = x_{{n + 1}}^{{(0)}} < {{x}_{{n + {{N}_{{{\text{int}}}}} + 1}}} = x_{{n + 1}}^{{(0)}} + 1\frac{{x_{{n + 2}}^{{(0)}} - x_{{n + 1}}^{{(0)}}}}{{{{N}_{{{\text{int}}}}}}} < \ldots $

$ \ldots \; < {{x}_{{n + {{N}_{{{\text{int}}}}} + {{N}_{{{\text{int}}}}} - 1}}} = x_{{n + 1}}^{{(0)}} + ({{N}_{{{\text{int}}}}} - 1)\frac{{x_{{n + 2}}^{{(0)}} - x_{{n + 1}}^{{(0)}}}}{{{{N}_{{{\text{int}}}}}}} < {{x}_{{n + {{N}_{{{\text{int}}}}} + {{N}_{{{\text{int}}}}}}}} = x_{{n + 2}}^{{(0)}} < $

$ < \;{{x}_{{n + 2{{N}_{{{\text{int}}}}} + 1}}} = x_{{n + 3}}^{{(0)}} < \ldots < {{x}_{{{{N}_{0}} - 3 + 2{{N}_{{{\text{int}}}}}}}} = x_{{{{N}_{0}} - 1}}^{{(0)}} < \left. {{{x}_{{{{N}_{0}} - 2 + 2{{N}_{{{\text{int}}}}}}}} = x_{{{{N}_{0}}}}^{{(0)}} = \mathop 1\limits_{}^{} } \right\}.$

Заметим, что для выполнения апостериорных асимптотически точных оценок погрешностей удобно, чтобы базовая сетка $X_{{{{N}_{0}}}}^{{(0)}}$ была равномерной. В этом случае для оценки погрешности достаточно использовать только те узлы базовой сетки, которые совпадают на каждой сетке семейства $\left\{ {X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n,n + 1}}} \right\}$.

Фиг. 2.

Построение динамически адаптированной сетки: $\square $ – узлы, используемые при вычислениях; o – узлы, в которых производится интерполяция; $ \times $ – узлы, исключенные из рассмотрения на последующих шагах.

3. Используя формулы (21) и (22) (или, в общем случае, уравнения (15) и (16)), получаем оценку положения точки перехода

(24)

${{x}_{{tr}}}(t,\varepsilon ) = {{x}_{0}}(t) + \varepsilon {{x}_{1}}(t) + O({{\varepsilon }^{2}}).$

Важно, что в соответствии с результатом теоремы 1, для аппроксимации решения с точностью $O(\varepsilon )$, необходимо использовать два первых слагаемых ряда (22) для ${{x}_{{tr}}}(t,\varepsilon )$, описывающего положение движущегося фронта. В противном случае положение точки перехода может быть потеряно в процессе численных расчетов.

Имея оценку локализации фронта, фиксируем интервал базовой сетки, который содержит точку ${{x}_{{tr}}}({{t}_{0}},\varepsilon )$. Не ограничивая общности, можно считать, что ${{x}_{{tr}}}({{t}_{0}},\varepsilon )$ находится на $n$-м интервале $\left[ {x_{{n - 1}}^{{(0)}},x_{n}^{{(0)}}} \right)$ сетки $X_{{{{N}_{0}}}}^{{(0)}}$. Если ${{x}_{{tr}}}({{t}_{0}},\varepsilon ) \geqslant \left( {x_{{n - 1}}^{{(0)}} + x_{n}^{{(0)}}} \right){\text{/}}2$, то используем сетку ${{X}_{N}}({{t}_{0}}) = X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n - 1,n}}$ и полагаем $n: = n + 1$; в противном случае ${{X}_{N}}({{t}_{0}}) = X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n - 2,n - 1}}$ и $n: = n$.Считаем, что $m = 0$.

4. Далее, положим $m: = m + 1$. Если $\left( {x_{{n - 2}}^{{(0)}} + x_{{n - 1}}^{{(0)}}} \right){\text{/}}2 \leqslant {{x}_{{tr}}}({{t}_{m}},\varepsilon ) < \left( {x_{{n - 1}}^{{(0)}} + x_{n}^{{(0)}}} \right){\text{/}}2$, то используем сетку ${{X}_{N}}({{t}_{m}}) = {{X}_{N}}({{t}_{{m - 1}}})$. Если ${{x}_{{tr}}}({{t}_{m}},\varepsilon ) > (x_{{n - 1}}^{{(0)}} + x_{n}^{{(0)}})/2$, то для дальнейших вычислений используем сетку ${{X}_{N}}({{t}_{m}}) = X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n - 1,n}}$. В этом случае мы исключаем из процесса значения сеточной функции ${\mathbf{u}}({{t}_{{m - 1}}})$ на $(n - 1)$-м интервале базовой сетки $X_{{{{N}_{0}}}}^{{(0)}}$ и интерполируем ее на сгущенный $(n + 1)$-й интервал базовой сетки. Для правильной интерполяции используется априорная информация (см. (13)) об экспоненциальном характере стремления интерполируемой функции к уровням ${{\varphi }^{r}}(x,t)$ и ${{\varphi }^{l}}(x,t)$. В случае ${{x}_{{tr}}}({{t}_{m}},\varepsilon ) < \left( {x_{{n - 2}}^{{(0)}} + x_{{n - 1}}^{{(0)}})} \right){\text{/}}2$ для дальнейших вычислений используется сетка ${{X}_{N}}({{t}_{m}}) = X_{{{{N}_{0}},{{N}_{{{\text{int}}}}}}}^{{n - 3,n - 2}}$. Из процесса исключаются значения сеточной функции ${\mathbf{u}}({{t}_{{m - 1}}})$ на $n$-м интервале базовой сетки $X_{{{{N}_{0}}}}^{{(0)}}$ и аналогичным образом производится интерполяция на сгущенный $(n - 2)$-й базовый интервал. Затем полагаем $n: = n + 1$.

Как было отмечено выше, для интерполяции используется априорная информация (13) об экспоненциальном характере стремления интерполируемой функции к уровням ${{\varphi }^{r}}(x,t)$ и ${{\varphi }^{l}}(x,t)$. Записав этот факт в виде формулы $\left| {{{\varphi }^{r}}(x) - y} \right| = a{{e}^{{bx}}}$ ($\left| {{{\varphi }^{l}}(x) - y} \right| = a{{e}^{{bx}}}$) (или в эквивалентном виде $loga + bx = log\left| {{{\varphi }^{r}}(x) - y} \right|$) для пары интерполируемых точек $(x,y)$, можно получить коэффициенты $a$ и $b$ интерполирующей функции $f(x) = {{\varphi }^{r}}(x) + a{{e}^{{bx}}}$ и $f(x) = {{\varphi }^{l}}(x) + a{{e}^{{bx}}}$.

5. Если $m = M$, то процесс завершается. В противном случае, переходим к п. 4.

3.1. Динамически адаптированная $\varepsilon $-независимая сетка

Заметим, что динамически адаптированная сетка, процесс построения которой описан выше, является $\varepsilon $-зависимой по переменной $x$, т.е число узлов сетки растет, если $\varepsilon $ стремится к нулю. Однако этот выбор сетки имеет ряд преимуществ для численных расчетов в практических задачах с фиксированным $\varepsilon $. Во-первых, используя асимптотическую информацию о положении фронта, мы уточняем только два интервала базовой сетки, что в любом случае обеспечивает более экономичные численные расчеты по сравнению с равномерной сеткой. Во-вторых, узлы базовой грубой сетки фиксированы, следовательно, базовая сетка такого типа позволяет упростить процедуру интерполяции для вновь вводимых узлов нестационарной сгущенной сетки: не нужно вычислять значения сеточной функции во всех узлах, а необходимо лишь выполнить интерполяцию в двух вновь вводимых сгущенных интервалах сетки. Это дает возможность также упростить процедуру апостериорного контроля точности с помощью экстраполяции Ричардсона, выполненной только на базовой сетке.

Теперь рассмотрим другой подход к построению динамически адаптированной сетки, $\varepsilon $-независимой по пространственной переменной $x$. Эта сетка в каждый момент времени имеет одинаковое (или пропорциональное) количество узлов внутри и вне области локализации внутреннего слоя, причем их число не зависит от малого параметра $\varepsilon $.

Для каждого временного слоя $m$ сетки ${{T}_{M}}$ вводим кусочно-равномерную сетку $X{{T}_{m}}$ по переменной $x$, с равномерным сгущением в $\varepsilon ln\varepsilon $-окрестности точки перехода ${{x}_{{tr}}}({{t}_{m}})$:

$\begin{gathered} X{{T}_{m}} = \left\{ {{{x}_{n}},\;0 \leqslant n \leqslant N + KN} \right.:{{x}_{n}} = 0\;{\text{д л я }}\;n = 0,\;{{x}_{n}} = {{x}_{{n - 1}}} + {{h}_{{{\text{left}}}}}\;{\text{д л я }}\;n = \overline {1,{{N}_{{{\text{left}}}}}} , \\ {{x}_{n}} = {{x}_{{n - 1}}} + {{h}_{{{\text{int}}}}}\;{\text{д л я }}\;n = \overline {{{N}_{{{\text{left}}}}} + 1,{{N}_{{{\text{left}}}}} + KN} , \\ {{x}_{n}} = {{x}_{{n - 1}}} + {{h}_{{{\text{right}}}}}\;{\text{д л я }}\;\left. {n = \overline {{{N}_{{{\text{left}}}}} + {{K}_{{{\text{int}}}}}N + 1,N + KN} } \right\}, \\ \end{gathered} $

где $N$ – управляющий параметр, который определяет число интервалов внутри области переходного слоя. При этом интервалы вне области фронта нормируются так, что ${{h}_{{{\text{left}}}}} \cong {{h}_{{{\text{right}}}}}$:

$\begin{gathered} {{N}_{{{\text{left}}}}} = \left[ {\frac{{\left( {{{x}_{{tr}}}({{t}_{m}}) - C\left| {ln\varepsilon } \right|} \right) - ( - 1)}}{{1 - ( - 1) - 2C\left| {\varepsilon ln\varepsilon } \right|}}N} \right],\quad {{N}_{{{\text{right}}}}} = N - {{N}_{{{\text{left}}}}}, \\ {{h}_{{{\text{left}}}}} = \frac{{\left( {{{x}_{{tr}}}({{t}_{m}}) - C\left| {\varepsilon ln\varepsilon } \right|} \right) - ( - 1)}}{{{{N}_{{{\text{left}}}}}}},\quad {{h}_{{{\text{right}}}}} = \frac{{1 - \left( {{{x}_{{tr}}}({{t}_{m}}) + C\left| {\varepsilon ln\varepsilon } \right|} \right)}}{{{{N}_{{{\text{right}}}}}}},\quad {{h}_{{{\text{int}}}}} = \frac{{2C\left| {\varepsilon ln\varepsilon } \right|}}{{{{K}_{{{\text{int}}}}}N}}. \\ \end{gathered} $

Здесь $C$ – управляющий параметр, позволяющий регулировать ширину области сгущения (по умолчанию $C = 1$); $K$ – управляющий параметр, задающий относительную плотность сгущения в области внутреннего переходного слоя (по умолчанию $K = 1$).

В результате получаем динамически адаптированную сетку $X{{T}_{M}}$ (см. фиг. 3), которая на каждом слое по времени ${{T}_{m}}$ имеет фиксированное число интервалов $N + KN$.

Фиг. 3.

Пример построения $\varepsilon $-независимой динамически адаптированной сетки $X{{T}_{M}}$.

4. ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ

В данном разделе предложенный метод проиллюстрирован на примере уравнения Бюргерса с периодическими коэффициентами

(25)

$\begin{gathered} \varepsilon \left( {\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{2}}}} - \frac{{\partial u}}{{\partial t}}} \right) + u\frac{{\partial u}}{{\partial x}} - b(t)u = 0,\quad x \in ( - 1,\;1),\quad t \in (0,\;1.5], \\ u( - 1,t) = {{u}_{{{\text{left}}}}}(t),\quad u(1,t) = {{u}_{{{\text{right}}}}}(t),\quad t \in (0,\;1.5], \\ u(x,0) = {{u}_{{{\text{init}}}}}(x),\quad x \in [ - 1,\;1], \\ \end{gathered} $

где

(26)

$\begin{gathered} {{u}_{{{\text{left}}}}}(t) = - 8 + sin(4\pi t), \\ {{u}_{{{\text{right}}}}}(t) = 8 - 2sin(4\pi t), \\ b(t) = 2 + cos(4\pi t). \\ \end{gathered} $

Как было показано в конце разд. 2, все условия применимости процедуры асимптотического анализа [18] для данного примера выполнены.

Зададим начальное условие ${{u}_{{{\text{init}}}}}(x)$ как решение стационарной задачи (см. фиг. 4):

$\begin{gathered} \varepsilon \frac{{{{\partial }^{2}}g}}{{\partial {{x}^{2}}}} + g\frac{{\partial g}}{{\partial x}} - b(0)g = 0,\quad x \in ( - 1,\;1), \\ g( - 1) = {{u}_{{{\text{left}}}}}(0),\quad g(1) = {{u}_{{{\text{right}}}}}(0). \\ \end{gathered} $

Фиг. 4.

Пример начального условия ${{u}_{{{\text{init}}}}}(x)$ для $\varepsilon = {{10}^{{ - 2}}}$ (в окрестности переходного слоя выполнено сгущение сетки).

Для численного решения уравнения (25) применяется жесткий метод прямых (SMOL), т.е. уравнение в частных производных сводится к системе обыкновенных дифференциальных уравнений, которая может быть решена по схеме Розенброка с комплексными коэффициентами [20]:

$\begin{gathered} \frac{{d{{u}_{n}}}}{{dt}} = \frac{2}{{{{x}_{{n + 1}}} - {{x}_{{n - 1}}}}}\left( {\frac{{{{u}_{{n + 1}}} - {{u}_{n}}}}{{{{x}_{{n + 1}}} - {{x}_{n}}}} - \frac{{{{u}_{n}} - {{u}_{{n - 1}}}}}{{{{x}_{n}} - {{x}_{{n - 1}}}}}} \right) + \frac{1}{\varepsilon }{{u}_{n}}\frac{{{{u}_{{n + 1}}} - {{u}_{{n - 1}}}}}{{{{x}_{{n + 1}}} - {{x}_{{n - 1}}}}} - \frac{1}{\varepsilon }b(t){{u}_{n}},\quad n = \overline {1,N - 1} , \\ {{u}_{0}} = {{u}_{{{\text{left}}}}}(t),\quad {{u}_{N}} = {{u}_{{{\text{right}}}}}(t),\quad t \in (0,\;1.5], \\ u({{x}_{n}},0) = {{u}_{{{\text{init}}}}}({{x}_{n}}),\quad n = \overline {0,N} . \\ \end{gathered} $

Данная система может быть переписана в виде

$\begin{gathered} \frac{{d{\mathbf{u}}}}{{dt}} = {\mathbf{f}}({\mathbf{u}},t),\quad t \in (0,\;1.5], \\ {\mathbf{u}}(0) = {{{\mathbf{u}}}_{{{\text{init}}}}}, \\ \end{gathered} $

где ${\mathbf{u}} = {{({{u}_{1}},{{u}_{2}},{{u}_{3}},\; \ldots \;{{u}_{{N - 1}}})}^{{\text{т }}}}$, ${\mathbf{f}} = {{({{f}_{1}},{{f}_{2}},{{f}_{3}},\; \ldots \;{{f}_{{N - 1}}})}^{{\text{т }}}}$ and ${{{\mathbf{u}}}_{{{\text{init}}}}} = {{({{u}_{{{\text{init}}}}}({{x}_{1}}),{{u}_{{{\text{init}}}}}({{x}_{2}}),{{u}_{{{\text{init}}}}}({{x}_{3}}),\; \ldots \;{{u}_{{{\text{init}}}}}({{x}_{{N - 1}}}))}^{{\text{т }}}}$.

Численное решение системы (27) строится по схеме Розенброка с комплексным коэффициентом (CROS1), которая является монотонной, устойчивой и обеспечивает порядок точности $O({{\tau }^{2}})$ [20]:

(28)

${\mathbf{u}}({{t}_{{m + 1}}}) = {\mathbf{u}}({{t}_{m}}) + ({{t}_{{m + 1}}} - {{t}_{m}})Re\,{\mathbf{w}},$

где ${\mathbf{w}}$ – решение системы линейных алгебраических уравнений

$\left[ {E - \frac{{1 + i}}{2}\left( {{{t}_{{m + 1}}} - {{t}_{m}}} \right){{{\mathbf{f}}}_{{{\mathbf{u}}({\mathbf{u}}({{t}_{m}}),t)}}}} \right]{\mathbf{w}} = {\mathbf{f}}\left( {{\mathbf{u}}({{t}_{m}}),\frac{{{{t}_{m}} + {{t}_{{m + 1}}}}}{2}} \right),$

$E$ – единичная матрица, а ${{{\mathbf{f}}}_{{\mathbf{u}}}}$ – матрица Якоби.

В численной схеме (28) используются значения сеточной функции ${\mathbf{u}}({{t}_{m}})$, вычисленные на сетке $X{{T}_{m}}$. После применения схемы (28) необходимо интерполировать сеточную функцию ${\mathbf{u}}({{t}_{{m + 1}}})$ на сетку $X{{T}_{{m + 1}}}$. Заметим, что процедуры указанной интерполяции различны для разных подходов построения динамически адаптированной сетки, предложенных выше: в первом случае ($\varepsilon $-зависимая сетка) выполняется только локальная интерполяция с использованием информации о структуре внутреннего переходного слоя; во втором случае ($\varepsilon $-независимая сетка) необходимо производить полную интерполяцию.

Пример результатов численных расчетов представлен на фиг. 5.

Фиг. 5.

Пример численного расчета для $\varepsilon = {{10}^{{ - 2}}}$. Параметр ${{N}_{0}} = 43$ (вычисляется автоматически), ${{N}_{{{\text{int}}}}} = 100$ (управляющий параметр, выбираемый вручную).

5. ЗАКЛЮЧЕНИЕ

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

На основе этих идей предложен аналитико-численный метод для решения некоторых классов сингулярно возмущенных уравнений типа реакция–диффузия–адвекция, позволяющий существенно повысить устойчивость и скорость численных расчетов по сравнению с классическими подходами. Разработанный алгоритм с использованием динамически адаптированных сеток был применен авторами и показал высокую эффективность для решения систем с малыми параметрами (см. [14]), а также оптимизации решения обратных задач для сингулярно возмущенных уравнений типа реакция–диффузия–адвекция (см. [19]).

Список литературы

Shishkin G. Grid approximation of a singularly perturbed quasilinear equation in the presence of a transition layer // Russian Acad. Sci. Dokl. Math. 1993. V. 47. № 1. P. 83–88.
O’Riordan E., Shishkin G. Singularly perturbed parabolic problems with non-smooth data // J. of Comput. and Appl. Math. 2004. 66. № 1. P. 233–245.
Franz S., Kopteva N. Green’s function estimates for a singularly perturbed convection-diffusion problem // J. Differential Equations. 2012. V. 252. № 2. P. 1521–1545.
Kopteva N. Numerical analysis of a 2d singularly perturbed semilinear reaction-diffusion problem // Lecture Notes in Computer Sci. 2009. V. 5434. P. 80–91.
Kopteva N., O’Riordan E. Shishkin meshes in the numerical solution of singularly perturbed differential equations // Internat. J. of Numerical Analysis and Modeling. 2010. V. 1. № 1. P. 1–18.
Farrell P., Hegarty A., Miller J., O’Riordan E., Shishkin G. Robust computational techniques for boundary layers // Chapman and Hall/CRC, Boca Raton, FL. 2000.
O’Riordan E., Quinn J. Numerical method for a nonlinear singularly perturbed interior layer problem // Lect. Notes in Comput. Sci. and Eng. 2011. V. 81. P. 187–195.
O’Riordan E., Quinn J. Parameter-uniform numerical method for some linear and nonlinear singularly perturbed convection-diffusion boundary turning point problems // BIT Numerical Math. 2011. V. 51. № 2. P. 317–337.
Kopteva N., Stynes M. Numerical analysis of a singularly perturbed nonlinear reaction-diffusion problem with multiple solutions // Applied Numerical Math. 2004. V. 51. P. 273–288.
Kopteva N., Stynes M. Stabilised approximation of interior-layer solutions of a singularly perturbed semilinear reaction-diffusion problem // Numerische Math. 2011. V. 119. № 2. P. 787–810.
Shishkin G.I., Shishkina L.P., Hemker P.W. A class of singularly perturbed convection-diffusion problems with a moving interior layer. An a Posteriori Adaptive Mesh Technique // Comput. Meth. Appl. Math. 2004. V. 4. № 1. P. 105–127.
Shishkin G. Grid approximation of a singularly perturbed parabolic equation on a composite domain in the case of a concentrated source on a moving interface // Comput. Math. and Math. Phys. 2003. V. 43. № 12. P. 1738–1755.
Shishkin G. Necessary conditions for $\varepsilon $-uniform convergence of finite difference schemes for parabolic equations with moving boundary layers // Comput. Math. and Math. Phys. 2007. V. 47. № 10. P. 1636–1655.
Melnikova A., Levashova N., Lukyanenko D. Front dynamics in an activator-inhibitor system of equations // Lect. Notes in Comput. Sci. 2017. V. 10187. P. 492–499.
Volkov V., Lukyanenko D., Nefedov N. Asymptotic-numerical method for the location and dynamics of internal layers in singular perturbed parabolic problems // Lect. Notes in Comput. Sci. 2017. V. 10187. P. 721–729.
Quinn J. A numerical method for a nonlinear singularly perturbed interior layer problem using an approximate layer location // Comput. and Appl. Math. 2015. V. 290. № 15. P. 500–515.
Lukyanenko D., Nefedov N., Nikulin E., Volkov V. Use of asymptotics for new dynamic adapted mesh construction for periodic solutions with an interior layer of reaction-diffusion-advection equations // Lect. Notes in Comput. Sci. 2017. V. 10187. P. 107–118.
Nefedov N., Recke L., Schneider K. Existence and asymptotic stability of periodic solutions with an interior layer of reaction-advection-diffusion equations // J. of Math. Analysis and Appl. 2013. V. 405. № 1. P. 90–103.
Lukyanenko D.V., Shishlenin M.A., Volkov V.T. Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data // Communicat. in Nonlinear Sci. and Numerical Simulat. 2018. V. 54. P. 233–247.
Al’shinn A.B., Al’shina E.A., Kalitkin N.N., Koryagina A.B. Rosenbrock schemes with complex coefficients for stiff and differential algebraic systems // Comput. Math. and Math. Phys. 2006. V. 46. № 8. P. 1320–1340.
Kalitkin N.N., Al’shin A.B., Al’shina E.A., Rogov B.V. Computations on Quasi-Uniform Grids. Moscow: Fizmatlit, 2005 [in Russian].

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал