Так же как и для «сосредоточенных» систем (описываемых обыкновенными дифференциальными уравнениями, см. § 5.2),
Рис. 6.2. Диаграмма стационарных решений задачи 17, $\mathrm{Re}=625$.

существует целый ряд методов и подходов, позволяющих находить зависимость стационарных решений от параметров в случае «распределенных» систем, описываемых уравнениями

с частными производными. Набор методов оказывается здесь даже более широким, поскольку здесь есть еще и проблема аппроксимации задачи в бесконечномерном пространстве посредством задачи в пространстве конечной размерности.

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

Результаты последовательного использования метода Ньютона для системы разностных уравнений (6.1.23), которая возникает в задаче 17, представлены на диаграмме стационарных решений, показанной на рис. 6.2. В данном случае уже само решение системы (6.1.23) для больших $n$ являлось достаточно серьезной проблемой.

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

Отметим, что методы, изложенные в п. 6.2 .1 и 6.2 .2 , могут использоваться в качестве предиктора в схеме типа «предиктор-корректор», описанной в п. 6.2.3, аналогично тому, как это делалось в § 5.2.

6.2.1. Метод дифференцирования по параметру

Рассматриваемая здесь методика аналогична описанной в п. 5.2.2 для конечномерных задач. Продемонстрируем здесь одну из возможных численных реализаций этого метода на примере продолжения по параметру $\alpha$ решения краевой задачи для (одного) дифференциального уравнения второго порядка $(‘=$ $=d / d z)$
\[
y^{\prime \prime}-f\left(z, y, y^{\prime}, \alpha\right)=0
\]

с граничными условиями
\[
\begin{array}{l}
a_{0} y(0)+b_{0} y^{\prime}(0)=c_{0}, \\
a_{1} y(1)+b_{1} y^{\prime}(1)=c_{1} .
\end{array}
\]

Подставив в уравнение (6.2.1) $y=y(z, \alpha)$ и дифференцируя по $\alpha$, мы получаем уравнение
\[
\frac{\partial^{3} y}{\partial z^{2} \partial \alpha}-\frac{\partial f}{\partial y} \frac{\partial y}{\partial \alpha}-\frac{\partial f}{\partial y^{\prime}} \frac{\partial^{2} y}{\partial z \partial \alpha}-\frac{\partial f}{\partial \alpha}=0
\]

с граничными условиями вида
\[
\begin{array}{l}
a_{0} y(0, \alpha)+b_{0} \frac{\partial y(0, \alpha)}{\partial z}=c_{0}, \\
a_{1} y(1, \alpha)+b_{1} \frac{\partial y(1, \alpha)}{\partial z}=c_{1} .
\end{array}
\]

Предположим, что известно решение краевой задачи (6.2.1)(6.2.2) при $\alpha=\alpha_{0}$
\[
y\left(z, \alpha_{0}\right)=\varphi(z) .
\]

Дифференциальное уравнение в частных производных третьего порядка (6.2.3) можно решать разностными методами, аналогично простейшим уравнениям параболического типа (см. $\S 6.4)$. Выберем прямоугольную сетку узловых точек $\left(z_{i}, \alpha_{j}\right)$, $. z_{i}=i h, i=0,1, \ldots, n ; \alpha_{j}=\alpha_{0}+j k, j=0,1, \ldots$, и обозначим $y_{i}^{j} \sim y\left(z_{i}, \alpha_{j}\right)$. Простейшая разностная аппроксимация уравнения (6.2.3) с использованием шеститочечной схемы имеет вид
\[
\begin{array}{c}
\frac{1}{k h^{2}}\left(y_{i-1}^{i+1}-2 y_{i}^{j+1}+y_{i+1}^{i+1}-y_{i-1}^{i}+2 y_{i}^{i}-y_{i+1}\right)-\bar{f}_{y} \frac{1}{k}\left(y_{i}^{i+1}-y_{i}^{i}\right)- \\
-\vec{f}_{y^{\prime}} \frac{1}{2 k h}\left(y_{i+1}^{i+1}-y_{i-1}^{i+1}-y_{i+1}^{i}+y_{i-1}^{i}\right)-\vec{f}_{\alpha}=0, \quad(6.2 .6) \\
i=1,2, \ldots, n-1 .
\end{array}
\]

Здесь $\bar{f}_{y}, \bar{f}_{y^{\prime}}, \bar{f}_{\alpha}$ – соответствующие частные производные функции $f$, вычисленные, например, для значений переменных $y_{i}$ на $j$-м (уже найденном) «слое»: так,
\[
\bar{f}_{y^{\prime}}=\frac{\partial f\left(z_{i}, y_{i}^{j},\left(y_{i+1}^{i}-y_{i-1}^{j}\right) / 2 h, \alpha_{i}\right)}{\partial y^{\prime}} .
\]

Используя подходящие разностные замены в граничных условиях (6.2.4) (см. п. 6.1.1 и § 6.4), мы получаем, вместе с соотношением (6.2.6), систему линейных алгебраических уравнений

относительно неизвестных $y_{0}^{i+1}, y_{1}^{j+1}, \ldots, y_{n}^{i+1} \quad$ (переменные с верхним индексом $j$ нам известны). Матрица системы является трехдиагональной, и поэтому мы можем воспользоваться: алгоритмом исключения Гаусса, работающим только с ненулевыми элементами.

Поскольку разностная замена (6.2.6) центрирована относительно точки ( $i, j+1 / 2)$, то уменьшения погрешностей аппроксимации можно достичь, заменяя коэффициенты $\bar{f}_{y}, \bar{f}_{y^{\prime}}, \bar{f}_{\alpha}$ в этой. точке, т. е. вводя в соотношения (6.2.7) аргументы
\[
\left(z_{i}, \frac{y_{i}^{i}+y_{i}^{j+1}}{2}, \frac{y_{i+1}^{i}+y_{i+1}^{j+1}-y_{i-1}^{i}-y_{i-1}^{j+1}}{4 h}, \frac{\alpha_{j}+\alpha_{i+1}}{2}\right) .
\]

