В этом параграфе мы рассмотрим методы нахождения периодических решений в распределенных системах. Эти методы требуют обычно больших затрат машинного времени и большого объема памяти ЭВМ, и поэтому здесь трудно проводить достаточно полное параметрическое исследование системы. Большей частью речь здесь идет о численном (машинном) эксперименте. Сами методы будут описываться лишь в краткой форме. Во второй части параграфа мы коротко рассмотрим решения волнового характера, которые также могут появляться в системах, описываемых нелинейными параболическими уравнениями.

6.5.1. Вычисление периодических решений и процесс продолжения
Периодическое (во времени) решение системы типа «реакция-диффузия» (4.3.7) с заданными граничными условиями по определению удовлетворяет тождествам
\[
x(z, t+T) \equiv x(z, t), \quad y(z, t+T) \equiv y(z, t)
\]

при $z \in[0,1]$. Здесь $T$ – период данного решения.

Рис. 6.13. Периодическое решение задачи 11 при ГУ $1 ; A=2, B=5,45$, $D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004, L=0,75 ;$ период $T=3,12$.
Простейшим подходом для нахождения устойчивого периодического решения является использование какого-либо из методов динамического моделирования (см. § 6.4) в течение достаточно большого промежутка времени $t$, за который решение системы стабилизируется, становясь (приближенно) периодическим. Эта стабилизация может зависеть, конечно, от начальных условий. В качестве примера того, как может выглядеть периодическое решение в распределенной системе, на рис. 6.13 при-
ведено периодическое решение задачи 11 для случая ГУ1. На рисунке показано несколько типичных профилей переменной $x$ для нескольких моментов времени в течение одного периода.

В п. 6.3.3.2 мы определили точку комплексной бифуркации для задачи 14 (см. табл. 6.11). Динамическое моделирование уравнений (P14-7), (P14-8) для значений параметра $\mathrm{Da}$ из окрестности $\mathrm{Da}^{+} \approx 0,1254\left(\mathrm{Da}<\mathrm{Da}^{+}\right)$показало, что в данном случае существует ветвь периодических решений для значений $\mathrm{Da} \in(0,122 ; 0,1254)$. Значения конверсии и температур на входе в реактор и на выходе из него, т. е. величины $y(0, t), \Theta(0, t)$ и $y(1, t), \Theta(1, t)$, представлены на фазовой диаграмме $y-\Theta$ (рис. 6.14).

Рис. 6.14. Устойчивое периодическое решение задачи 14 (две проекции); $\gamma \rightarrow \infty, B=12, \beta=2, \Theta_{c}=0, \mathrm{Da}=0,123, \mathrm{Le}=1$.
Продолжение периодических решений по параметру для указанного случая можно проводить, применяя вышеописанный подход (т. е. численно моделируя процесс стабилизации) для некоторой последовательности значений данного параметра. Другая возможность состоит в использовании квазистационарного описания (см. §6.6).

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

6.5.1.1. Метод стрельбы в рамках метода прямых

Если для аппроксимации параболических уравнений воспользоваться методом прямых (п. 6.4.5), то условия принимают вид ( $t=0$ )
\[
x_{i}(T)=x_{i}(0), \quad y_{i}(T)=y_{i}(0), \quad i=1,2, \ldots, n-1 .
\]

Тем самым мы получаем нелинейную краевую задачу (см., например, (6.4.24)) для $2(n-1)$ дифференциальных уравнений первого порядка со смешанными граничными условиями (6.5.2).1) Ситуация здесь полностью аналогична той, которая имела место для задач, прямо приводящих к обыкновенным дифференциальным уравнениям (ср. (5.8.1), (5.8.2)), однако размерность изучаемой системы оказывается большой. Таким образом, мы можем использовать методику $\S 5.8$, которая основывалась на методе стрельбы, а также на идее алгоритма продолжения DERPER (см. п. 5.8.4).