В результате мы получаем систему нелинейных уравнений относительно неизвестных с верхним индексом $j+1$, для решения которой можно использовать какой-либо итерационный метод, например, метод Ньютона.

Уравнение (6.2.3) можно также аппроксимировать методом прямых (см. §6.4). Более подробный анализ метода дифференцирования по параметру можно найти в книге [6.8].

Выше мы сначала проводили дифференцирование по параметру $\alpha$, а затем решали задачу численно (дискретизировали ее). Обратный порядок действий также возможен. Так, сначала мы можем преобразовать задачу (6.2.1), (6.2.2) с помощью соответствующих разностных формул, например, приводя еек виду
\[
\begin{array}{c}
\frac{y_{i-1}-2 y_{i}+y_{i+1}}{h^{2}}-f\left[z_{i}, y_{i}, \frac{y_{i+1}-y_{i-1}}{2 h}, \alpha\right]=0, \\
i=0,1, \ldots, n, \\
a_{0} y_{0}+b_{0}\left(y_{1}-y_{-1}\right) / 2 h=c_{0}, \\
a_{1} y_{n}+b_{1}\left(y_{n+1}-y_{n-1}\right) / 2 h=c_{1} .
\end{array}
\]

Переменные $y_{-1}$ и $y_{n+1}$ можно найти из условий (6.2.10), а затем подставить их в уравнение (6.2.9). К полученной таким образом системе $n+1$ нелинейных уравнений $f_{i}\left(y_{0}, y_{1}, \ldots, y_{n}\right)=0$, $i=0,1, \ldots, n$ можно затем применить метод п. 5.2.2. Интегрируя дифференциальные уравнения (5.2.8) методом Эйлера, мы получаем алгоритм, сходный с тем, который описывался выше.

Точно так же, как и в случае задач конечной размерности, здесь метод дифференцирования по параметру $\alpha$ не позволяет преодолевать точку поворота на диаграмме решений. Метод. дифференцирования по длине дуги в том виде, как он был описан в п. 5.2.2, позволяет спокойно проходить точку поворота,

хотя для задач типа (6.2.9) при большом $n$ без преобразований, использующих специальные структуры матрицы Якоби, этот метод применять невозможно. Использование указанного метода будет обсуждаться ниже, в п. 6.2.3, а также в связи с методом дифференцирования по граничному условию (п. 6.2.2).

Продемонстрируем теперь применение метода дифференцирования по параметру на примере задачи 16. В предыдущем пункте мы вывели уравнение (6.1.21), которое можно представить в виде (6.2.1), полагая $z=r$ и $\alpha=\Phi$. Если теперь обо:значить
\[
E(y)=\exp \frac{\gamma \beta(1-y)}{1+\beta(1-y)},
\]
‘то уравнение (6.2.3) для (6.1.21) принимает вид
\[
\begin{array}{c}
\frac{\partial^{3} y}{\partial r^{2} \partial \Phi}+\frac{a}{r} \frac{\partial^{2} y}{\partial r \partial \Phi}-\Phi^{2}\left[1-\frac{\gamma \beta y}{(1+\beta(1-y))^{2}}\right] E(y) \frac{\partial y}{\partial \Phi}- \\
-2 \Phi y \cdot E(y)=0 .
\end{array}
\]

Граничные условия при $\mathrm{Nu} \rightarrow \infty, \mathrm{Sh} \rightarrow \infty$ записываются как
\[
\frac{\partial y(0, \Phi)}{\partial r}=0, \quad y(1, \Phi)=1 .
\]

Начальное условие (6.2.5) легко получается аналитически в виде
\[
\Phi=0: y(r, 0) \equiv 1 .
\]

При этом разностная аппроксимация уравнения (6.2.11) проводилась в соответствии с формулами (6.2.6) и (6.2.7), а замена граничного условия в точке $r=0$ осуществлялась точно так же, как в уравнении (6.1.22b). Влияние выбора шагов сетки $h$ и $k$ иллюстрируется табл. 6.6. Заметим, что для получения значения $y(0,1)$ при $h=0,025$ и $k=0,002$ нам пришлось 500 раз:

Таблица 6.6. Метод дифференцирования по параметру Ф для задачи 16, $a=0, \gamma=20, \beta=0,05$. Приведены значения $y(0,1)$, т. е. значения концентрации в центре частицы при $\Phi=1$ (точное значение равно 0,5521 ). Видно влияние шагов $h$ и $k$.

решить систему 40 линейных алгебраических уравнений с трехдиагональной матрицей. Тем самым было найдено 500 точек на диаграмме, описывающей зависимость стационарного решения от параметра.

Применение метода дифференцирования по параметру мы продемонстрировали на примере одного уравнения второго порядка. Переход к системе таких уравнений также не представляет особых трудностей – при этом лишь увеличится размерность решаемых задач и вместо трехдиагональной мы будем иметь дело с ленточной матрицей. Читатель может легко вывести соответствующие уравнения в частных производных для задач 11, 12, 13, 14. Задача 17 была решена этим методом в работе [6.12]. В следующем пункте описывается модификация данного метода, с помощью которой мы можем проходить точки поворота на диаграмме стационарных решений.

6.2.2. Метод дифференцирования по граничному условию

Обратимся вновь к задаче (6.2.1), (6.2.2), предполагая, что $b_{0}
eq 0$. В качестве нового «параметра» задачи $\xi$ введем значение решения в точке $z=0$ :
\[
y(0)=\xi \text {. }
\]