Продемонстрируем использование метода прямых на примере задачи 11 с граничными условиями первого рода. Используя трехточечную замену по переменной $z$, получим
\[
\begin{array}{l}
\frac{d x_{i}}{d t}=\frac{D_{\mathrm{x}}}{L^{2} h^{2}}\left(x_{i-1}-2 x_{i}+x_{i+1}\right)+f\left(x_{i}, y_{i}\right), \\
\frac{d y_{i}}{d t}=\frac{D_{\mathrm{y}}}{L^{2} h^{2}}\left(y_{i-1}-2 y_{i}+y_{i+1}\right)+g\left(x_{i}, y_{i}\right)
\end{array}
\]

при $i=1,2, \ldots, n-1 ; h=1 / n, x_{0}=x_{n}=\bar{x}, y_{0}=y_{n}=\bar{y}$.
Для численного решения были выбраны следующие параметры задачи: $A=2, B=5,45, D_{x}=0,008, D_{y}=0,004$. В случае ГУ1 мы имеем тривиальное стационарное решение $x(t) \equiv$ $\equiv \bar{x}=A, y(t) \equiv \bar{y}=B / A$, устойчивое при $L \rightarrow 0$. С помощью методики, описанной в 6.3.1, находим критическую длину, отвечающую комплексной бифуркации $L_{1}^{+}=0,5130$ (вещественные бифуркации при изменении параметра $L$ не наблюдаются). Следовательно, тривиальное решение является устойчивым при $L \in\left(0, L_{1}^{+}\right)$и неустойчивым при $L>L_{1}^{+}$. При значении $L=L_{1}^{+}$ от ветви тривиальных решений отходит ветвь устойчивых периодических решений. Диаграмма периодических решений, построенная для вышеприведенной системы с помощью алгоритма продолжения DERPER, показана на рис. 6.15.

Добавим еще несколько замечаний по поводу этого рисунка. Как уже сказано, при значении параметра $L=L_{1}^{+}$от ветви стационарных решений отходит ветвь периодических решений (существующая при $L>L_{1}^{+}$). Этой ветви принадлежит, в частности, устойчивое периодическое решение при $L=0,75$, показанное на рис. 6.13. Периодические решения с указанной ветви, как это следует из рис. 6.13, являются пространственно симметричными: для всякого $t$ и $z \in[0,1]$ имеют место соотношения
1) «Смешанное» условие включает значения решения на обоих концах отрезка $0 \leqslant t \leqslant T$. – Прим. ред.

$x(z, t)=x(1-z, t) ; y(z, t)=y(1-z, t)$. В точке, обозначенной на рис. 6.15 буквами ТПС, периодические решения утрачивают устойчивость (происходит бифуркация типа «+1») и от исходной ветви симметричных решений отделяется ветвь устойчивых пространственно несимметричных решений. Заметим, что в данном случае имеются две такие ветви; они содержат периодические решения, которые (при заданном $L$ ) обладают одинаковым

Рис. 6.15. Диаграмма периодических решений для задачи 11 при ГУ1; зависимость периода $T$ ( $\left.A=2, B=5,45, D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004\right)$. Метод прямых (6.5.3), $n=20$. Сплошные линии – устойчивые периодические решения, штриховые – неустойчивые периодические решения, $T_{1}, T_{2}, T_{3}$ – точки бифуркации тила тора, ТПС – точка бифуркации с потерей симметрии.

периодом $T$. При этом для любого $t$ их пространственные профили являются взаимно симметричными, т. е. при произвольном $z \in[0,1] \quad$ выполняется $\quad x^{(1)}(z, t)=x^{(2)}(1-z, t), \quad y^{(1)}(z, t)=$ $=y^{(2)}(1-z, t)$, где верхний индекс определяет ветвь решений. В табл. 6.13 приведены 14 наибольших по абсолютной величине мультипликаторов, а также величины периода $T$, для неустойчивого периодического решения в случае $L=1,5$ (обозначенного на рис. 6.15 кружком) при двух различных значениях $n=20$ и $n=40$. Напомним, что при $n=20$ матрица монодромии имеет 38 собственных чисел, а при $n=40$ у нее 78 собственных чисел. Остальные собственные числа, которые не приведены в таблице, оказываются близкими к нулю. Из таблицы

Таблица 6.13. Собственные числа матрицы монодромии и период $T$ в зависимости от разбиения отрезка $0 \leqslant z \leqslant 1$. Периодическое решение при $L=1,5$ (на рис. 6.15 оно обозначено кружком).

видно, что главные мультипликаторы (по абсолютной величине большие $1 / 2$ ) для обоих разбиений оказываются практически одними и теми же. Поскольку $\left|\lambda_{2,3}\right|>1$, данное периодическое решение оказывается неустойчивым. Величина $\left|\lambda_{2,3}\right|$ невелика, и можно говорить о слабой неустойчивости: при динамическом моделировании с начальным условием вблизи этого решения соответствующая траектория будет очень медленно «разматываться», уходя от этого периодического решения.

На каждой из ветвей пространственно несимметричных решений была найдена точка бифуркации, при которой рождается инвариантный тор (обозначенная как $T_{1}$ на рис. 6.15 ; обе ветви на изображенной здесь диаграмме решений сливаются).
Таблица 6.14. Сравнение бифуркационных значений $L^{*}$ точки $T_{1}$ (рис. 6.15) при различных $n$.

В табл. 6.14 приведены бифуркационные значения $L^{*}$ для точки $T_{1}$, найденные при различных $n$. Значение $L^{*}$, подсчитанное с помощью экстраполяции по Ричардсону для $n=20$ и
Рис. 6.16а. Проекция орбиты отображения Пуанкаре при $L=1,4050$.
$n=40$, иллюстрирует точность определения этой величины в зависимости от $n$.
Рис. 6.16b. Проекция орбиты отображения Пуанкаре при $L=1,4070$.
На рис. 6.15 показаны, кроме того, еще две другие точки бифуркации рождения инвариантного тора $\left(T_{2}, T_{3}\right)$, лежащие на другой ветви периодических решений.

В окрестности точки бифуркации $T_{1}$ при $L>L^{*}$ мы с помощью динамического моделирования нашли устойчивый инвариантный тор. На рис. 6.16а изображена проекция орбиты отображения Пуанкаре для траектории на этом инвариантном торе. При возрастании величины $L$ указанный тор теряет устойчивость и рождается «удвоенный» устойчивый тор (рис. 6.16b). При этом отображение Пуанкаре в обоих случаях определяется уравнением $x(0,3 ; \tilde{t})=2$. При возрастании значений параметра $L$ можно наблюдать последовательность нескольких дальнейших «удвоений» тора. Учитывая слабую неустойчивость существующих одновременно периодических решений, можно было наблюдать медленное приближение траектории к устойчивому инвариантному тору. Если начальная точка выбиралась в окрестности периодического решения, то для того чтобы решение стабилизировалось на устойчивом инвариантном торе, нам приходилось интегрировать в течение времени порядка нескольких тысяч периодов.
6.5.1.2. Метод конечных разностей

Заменим дифференциальные уравнения (4.3.7) конечно-разностными по схеме Кранка-Николсона (6.4.3) при $w=0,5$ и $\tau=T / m$. Мы получим систему
\[
\begin{array}{c}
\frac{m}{T}\left(x_{i}^{j+1}-x_{i}^{j}\right)-\frac{D_{\mathrm{x}}}{L^{2} h^{2}}\left[x_{i-1}^{j+1 / 2}-2 x_{i}^{j+1 / 2}+x_{i+1}^{j+1 / 2}\right]- \\
-f\left(x_{i}^{i+1 / 2}, y_{i}^{i+1 / 2}\right)=0 \\
x_{i}^{j+1 / 2}=\frac{1}{2}\left(x_{i}^{j+1}+x_{i}^{j}\right) ; \quad y_{i}^{j+1 / 2}=\frac{1}{2}\left(y_{i}^{j+1}+y_{i}^{j}\right) \\
\frac{m}{T}\left(y_{i}^{j+1}-y_{i}^{j}\right)-\frac{D_{\mathrm{y}}}{L^{2} h^{2}}\left[y_{i-1}^{j+1 / 2}-2 y_{i}^{j+1 / 2}+y_{i+1}^{i+1 / 2}\right]- \\
-g\left(x_{i}^{j+1 / 2}, y_{i}^{j+1 / 2}\right)=0
\end{array}
\]