Решение $y$ задачи (6.2.1), (6.2.2) будет зависеть от пространственной переменной $z$ и параметра $\xi$, т. е. $y=y(z, \xi)$. Дифференцируя уравнение (6.2.1) по $\xi$, находим
\[
\frac{\partial^{3} y}{\partial z^{2} \partial \xi}-\frac{\partial f}{\partial y} \frac{\partial y}{\partial \xi}-\frac{\partial f}{\partial y^{\prime}} \frac{\partial^{2} y}{\partial z \partial \xi}-\frac{\partial f}{\partial \alpha} \frac{d \alpha}{d \xi}=0 .
\]

Граничные условия (6.2.2) принимают вид
\[
a_{0} y(0, \xi)+b_{0} \frac{\partial y(0, \xi)}{\partial z}=c_{0}, \quad a_{1} y(1, \xi)+b_{1} \frac{\partial y \cdot(1, \xi)}{\partial z}=c_{1},
\]

причем из формулы (6.2.14) следует, что
\[
y(0, \xi)=\xi \text {. }
\]

Теперь $\alpha$ есть функция от $\xi$, и в уравнение (6.2.15) входит неизвестная величина
\[
C=d \alpha / d \xi .
\]

Поэтому три граничных условия (6.2.16) не делают задачу переопределенной. Начальное условие для продолжения по параметру мы задаем в форме
\[
\xi=\xi_{0}: \alpha=\alpha_{0}, \quad y\left(z, \xi_{0}\right)=\varphi(z),
\]

тде функция $\varphi(z)$ есть решение исходной краевой задачи $(6.2 .1),(6.2 .2)$ при $\alpha=\alpha_{0}$ и $\varphi(0)=\xi_{0}$. В этом случае функцию $\varphi(z)$ часто можно найти в аналитическом виде для некоторого промежуточного значения параметра $\alpha$.

Дифференциальное уравнение в частных производных (6.2.15), точно так же, как и в предыдущем пункте, нетрудно аппроксимировать соответствующими разностными уравнениями. Обозначим $z_{i}=i h, i=0,1, \ldots, n, \xi_{j}=\xi_{0}+j k, j=1,2, \ldots$, $y\left(z_{i}, \xi_{j}\right) \sim y_{i}^{i}, \quad \alpha\left(\xi_{j}\right) \sim \alpha_{j}, \quad h=1 / n$. Тогда разностная аппроксимация уравнений (6.2.15) и (6.2.16) приводит к соотношениям
\[
\begin{array}{c}
s_{i-1}^{i} y_{i-1}^{i+1}+s_{i}^{i} y_{i}^{j+1}+s_{i+1}^{j} y_{i+1}^{i+1}+r_{i}^{i} C=q_{i}^{i}, \quad i=0,1, \ldots, n,(6.2 .19 \mathrm{a}) \\
a_{0} y_{0}^{j+1}+b_{0}\left(y_{1}^{j+1}-y_{-1}^{j+1}\right) / 2 h=c_{0}, \\
a_{1} y_{n}^{j+1}+b_{1}\left(y_{n+1}^{i+1}-y_{n-1}^{j+1}\right) / 2 h=c_{1}, \\
y_{0}^{j+1}=\xi_{j}+k .
\end{array}
\]

Коэффициенты $s_{i-1}^{i}, s_{i}^{i}, s_{i+1}^{i}, r_{i}^{i}, q_{i}^{i}$ вычисляются с помощью известных значений решения па предыдущем «слое» при $\xi=\xi_{j}$ (т. е. при $\alpha=\alpha_{j}$ ). Соотношения (6.2.19) представляют собой систему $n+4$ линейных алгебраических уравнений относительно $n+4$ неизвестных $y_{-1}^{i+1}, \ldots, y_{n+1}^{i+1}, C$ с почти ленточной матрицей. Решив эту систему, мы получим значения решения $y^{j+1}$ для $\xi=\xi_{j+1}=\xi_{j}+k$. Для нахождения $\alpha$ воспользуемся методом Эйлера, т. е. формулой
\[
\alpha_{i+1}=\alpha_{j}+k C,
\]

которая следует из уравнения (6.2.17).
Продемонстрируем применение метода дифференцирования по граничному условию на примере задачи 16. По аналогии с уравнением (6.2.11) найдем производную уравнения (6.1.21) по $\xi=y(0)$, рассматривая его как уравнение относительно $y(r, \xi), \quad \delta(\xi) \quad\left(\delta=\Phi^{2}\right.$, при сравнении с общими формулами (6.2.1) и (6.2.15) $z=r$ и $\alpha=\delta$ ):
\[
\begin{array}{c}
\frac{\partial^{3} y}{\partial r^{2} \partial \xi}-\frac{a}{r} \frac{\partial^{2} y}{\partial r \partial \xi}-\delta\left[1-\frac{\gamma \beta y}{(1+\beta(1-y))^{2}}\right] E(y) \frac{\partial y}{\partial \xi}- \\
-y E(y) \frac{d \delta}{d \xi}=0
\end{array}
\]

Граничные условия записываются в виде
\[
\frac{\partial y(0, \xi)}{\partial r}=0, \quad y(1, \xi)=1, \quad y(0, \xi)=\xi .
\]

При $\delta_{0}=0$ решение задачи (6.2.1) и (6.2.2) (здесь-уравнения (6.1.21) при условиях $y^{\prime}(0)=0, \quad y(1)=1$ ) имеет вид. $y(r) \equiv 1$ и, следовательно, можно положить
\[
\xi_{0}=1: y(r, 1) \equiv 1, \quad \delta(1)=0 .
\]

В табл. 6.7 показано влияние шагов $h$ и $k$ на точность получаемых результатов. Ясно, что найденные значения решения даже при сравнительно большом числе шагов по переменной $\xi$ не: слишком точны. Это лишь подтверждает наши выводы, сделанные в п. 5.2.2, о том, что сам по себе метод дифференцирования по параметру (или по граничному условию) не доставляет эффективного алгоритма для продолжения решения по параметру. Эффективный алгоритм, в котором вышеприведенные методы. используются в качестве предиктора, а метод Ньютона – в качестве корректора, будет рассмотрен нами в следующем пункте.
Таблица 6.7. Метод дифференцирования по граничному условию для задачи $16(a=0, \gamma=20, \beta=0,05)$. Приведены значения параметра $\delta=\Phi^{2}$ при $\xi=y(0)=0,5$ (точное значение равно $\delta=1,1546$ ). Видно влияние шагов $h$ и $k$.

В заключение этого пункта отметим, что метод дифференцирования по граничному условию позволяет проходить точки поворота на диаграмме решений (построенной по параметру $\alpha$ ), но не способен, однако, преодолевать точки, где $d \xi / d \alpha=0$ (т. е. $C=\infty$ ). Более подробный анализ метода дифференцирования по граничному условию читатель может найти в работе [6.8].

6.2.3. Алгоритм продолжения типа предиктор – корректор

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

щим образом сформулированной задачи, а функцию предиктора возьмут на себя описанные выше методы. Это и будет первый вариант алгоритма продолжения типа «предиктор-корректор», основанный на соответствующих разностных заменах исходной задачи (п. п. 6.2.3.1). Другой вариант этого алгоритма (п. п. 6.2.3.2) основывается на использовании метода стрельбы.
6.2.3.1. Алгоритм, основанный на разностных методах

В качестве примера рассмотрим вновь задачу (6.2.1) и (6.2.2). Точно так же как и в п. 5.2.2, введем новый параметр $u$, полагая при этом
\[
\Omega(z, u)=\frac{\partial y}{\partial u}, \quad C(u)=\frac{d \alpha}{d u} .
\]

Дифференцируя уравнения (6.2.1) и (6.2.2) по $u$, получаем
\[
\begin{array}{c}
\Omega^{\prime \prime}-\frac{\partial f}{\partial y} \Omega-\frac{\partial f}{\partial y^{\prime}} \Omega^{\prime}-\frac{\partial f}{\partial \alpha} C=0, \\
a_{0} \Omega(0)+b_{0} \Omega^{\prime}(0)=0, \quad a_{1} \Omega(1)+b_{1} \Omega^{\prime}(1)=0 .
\end{array}
\]

Заметим, что $u$ играет здесь роль переменной $z$ из уравнений (5.2.12). Для того чтобы $и$ было длиной дуги в пространстве $\mathrm{B} \times \mathrm{R}^{1}$, где $y \in \mathrm{B}$ и $\alpha \in \mathrm{R}^{1}$, должно выполняться условие, аналогичное соотношению (5.2.13), а именно ${ }^{1)}$
\[
\|\Omega\|^{2}+C^{2}=1 \text {. }
\]

Численная реализация этого равенства при разностной аппроксимации уравнений (6.2.25), (6.2.26) приводит к алгоритму, почти идентичному алгоритму DERPAR, описанному в п. 5.2.3. Размерность решаемых задач, однако, может оказаться значительной (при большом числе узловых точек или в случае, когда вместо одного уравнения (6.2.1) у нас имеется система дифференциальных уравнений). Кроме того, здесь всегда приходится использовать специальные методы решения возникающих разреженных систем линейных алгебраических уравнений со специальной (почти ленточной) структурой.

Ниже будет рассмотрена упрощенная схема данного алгоритма. Используя обозначения (6.2.24), выберем предиктор в форме метода Эйлера
\[
y(z, u+\Delta u)=y(z, u)+\Delta u \Omega(z, u), \quad \alpha(u+\Delta u)=\alpha(u)+\Delta u C(u),
\]

где $\Omega$ и $C$ вычисляются из задачи (6.2.25), (6.2.26) при $y=$ $=y(z, u)$ и $\alpha=\alpha(u)$. При
\[
C(u)= \pm 1
\]

мы получаем метод, формально аналогичный методу дифференцирования по параметру (п. 6.2.1), причем $u$ играет здесь. роль параметра $\alpha$ (на том участке, где $C=-1$, роль параметра $\alpha$ выполняет величина – $u$ ). Если теперь положить
\[
\Omega(0, u)= \pm 1,
\]

то мы получаем вариант метода дифференцирования по граничному условию, описанного в п. 6.2.2. Параметр $u$ играет здесь роль задаваемого граничного значения $\xi=y(0)$, или, соответственно, 一 . Шаг $\Delta u$ в (6.2.28) можно выбирать в зависимости от того, используем ли мы соотношения (6.2.29) или (6.2.30).

Очевидно, что выбор типа (6.2.29) будет невыгоден в окрестности точки поворота и невозможен в самой этой точке ${ }^{1)}$. И, наоборот, выбор типа (6.2.30) будет непригодным в окрестности экстремума зависимости $y(0)$ от параметра $\alpha$ на диаграмме решений. Из этих соображений вытекает алгоритм выбора между двумя указанными возможностями: если последнее вычисленное значение удовлетворяет условию
\[
\Omega(0, u)<M C,
\]

то для следующего шага принимается условие (6.2.29). В противоположном случае выбирается условие (6.2.30). Конечно, в обоих этих случаях мы должны соблюдать ориентацию вдоль кривой ${ }^{2)}$, т. е. знаки «十» или «-》 определяются с учетом предыдущих значений $C$ или $\Omega(0)$. Постоянную $M$ в условии (6.2.31) можно полагать равной 1 , если $y(0)$ и $\alpha$ соизмеримы, в противном же случае она выбирается с учетом того, какие относительные изменения $y(0)$ или $\alpha$ ожидаются вдоль ветви соответствующей диаграммы решений.