при $i=1,2, \ldots, n-1$ и $j=0,1, \ldots, m-1$. Граничные условия 1 -го рода, дают уравнения (6.4.4) для $j=-1,0,1, \ldots$ $\ldots, m-1$, а из условия периодичности (6.5.1) находим
\[
x_{i}^{0}-x_{i}^{m}=0, \quad y_{i}^{0}-y_{i}^{m}=0, \quad i=1,2, \ldots, n-1 .
\]

Тем самым в целом у нас имеется $N=2(n+1)(m+1)$ уравнений $(6.5 .4 \mathrm{a}),(6.5 .4 \mathrm{~b}),(6.5 .4 \mathrm{c}),(6.4 .4)$ относительно $N+1$ неизвестных $x_{i}^{i}, y_{i}^{j}, T$. Точно так же, как и в п. 5.8.2, зафиксируем какую-либо одну компоненту решения, например $x_{1}^{0}$, оста-

вив период $T$ в качестве неизвестной. Полученную систему нелинейных уравнений можно решать, например, с помощью метода Ньютона, причем системы линейных алгебраических уравнений, возникающие на каждой итерации, решаются с помощью подходящего алгоритма для почти ленточных матриц. Получающиеся при этом системы имеют большие размеры, и поэтому этот метод удобно реализовывать лишь на достаточно мощных ЭВМ и при не слишком больших размерах $m$ и $n$. Периодическое решение задачи 11, найденное таким методом для относительно редкой сетки узловых точек, представлено в табл. 6.15.

Таблица 6.15. Периодическое решение задачи 11 с граничными условиями типа ГУ 1 , найденное из конечно-разностных уравнений $(6.5 .4) ; A=2$, $B=5,45, D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004, L=0,75, m=19, n=5$. Окончательное значение $T=3,3073$; в качестве фиксированного значения одной из неизвестных выбиралось $y_{4}^{0}=3,4$. Найденное решение является симметричным по переменной $z$ относительно точки $z=0,5$; поэтому $x_{3}^{j}=x_{2}^{j}, x_{4}^{j}=x_{1}^{j}, y_{3}^{j}=y_{2}^{j}$ $y_{4}^{j}=y^{j}$.

На рис. 6.17 показан еще один способ изображения периодического решения в распределенных системах (ср. с рис. 6.13). Этот способ удобен в тех случаях, когда профили решения качественно сохраняют свою форму и изменяются не слишком сильно (так, чтобы на рисунке не перекрывать друг друга).

Для продолжения периодического решения при изменении параметра, например параметра $L$, можно разработать алгоритм, основанный на конечно-разностной схеме (6.5.4). Вблизи точек поворота величину $L$ необходимо рассматривать как неизвестную и последовательно изменять значения какой-либо другой переменной, например $T$.

6.5.2. Решения волнового характера

Рассмотрим сначала уравнения с частными производными (4.3.7) без граничных условий на бесконечном интервале изменения независимой переменной $z$. Подстановка
\[
x(z, t)=\varphi(z-c t), \quad y(z, t)=\psi(z-c t)
\]

в уравнения (4.3.7) приводит к системе двух обыкновенных диф-

Рис. 6.17. Периодическое решение задачи 11 при ГУ1, найденное с помощью разностных уравнений (6.5.4): $m=18, n=5, A=2, B=5,45, D_{\mathrm{x}}=0,008$, $D_{y}=0,004, T=3,3073$.

ференциальных уравнений 2-го порядка для функций $\varphi(\xi)$ и $\psi(\xi):$
\[
\begin{array}{l}
-c \varphi^{\prime}=\frac{D_{\mathrm{x}}}{L^{2}} \varphi^{\prime \prime}+f(\varphi, \psi),\left(^{\prime}=\frac{d}{d \xi}\right), \\
-c \psi^{\prime}=\frac{D_{\mathrm{y}}}{L^{2}} \psi^{\prime \prime}+g(\varphi, \psi) .
\end{array}
\]

Периодические решения, а также гомоклинические и гетероклинические траектории системы (6.5.6) соответствуют решениям волнового типа исходной системы.