Таким образом, мы разрешили проблему предиктора. В качестве корректора, как и в § 5.2, можно использовать метод Ньютона, применяя его к системе разностных уравнений (6.2.9), (6.2.10). При этом вновь необходимо различать два случая. В первом из них выполнялось условие (6.2.31), и тогда предиктор выбирался в соответствии с формулой (6.2.29), во

втором – имела место противоположная ситуация. В первом случае мы решаем уравнения (6.2.9), (6.2.10) при фиксированном (предсказанном) значении параметра $\alpha$; во втором случае – при фиксированном (предсказанном) значении $y_{0}=y(0)$, а параметр $\alpha$ считается неизвестным. Блок-схему предложенного алгоритма продолжения типа «предиктор-корректор» читатель легко может построить самостоятельно.

6.2.3.2. Алгоритм, основанный на методе стрельбы

Идею этого метода мы теперь продемонстрируем на примере задачи (6.1.1), (6.1.2) с граничными условиями (6.1.4). Зададим снова два недостающих начальных условия в точке $z=0$ (см. (6.1.24)): $x(0)=\eta_{1}, y(0)=\eta_{2}$. Интегрируя уравнения (6.1.1), (6.1.2) с начальными условиями (6.1.4a), (6.1.24) при некотором значении параметра $L$, мы должны получить в точке $z=1$
\[
F_{1}\left(\eta_{1}, \eta_{2}, L\right)=x^{\prime}(1, \eta, L)=0, \quad F_{2}\left(\eta_{1}, \eta_{2}, L\right)=y^{\prime}(1, \eta, L)=0 .
\]

Таким образом, у нас имеются два уравнения относительно трех неизвестных $\eta_{1}, \eta_{2}, L$. Для нахождения соответствующей кривой в пространстве $\left(\eta_{1}, \eta_{2}, L\right)$ мы можем использовать алгоритм DERPAR, описанный в § 5.2. Для вычисления матрицы Якоби
\[
\left[\begin{array}{l}
\frac{\partial F_{1}}{\partial \eta_{1}}, \frac{\partial F_{1}}{\partial \eta_{2}}, \frac{\partial F_{1}}{\partial L} \\
\frac{\partial F_{2}}{\partial \eta_{1}}, \frac{\partial F_{2}}{\partial \eta_{2}}, \frac{\partial F_{2}}{\partial L}
\end{array}\right]
\]

используем вариационные переменные $p_{x i}, p_{y i}$, определяемые соотношениями (6.1.26) – (6.1.28), а также переменные
\[
p_{x L}=\partial x / \partial L, \quad p_{y L}=\partial y / \partial L,
\]

для которых нетрудно получить следующие уравнения в вариациях:
\[
\begin{array}{l}
p_{x L}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{x}}}\left(\frac{\partial f}{\partial x} p_{x L}+\frac{\partial f}{\partial y} p_{y L}\right)+\frac{2 L}{D_{\mathrm{x}}} f(x, y)=0, \\
p_{y L}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{y}}}\left(\frac{\partial g}{\partial x} p_{x L}+\frac{\partial g}{\partial y} p_{y L}\right)+\frac{2 L}{D_{\mathrm{y}}} g(x, y)=0
\end{array}
\]

с начальными условиями
\[
p_{x L}(0)=p_{x L}^{\prime}(0)=p_{y L}(0)=p_{y L}^{\prime}(0)=0 .
\]
Частные производные в уравнении (6.2.33) задаются формулами (6.1.29), а также соотношениями
\[
\frac{\partial F_{1}}{\partial L}=p_{x L}^{\prime}(1), \quad \frac{\partial F_{2}}{\partial L}=p_{y L}^{\prime}(1) .
\]

Таким образом, мы имеем все необходимое для использования алгоритма продолжения DERPAR (см. п. 5.2.3): можем вычислить функции $F_{1}$ и $F_{2}$ по формулам (6.2.32), а также матрицу Якоби (6.2.33). Для одного такого вычисления нам необходимо решить задачу Коши для восьми дифференциальных уравнений второго порядка $(6.1 .1),(6.1 .2),(6.1 .27)$ и (6.2.35), т. е. для системы 16 -го порядка.

В случае ГУ1 алгоритм изменяется в соответствии с формулами (6.1.30), (6.1.31), (6.1.32); при этом частные производные, входящие в матрицу Якоби (6.2.33), имеют вид
\[
\begin{array}{lll}
\frac{\partial F_{1}}{\partial \eta_{1}}=p_{x 1}(1), & \frac{\partial F_{1}}{\partial \eta_{2}}=p_{x 2}(1), & \frac{\partial F_{1}}{\partial L}=p_{x L}(1), \\
\frac{\partial F_{2}}{\partial \eta_{1}}=p_{y 1}(1), & \frac{\partial F_{2}}{\partial \eta_{2}}=p_{y 2}(1), & \frac{\partial F_{2}}{\partial L}=p_{y L}(1) .
\end{array}
\]

Приведенный подход, основанный на методе стрельбы и вариационных дифференциальных уравнениях, называется в литературе методом GPM $[6.8,6.13,6.14]$.

Диаграмма решений, найденная с помощью метода стрельбы и алгоритма продолжения DERPAR для задачи 11 в случае ГУ2, представлена на рис. 6.3 [6.15]. Здесь же приведены решения, полученные с помощью «суммирования» (профилей) решений (см. формулу (4.3.17)). Для некоторых ветвей указана устойчивость соответствующих решений. Из рисунка видно, что при больших $L$ существует значительное число различных стационарных решений. Попробуем оценить это число. Рассмотрим ветвь элементарных решений, т. е. решений, которые не могут быть получены путем «сложения» решений с меньшими $L$. Пусть такая ветвь существует при $L \in\left(L_{i}, L_{i}+\Delta L_{i}\right)$. Если мы выберем максимально широкие промежутки существования, то тогда $L_{i}$ и $L_{i}+\Delta L_{i}$ представляют собой координаты точек бифуркации. Фиксируем теперь длину $L$ и исследуем, сколько различных решений, полученных «сложением» из решений этой ветви, будет существовать для такого $L$. Это число равняется разности
\[
N_{i}=\left[\frac{L}{L_{i}}\right]-\left[\frac{L}{L_{i}+\Delta L_{i}}\right],
\]

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

длины $L$ будет равно
\[
N=1+\sum_{i} N_{i},
\]

где единица представляет собой (при всех $L$ ) существующее тривиальное решение $x \equiv A, y \equiv B / A$. Эти рассуждения остаются справедливыми для любых систем типа «реакция-диффузия» с граничными условиями типа ГУ2. Если нам известны не

Рис. 6.3. Диаграмма стационарных решений для задачи 11 при $Г У 2, A=2$, $B=4,6, D_{\mathrm{x}}=0,0016, D_{\mathrm{y}}=0,008, \eta_{1}=x(0) ; \mathrm{s}$ – устойчивые, $\mathrm{n}$ – неустойчивые решения.

все ветви элементарных решений, то тогда выражение (6.2.40) дает нижнюю оценку числа решений для данного $L$. Так, при $L=1$ на основании рис. 6.3 мы находим, что существует минимум 35 стационарных решений задачи 11 для случая ГУ2 (см. табл. 6.8 с учетом того, что существуют еще тривиальные решения задачи (при любых $L$ )).

Пример диаграммы стационарных решений для задачи 11 в случае ГУ 1 представлен на рис. 6.4 , а более полную картину решений наряду с профилями $x(z), y(z)$ читатель может найти в [6.15]. У некоторых ветвей на рис. 6.4 указан характер их

Таблица 6.8. Число $N_{i}$ стационарных решений задачи 11 для случая ГУ2 $(L=1)$, полученных $M_{i}$-кратным «сложением» элементарных решений при $L=L_{i}, A=2, B=4,6, D_{x}=0,0016, D_{y}=0,008$.

устойчивости. Отметим также, что операция «сложения решений» в случае ГУ1 оказывается невозможной.

Рис. 6.4. Диаграмма рещений для задачи 11 при ГУ $1, A=2, B=4,6, D_{\mathrm{x}}=$ $=0,0016, D_{\mathrm{y}}=0,008, \eta_{1}=x^{\prime}(0), \eta_{2}=y^{\prime}(0) ; \mathrm{s}-$ устойчивые, $\mathrm{n}-$ неустойчивые решения.
Часть диаграммы решений для задачи 12 в случае ГУ2, содержащая только элементарные решения, представлена на рис. 6.5. При выбранных параметрах задача имеет три тривиальных решения (они указаны в описании к рисунку, так что читатель может сравнить эти данные с данными рис. $5.6 \mathrm{~b}$ и $5.20 \mathrm{e}$ ). Точно так же, как и в случае рис. 6.3, мы можем под-

считать число решений при больших $L$, полученных путем сложения элементарных решений. Диаграммы решений с использованием условий типа ГУ1 представлены на рис. 6.6 для случая граничных условий, определяемых тривиальным решением № 3 , указанным в описании к рис. 6.5.

На рис. 6.7 изображена диаграмма решений задачи 13 в случае граничных условий ГУ2. Здесь же частично представлены решения, полученные сложением решений, и указан характер устойчивости отдельных ветвей.

Рис. 6.5. Диаграмма решений для задачи 12 при ГУ $2, \alpha=12, \beta=1,5, \gamma=$ $=3, \delta=1, v_{0}=0,01, D_{x}=0,008, D_{y}=0,004 ; \eta_{1}=x(0), \eta_{2}=y(0)$.
Тривиальные решения в системе без пространственных градиентов $\left(D_{\mathrm{x}}=\right.$ $\left.=D_{\mathrm{y}}=0\right)$ :

Покажем теперь, как изменится схема алгоритма в случае задачи 14. Будем искать зависимость решения уравнений (6.1.35a,b) с граничными условиями $(6.1 .36 \mathrm{a}, \mathrm{b})$ от параметра Da. Выберем два недостающих начальных условия в форме (6.1.37) и положим в точке $z=0$
\[
\begin{array}{l}
F_{1}(\boldsymbol{\eta}, \mathrm{Da})=\mathrm{Pe}_{\mathrm{M}} y(0)-y^{\prime}(0)=0, \\
F_{2}(\boldsymbol{\eta}, \mathrm{Da})=\mathrm{Pe}_{\mathrm{H}} \Theta(0)-\Theta^{\prime}(0)=0 .
\end{array}
\]

Вариационные уравнения для переменных $p_{y 1}=\partial y / \partial \eta_{1}, p_{y 2}=$ $=\partial y / \partial \eta_{2}, \quad p_{y \mathrm{Da}}=\partial y / \partial \mathrm{Da}, p_{\Theta 1}=\partial \Theta / \partial \eta_{1}, p_{\Theta 2}=\partial \Theta / \partial \eta_{2}, p_{\Theta \mathrm{Da}}=$ $=\partial \Theta / \partial \mathrm{Da}$ читатель может легко получить дифференцированием уравнений $(6.1 .35 \mathrm{a}, \mathrm{b})$ по $\eta_{1}, \eta_{2}$ и Da. Начальные условия
20 М. Холодннок и др.