Обратимся теперь к случаю конечного промежутка изменения переменной $z$ и рассмотрим вновь систему (4.3.7) с соответствующими граничными условиями.

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

В первой части этого пункта мы рассмотрим решения типа импульса, а во второй части исследуем волну типа фронта (см. $\S 2.6$ ).

6.5.2.1. Волна типа импульса
Волновые решения уравнений (4.3.7) зависят от свойств системы, получающейся при $D_{\mathrm{x}} \rightarrow 0, \quad D_{\mathrm{y}} \rightarrow 0$. Рассмотрим задачу 12, т. е. зададим $f$ и $g$ формулами (Р12-1) и (Р12-2). Выберем следующие значения параметров: $\gamma=3, v_{0}=0,01, \beta=$

Рис. 6.18. Фазовый портрет системы (P4-1), (P4-2): $\gamma=3, v_{0}=0,01, \beta=$ $=1,5, \delta=1, \alpha=12 ; \mathrm{S}_{1}, \mathrm{~S}_{2}, \mathrm{~S}_{3}$ – стационарные решения.
$=1,5, \quad \delta=1, \alpha=12$. Соответствующая система при $D_{\mathrm{x}}=$ $=D_{\mathrm{y}}=0$ (см. (P4-1), (Р4-2)) будет иметь три стационарных решения:
\[
\begin{array}{l}
\mathbf{S}_{1}: x_{1}=0,11566, y_{1}=0,19618 \text { (устойчивый фокус); } \\
\mathbf{S}_{2}: x_{2}=0,37278, y_{2}=0,89151 \text { (седло); } \\
\mathbf{S}_{3}: x_{3}=0,70063, y_{3}=3,51059 \text { (неустойчивый узел). }
\end{array}
\]

На рис. 6.18 кроме стационарных состояний показаны нулевые изоклины функций $f$ и $g$. Пунктиром изображена одна из траекторий. Из рисунка видно, что относительно малое возмущение стационар ного состояния $\mathbf{S}_{1}$ может вызвать изменение переменных состояния по «циклу возбуждения»; система уходит при этом далеко от устойчивого стационарного состояния.

(Если выбрать начальное условие справа от $S_{1}$ за нулевой изоклиной функции $f$, то фазовая траектория закручивается вокруг состояния $\mathbf{S}_{3}$, и лишь потом устремляется $\mathrm{K}_{1}$.)

Выберем теперь для уравнений (4.3.7) в задаче 12 следующие граничные условия:
\[
\begin{array}{c}
z=0: x(0, t)=x_{\mathrm{p}}, \quad y(0, t)=y_{1}, \\
z=1: \frac{\partial x(1, t)}{\partial z}=0, \quad \frac{\partial y(1, t)}{\partial z}=0 .
\end{array}
\]

Если $x_{\mathrm{p}}=x_{1}$ ( $x_{1}$ и $y_{1}$ представляют собой значения $x$ и $y$ в стационарном состоянии $\mathbf{S}_{1}$, см. выше), то задача имеет тривиальное решение $x(z, t) \equiv x_{1}, y(z, t) \equiv y_{1}$, которое является устойчивым. Таким образом, величина $x_{\mathrm{p}}$ характеризует собой

Рис. 6.19. Бегущая волна в задаче 12 при ГУ2: $\gamma=3, v_{0}=0,01, \beta=1,5$, $\delta=1, \alpha=12, D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004, L=2,5, x_{\mathrm{p}}=2 ; v=0,108, T=8,0$.

возмущение переменной $x$ на левом конце промежутка. На правом конце заданы граничные условия второго рода, которые описывают непроницаемую для исследуемых веществ стенку. Уравнения (4.3.7) решались посредством конечно-разностной схемы типа Кранка – Николсона с заменой нелинейности с помощью отрезка ряда Тейлора (см. § 6.4) при $n=160, \tau=0,1$. Если выбирать $x_{\mathrm{p}}$ не слишком близко к $x_{1}$, то в системе начинают возникать волны концентрации, бегущие от $z=0$ до $z=1$. Одна из таких волн изображена на рис. 6.19 (для некоторого момента времени $t$ ). Эта волна возникла вблизи левого конца промежутка, в точке $z=0,2$ ее форма оказалась уже достаточно развитой, далее она сохраняет свой вид до значения ${ }^{1)} z=0,8$, после чего на правом конце промежутка она исчезает в силу действия граничных условий (6.5.7b). Определим скорость движения волны как $v=d z_{v} / d t$, где $z_{v}$ есть, например, координата вершины волны. Далее, определим период