Рис. 6.6. Диаграмма решений для задачи 12 при ГУ1. $\bar{x}=0,7006, \bar{y}=3,5106$; $\eta_{1}=x^{\prime}(0), \eta_{2}=y^{\prime}(0)$. Остальные параметры указаны на рис. $6.5 ; \mathrm{s}-$ устойчивые, $\mathrm{n}$ – неустойчивые решения.

для них выбираются нулевыми, кроме условий $p_{y_{1}}(1)=1$, $p_{\theta 2}(1)=1$. При этом матрица Якоби системы (6.2.41) будет иметь вид
\[
\left[\begin{array}{lll}
\mathrm{Pe}_{\mathrm{M}} p_{y 1}(0)-p_{y 1}^{\prime}(0), & \mathrm{Pe}_{\mathrm{M}} p_{y 2}(0)-p_{y 2}^{\prime}(0), \mathrm{Pe}_{M} p_{y \mathrm{Da}}(0)-p_{y \mathrm{Da}}^{\prime}(0) \\
\mathrm{Pe}_{\mathrm{H}} p_{\Theta 1}(0)-p_{\Theta 1}^{\prime}(0), \mathrm{Pe}_{\mathrm{H}} p_{\mathrm{\theta} 2}(0)-p_{\Theta 2}^{\prime}(0), & \mathrm{Pe}_{\mathrm{H}} p_{\Theta \mathrm{Da}}(0)-p_{\Theta \mathrm{Da}}^{\prime}(0)
\end{array}\right] .
\]

Диаграмма стационарных решений, полученная этим методом ${ }^{1)}$, приведена на рис. 6.8. При $\mathrm{Da}=0$ начальную точку зависимости можно найти аналитически. Сравните число решений для различных значений параметра $\mathrm{Da}$, приведенных в табл. 6.4 и 6.5 .

Рис. 6.7. Диаграмма решений для задачи 13 при ГУ2; $\mu=0,0035,
u=0,0045$, $\rho_{0}=6 \cdot 10^{-4}, c=0,05, c^{\prime}=0,025, \rho=\rho^{\prime}=3,2, D_{\mathrm{x}}=0,01, D_{\mathrm{y}}=0,45, \eta_{1}=$ $=x(0), \quad \eta_{2}=y(0) ;$ сплошные линии – устойчивые, прерывистые – неустойчивые решения.
В заключение этого пункта проиллюстрируем, как найти зависимость решения от параметра в случае задачи 15. Здесь годится подход, вполне аналогичный тому, который использовался нами в задаче 14 (в чем, впрочем, мы могли убедиться еще в предыдущем параграфе). Итак, будем искать зависимость решения уравнений (Р15-6), (P15-7) совместно с нелинейными (алгебраическими) уравнениями (Р15-8), (Р15-9) от параметpa Da. Граничные условия вновь имеют форму (6.1.36a, b). Два недостающих начальных условия в точке $z=1$ выбираем в
1) То есть применением алгоритма DERPAR к нахождению кривой $F_{1}\left(\eta_{1}, \eta_{2}, \mathrm{Da}\right)=F_{2}\left(\eta_{1}, \eta_{2}, \mathrm{Da}\right)=0 .-\Pi$ рим. ред.
$9 *$

Рис. 6.8. Диаграмма решений для задачи 14. $\mathrm{P}_{\mathrm{M}}=10, \mathrm{Pe}_{\mathrm{H}}=5, B=15$ 。 $\boldsymbol{\beta}=2, \mathbf{\theta}_{c}=0, \gamma=20, \eta_{1}=y(1), \eta_{2}=\mathbf{\theta}(1)$.

Рис. 6.9. Диаграмма решений задачи 15. Зависимость температуры на выходе $\eta_{2}=\boldsymbol{\theta}(1)$ от параметра $\mathrm{D}$. $\mathrm{P}_{\mathrm{M}}=20, \mathrm{Pe}_{\mathrm{H}}=10, \beta=1, B=10, \gamma=20$, $\boldsymbol{\theta}_{c}=0, J_{\mathrm{M}}=J_{\mathrm{H}}=25$.

соответствии с формулами (6.1.37), а в точке $z=0$ имеем два нелинейных уравнения вида (6.2.41). Вариационное уравнение для $p_{y 1}$ имеет вид (6.1.40). Вариационные уравнения для остальных переменных (6.1.39), а также для переменных $p_{y \text { Da }}=\partial y / \partial \mathrm{Da}$, $p_{\Theta \mathrm{Da}}=\partial \Theta / \partial \mathrm{Da}$ выводятся аналогично. При этом матрица системы (6.2.41) имеет вид (6.2.42).

Напомним, что на каждом шаге интегрирования уравнений (P15-6), (P15-7), а также соответствующих вариационных уравнений нам необходимо вычислять значения переменных (6.1.41), т. е. переменных $\partial \omega / \partial y, \partial \omega / \partial \Theta, \partial \theta / \partial y, \partial \theta / \partial \Theta$. Эти значения мы находим по формулам (6.1.43) так же, как это делалось в $\S 6.1$.

На рис. 6.9 приведена диаграмма решений рассматриваемой задачи в зависимости от параметра Da. Қак видно из рисунка, в узком интервале значений параметра данная задача имеет 5 решений. Для указанных значений параметров уравнение (P15-11) имело только одно решение $\omega$ в промежутке $[y(z), 1]$ при любом $z \in[0,1]$.

Построение диаграммы решений несколько осложняется в тех случаях, когда уравнение (Р15-11) имеет несколько корней. Подробный анализ задачи 15 читатель может найти в сборнике [6.16].

6.2.4. Метод отображения параметра

В некоторых случаях оказывается возможным достаточно быстро построить диаграмму стационарных решений, если к выбранному значению некоторого начального условия с помощью интегрирования задачи Коши мы добавим соответствующее значение параметра. Аналогичный подход был описан для «сосредоточенных» систем в п. 5.2.1. Так же, как и там, возможность применения этого подхода зависит от конкретного вида дифференциальных уравнений и граничных условий, а также от того, как входит в уравнения выбранный параметр. Поэтому мы продемонстрируем применение указанного метода на двух конкретных задачах из гл. 4. Более общие соображения читатель. может найти в книгах $[6.8,6.9]$ и в приведенной в них библиографии.

Будем строить зависимость решения задачи 16 от параметра Ф. Для стационарного случая мы свели эту задачу к одному дифференциальному уравнению второго порядка (6.1.21) с граничными условиями ( $6.1 .16 \mathrm{a}, \mathrm{c})$. Введем новую независимую переменную
\[
z=\Phi r
\]

Тогда дифференциальное уравнение (6.1.21) приведется к виду
\[
\frac{d^{2} y}{d z^{2}}+\frac{a}{z} \frac{d y}{d z}-y \exp \frac{\gamma \beta(1-y)}{1+\beta(1-y)}=0,
\]
т. е. не будет содержать Ф.
Граничные условия приобретут вид
\[
\begin{array}{l}
z=0: \frac{d y}{d z}=0, \\
z=\Phi: y=1 .
\end{array}
\]

Выберем теперь недостающее начальное условие в виде
\[
y(0)=\eta, \quad 0<\eta<1
\]

и проинтегрируем полученную задачу Коши (6.2.44), (6.2.45а), (6.2.46) от $z=0$ до точки $z=z_{1}$, в которой
\[
y\left(z_{1}\right)=1 \text {. }
\]

Значение $z_{1}$ определяет значение параметра $\Phi$, соответствующее выбранному начальному условию $\eta$ (и найденному решению $y(z))$ :
\[
z_{1}=\Phi \text {. }
\]

Для выполнения равенства (6.2.47) можно воспользоваться каким-либо методом последовательных приближений, например методом деления промежутка пополам. Процесс можно, например, реализовать так. Когда в процессе интегрирования становится $y(z)>1$, мы возвращаемся на один шаг назад и продолжаем интегрирование с более коротким шагом (например, уменьшенным вдвое). По достижении достаточно короткого шага процесс прекращается.

Опишем коротко неитерационный способ нахождения величины $z_{1}$. Проинтегрируем уравнение (6.2.44), переписав его в виде системы двух дифференциальных уравнений первого порядка,
\[
\frac{d y}{d z}=w, \quad \frac{d w}{d z}=-\frac{a}{z} w+y \exp \frac{\gamma \beta(1-y)}{1+\beta(1-y)},
\]

с начальными условиями (6.2.45a), (6.2.46), т. е. с условиями $y(0)=\eta, w(0)=0$. В некоторой точке $\bar{z}>0$, где $y(\bar{z})=\bar{y}<1$ и $d y / d z=\bar{y}^{\prime}$, мы переходим к интегрированию дифференциальных уравнений вида
\[
\frac{d w}{d y}=-\frac{a}{z}+\frac{y}{w} \exp \frac{\gamma \beta(1-y)}{1+\beta(1-y)}, \quad \frac{d z}{d y}=\frac{1}{w}
\]

от $y=\bar{y}$, где $w(\bar{y})=\bar{y}^{\prime}, z(\bar{y})=\bar{z}$ до $y=1$. После этого полагаем $z(1)=\Phi$ (ср. аналогичный процесс вычисления отображения Пуанкаре, п. 5.9.2). Полученные таким образом значения $\Phi$ для последовательности значений $\eta$ приведены в табл. 6.9. Из этих данных видно, что при $\beta=0,4$ для определенного интервала значений $\Phi$ (например для $\Phi=0,3$ ) существует три решения исходной краевой задачи. Читатель может сравнить результаты, полученные при $\beta=0,05$, с результатами, приведенными в табл. 6.1.

Таблица 6.9. Результаты применения метода отображения параметра для задачи $16(\gamma=20, a=0)$. Значения Ф.

Другим примером, на котором мы рассмотрим возможности метода отображения параметра, является система типа «реакция-диффузия» (4.3.16) при $D_{\mathrm{y}} \rightarrow 0$. В этом предельном случае второе уравнение системы (4.3.16) сводится к конечному уравнению вида
\[
g(x, y)=0,
\]

решая которое, мы находим зависимость $y(x)$ (часто ваналитическом виде). Подстановка в уравнение (6.1.1) дает
\[
x^{\prime \prime}=-\frac{L^{2}}{D_{\mathrm{x}}} f(x, y(x))=-\frac{L^{2}}{D_{\mathrm{x}}} \bar{f}(x) .
\]

Положив
\[
\xi=L z,
\]

получим уравнение, не содержащее параметра $L$ :
\[
\frac{d^{2} x}{d \xi^{2}}=-\frac{1}{D_{\mathrm{x}}} \bar{f}(x) .
\]

Рассмотрим для примера граничные условия 1 рода (4.3.9). Для уравнения (6.2.54) они принимают вид
\[
\xi=0: x=\bar{x} ; \quad \xi=L: x=\bar{x} .
\]

Зададим дополнительно при $\xi=0$
\[
\frac{d x(0)}{d \xi}=\eta,
\]

после чего проинтегрируем уравнение (6.2.54) от $\xi=0$, где $x(0)=\bar{x}, x^{\prime}(0)=\eta$, до такого значения $\xi=L$, где $x(L)=\bar{x}$. Это значение $L$ мы и сопоставим выбранному $\eta$.