серии импульсов как время между двумя последовательными прохождениями импульса через заданную точку $z$, например $z=0,5$. Значения $v$ и $T$ для волны, изображенной на рис. 6.19, указаны в подписи к рисунку (они были найдены вычислениями из результатов динамического моделирования). При более высоких значениях $x_{\text {p }}$ волны уже не возникают и система стабилизируется на устойчивом неоднородном стационарном решении. Более подробные результаты можно найти, например, в работе [6.31].

6.5.2.2. Волна типа фронта
Фронт представляет собой переход решения системы типа «реакция – диффузия» от одного устойчивого стационарного решения к другому, причем в переходном режиме одна часть системы характеризуется значениями одного решения, а другая – значениями другого решения. Между этими двумя частями располагается узкая область (фронт), где наблюдается резкая зависимость физических переменных задачи от пространственной координаты. Указанный фронт перемещается, в результате чего размеры обеих областей изменяются до тех пор, пока система не достигнет конечного стационарного состояния, и фронт исчезнет.

Волну типа фронта мы снова продемонстрируем на примере задачи 12 при значениях параметров $\gamma_{3}=3, v_{0}=0,01, \beta=1,5$, $\delta=1,7, \alpha=12$. Стационарные решения системы (P4-1), (P4-2) в данном случае имеют вид:
\[
\begin{array}{l}
\mathbf{S}_{1}: x_{1}=0,12985, y_{1}=0,12404 \text { (устойчивый узел); } \\
\mathbf{S}_{2}: x_{2}=0,26666, y_{2}=0,27906 \text { (седло); } \\
\mathbf{S}_{3}: x_{3}=1,32824, y_{3}=5,33927 \text { (устойчивый фокус). }
\end{array}
\]

Если для распределенной системы (4.3.7) выбрать граничные условия второго рода, т. е. условия (4.3.12), то тривиальные стационарные решения
\[
\begin{array}{l}
x(z) \equiv x_{1}, \quad y(z) \equiv y_{1}, \\
x(z) \equiv x_{3}, \quad y(z) \equiv y_{3}
\end{array}
\]

оказываются устойчивыми при $D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004, L=2,5$, причем от начальных условий зависит, на каком из этих решений система стабилизируется.

Переход системы из состояния (6.5.8) в состояние (6.5.9) моделировался путем изменения начального условия для переменных $x$ и $y$ в точке $z=0$ в течение промежутка времени $t_{1}$ :
\[
t \in\left[0, t_{1}\right]: x(0, t)=x_{3}, \quad y(0, t)=y_{3} .
\]

Для моментов времени $t>t_{1}$ вновь действуют граничные условия второго рода. Состояние системы (в некоторый момент времени) для $t_{1}=0,1$ показано на рис. 6.20. Скорость движения фронта при этом составляла $v \sim 0,115$; весь переход от (6.5.8) к (6.5.9) продолжался приблизительно 22 временны́х единицы. При очень малых значениях $t_{1}$ система возвращается к стационарному решению (6.5.8). Более подробно эти результаты обсуждаются в работе [6.31]. Отметим также, что
Рис. 6.20. Волна типа фронта для задачи 12 при ГУ2: $\gamma=3, v_{0}=0,01, \beta=$ $=1,5, \delta=1,7, \alpha=12, D_{\mathrm{x}}=0,008, D_{\mathrm{y}}=0,004, L=2,5$; штриховая линия – решение (6.5.8), штрихпунктирная – решение (6.5.9).
переходы из одного стационарного состояния в другое могут происходить либо легко (под влиянием малого возмущения, действующего в течение короткого времени), либо с трудом, если на границе системы вообще не реализуемы физически допустимые возмущения.