Для иллюстрации методов нахождения стационарных решений параболических дифференциальных уравнений с частными производными рассмотрим системы типа «реакция-диффузия» (см. задачи 11-13). Левые части уравнений (4.3.7) при этом мы полагаем равными нулю. Таким образом, стационарное решение удовлетворяет уравнениям (‘ $=d / d z$ )
\[
\begin{array}{l}
x^{\prime \prime}+\frac{L^{2}}{D_{x}} f(x, y)=0, \\
y^{\prime \prime}+\frac{L^{2}}{D_{y}} g(x, y)=0 .
\end{array}
\]

Запишем граничные условия (4.3.8), (4.3.12), (4.3.13) (для граничных условий первого рода мы рассмотрим только симметричный случай (4.3.9)):
ГУ1:
\[
\begin{array}{ll}
x(0)=\bar{x}, & y(0)=\bar{y}, \\
x(1)=\bar{x}, & y(1)=\bar{y} .
\end{array}
\]

ГУ2:
\[
x^{\prime}(0)=0, \quad y^{\prime}(0)=0,
\]
\[
x^{\prime}(1)=0, \quad y^{\prime}(1)=0 .
\]

ГУ 3: $\quad \alpha_{\mathrm{x} 0} x^{\prime}(0)+\beta_{\mathrm{x} 0} x(0)=\gamma_{\mathrm{x} 0}, \alpha_{\mathrm{y} 0} y^{\prime}(0)+\beta_{\mathrm{y} 0} y(0)=\gamma_{\mathrm{y} 0}$,
\[
\alpha_{\mathrm{x} 1} x^{\prime}(1)+\beta_{\mathrm{x} 1} x(1)=\gamma_{\mathrm{x} 1}, \alpha_{\mathrm{y} 1} y^{\prime}(1)+\beta_{\mathrm{y} 1} y(1)=\gamma_{\mathrm{y} 1} .
\]

Если для граничных условий первого рода $f(\bar{x}, \bar{y})=g(\bar{x}, \bar{y})=0$, то имеется однородное по пространству решение системы $(6.1 .1-2): x(z) \equiv \bar{x}, y(z) \equiv \bar{y}$. Это решение, очевидно, удовлетворяет и граничным условиям второго рода. Численными решениями нелинейных краевых задач занимался целый ряд авторов, среди публикаций которых можно найти самые разные работы – от чисто теоретических статей до сугубо прикладных исследований. Здесь мы рассмотрим указанную проблему сравнительно кратко; читателей же, которые заинтересуются этой проблемой более глубоко, мы отсылаем к монографической литературе $[6.5,6.6,6.7,6.8,6.9,6.10]$.

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

6.1.1. Разностные методы

На промежутке $z \in[0,1]$ выберем сетку узловых точек $z_{0}=0, z_{1}, \ldots, z_{n}=1$. Сетка, как правило, выбирается эквидистантной, $z_{i}=i h, i=0,1, \ldots, n$, где $h=1 / n$ есть шаг сетки. На этой сетке узлов величины $x_{i}, y_{i}$ мы рассматриваем как аппроксимации значений решения $x\left(z_{i}\right), y\left(z_{i}\right), i=0,1, \ldots, n$. Далее, производные в уравнениях (6.1.1) и (6.1.2) заменяются соответствующими разностными формулами, например, трехточечными центральными разностями:
\[
\begin{array}{c}
x_{i}^{\prime \prime} \sim x^{\prime \prime}\left(z_{i}\right) \sim \frac{x_{i-1}-2 x_{i}+x_{i+1}}{h^{2}}, \\
y_{i}^{\prime \prime} \sim y^{\prime \prime}\left(z_{i}\right) \sim \frac{y_{i-1}-2 y_{i}+y_{i+1}}{h^{2}} .
\end{array}
\]

Подстановка этих выражений в уравнения (6.1.1) и (6.1.2) приводит к следующей системе нелинейных уравнений:
\[
\begin{array}{ll}
x_{i-1}-2 x_{i}+x_{i+1}+\frac{L^{2} h^{2}}{D_{x}} f\left(x_{i}, y_{i}\right)=0, & i=1,2, \ldots, n-1, \\
y_{i-1}-2 y_{i}+y_{i+1}+\frac{L^{2} h^{2}}{D_{y}} g\left(x_{i}, y_{i}\right)=0, & i=1,2, \ldots, n-1 .
\end{array}
\]

Добавляя к ним аппроксимации граничных условий, мы получаем систему $2 n+2$ уравнений. Так, для случая ГУ1 имеем
\[
\begin{array}{ll}
x_{0}=\bar{x}, & y_{0}=\bar{y}, \\
x_{n}=\bar{x}, & y_{n}=\bar{y} .
\end{array}
\]

Для ГУ2 (или ГУЗ) мы должны, кроме того, заменить производные в крайних точках. Если для аппроксимации первой производной в ГУ 2 использовать двухточечную замену типа
\[
x_{0}^{\prime} \sim \frac{x_{1}-x_{0}}{h}=0,
\]

то порядок аппроксимации в формулах (6.1.7), (6.1.8) понизится с $O\left(h^{2}\right)$ до $O(h)$. Одна из возможностей сохранить порядок аппроксимации равным $O\left(h^{2}\right)$ – использовать несимметричную разностную замену первой производной, включающую три узловые точки:
\[
x_{0}^{\prime} \sim \frac{-3 x_{0}+4 x_{1}-x_{2}}{2 h} .
\]

Другая возможность состоит в использовании аппроксимации, основанной на виртуальных точках с индексами $i=-1$ и $i=$ $=n+1$. Именно, распространим соотношения (6.1.7) и (6.1.8) на случай индексов $i=0$ и $i=n$ и заменим граничные условия (6.1.4) трехточечными центральными разностями:
\[
\begin{array}{l}
x_{0}^{\prime} \sim \frac{x_{1}-x_{-1}}{2 h}=0, \quad y_{0}^{\prime} \sim \frac{y_{1}-y_{-1}}{2 h}=0, \\
x_{n}^{\prime} \sim \frac{x_{n+1}-x_{n-1}}{2 h}=0, \quad y_{n}^{\prime} \sim \frac{y_{n+1}-y_{n-1}}{2 h}=0 .
\end{array}
\]

Подставляя значения $x_{-1}, y_{-1}, x_{n+1}$ и $y_{n+1}$ из этих соотношений в формулы (6.1.7) и (6.1.8), находим
\[
\begin{array}{l}
2 x_{1}-2 x_{0}+\frac{L^{2} h^{2}}{D_{x}} f\left(x_{0}, y_{0}\right)=0, \\
2 y_{1}-2 y_{0}+\frac{L^{2} h^{2}}{D_{y}} g\left(x_{0}, y_{0}\right)=0, \\
2 x_{n-1}-2 x_{n}+\frac{L^{2} h^{2}}{D_{x}} f\left(x_{n}, y_{n}\right)=0, \\
2 y_{n-1}-2 y_{n}+\frac{L^{2} h^{2}}{D_{y}} g\left(x_{n}, y_{n}\right)=0 .
\end{array}
\]

Так же, как и в случае ГУ 1 , мы получили систему $2 n+2$ нелинейных уравнений (6.1.7), (6.1.8), (6.1.11), (6.1.12) относительно $2 n+2$ неизвестных. Упорядочим эти неизвестные

следующим образом:
\[
\mathbf{X}=\left(x_{0}, y_{0}, x_{1}, y_{1}, x_{2}, \ldots, x_{n-1}, y_{n-1}, x_{n}, y_{n}\right) .
\]

Перепишем теперь соответствующие уравнения для случая ГУ2 в виде последовательности:
1. (6.1.11a),
2. $(6.1 .11 \mathrm{~b})$,
3. (6.1.7) при $i=1$,
4. (6.1.8) при $i=1$,
5. (6.1.7) при $i=2$,
6. (6.1.8) при $i=2$,
\[
\begin{array}{l}
2 n-1 .(6.1 .7) \text { при } i=n-1, \\
2 n . \quad(6.1 .8) \text { при } i=n-1, \\
2 n+1 .(6.1 .12 \mathrm{a}), \\
2 n+2 .(6.1 .12 \mathrm{~b}) .
\end{array}
\]

Мы получили систему нелинейных уравнений со специальной структурой вхождения неизвестных в уравнениях. Для описания этой структуры вводится специальная схема (матрица) размещения, имеющая столько строк, сколько исходная система уравнений, и столько столбцов, сколько неизвестных имеется в системе. Элементами матрицы размещения служат либо нули (неизвестная в соответствующем уравнении не фигурирует), либо крестики (неизвестная входит в соответствующее уравнение). Таким образом, первые шесть строк матрицы размещения для системы (6.1.13) – (6.1.14) имеют вид
\[
\begin{array}{l}
\mathrm{x} \times \mathrm{x} 0 \text {. . . . . . . . } 0 \\
\mathrm{x} \times 0 \mathrm{x} 0 \text {. . . . . } 0 \\
\mathrm{x} 0 \mathrm{x} \times \mathrm{x} 0 \text {….. } 0 \\
0 \times \times \times 0 \times 0 \text {…. } 0 \\
00 \times 0 \times \times \times 0 \text {. . . } 0 \\
000 \times \times \times 0 \times 0 \ldots 0 \\
\end{array}
\]

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

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

Применим теперь метод конечных разностей для нахождения стационарного решения задачи 16. Это стационарное решение описывается соотношениями (см. уравнения (P16-11), (P16-12))
\[
\begin{array}{l}
y^{\prime \prime}+\frac{a}{r} y^{\prime}-\Phi^{2} y \exp \frac{\Theta}{1+\Theta / \gamma}=0, \\
\Theta^{\prime \prime}+\frac{a}{r} \Theta^{\prime}+\gamma \beta \Phi^{2} y \exp \frac{\Theta}{1+\Theta / \gamma}=0,
\end{array}
\]

а граничные условия (P16-13), (P16-14) принимают вид
\[
\begin{array}{l}
y^{\prime}(0)=0, \\
\Theta^{\prime}(0)=0, \\
y(1)+\frac{1}{\mathrm{Sh}} y^{\prime}(1)=1, \\
\Theta(1)+\frac{1}{\mathrm{Nu}} \Theta^{\prime}(1)=0 .
\end{array}
\]

Иногда при вычислении стационарного решения удается понизить размерность исходной задачи. Продемонстрируем эту возможность на данном примере, предполагая $\mathrm{Nu}=\mathrm{Sh}$. Умножая уравнение (6.1.15a) на множитель $\gamma \beta$, складывая результат с уравнением (6.1.15b) и вводя обозначение
\[
u=\gamma \beta y+\Theta,
\]

мы получаем уравнение
\[
u^{\prime \prime}+\frac{a}{r} u^{\prime}=0
\]

Граничные условия (6.1.16a,b, c, d) переписываются в виде
\[
\begin{array}{c}
u^{\prime}(0)=0, \\
u(1)+\frac{1}{\mathrm{Nu}} u^{\prime}(1)=\gamma \beta .
\end{array}
\]

Левая часть уравнения (6.1.18) равняется $\left[r^{a} u^{\prime}\right]^{\prime} / r^{a}$, откуда с использованием (6.1.19a) находим $u^{\prime}=0$. Интегрируя, имеем $u=C$, и из условия (6.1.19b) следует, что $u(r) \equiv \gamma \beta$. Теперь из формулы (6.1.17) вытекает зависимость между $\Theta$ и $y$
\[
\Theta=\gamma \beta(1-y) \text {. }
\]

Заметим, что при $\mathrm{Nu}
eq \mathrm{Sh}$ это соотношение оказывается несправедливым. Таким образом, система (6.1.15) – (6.1.16) сводится к краевой задаче для одного дифференциального уравнения второго порядка
\[
y^{\prime \prime}+\frac{a}{r} y^{\prime}-\Phi^{2} y \exp \frac{\gamma \beta(1-y)}{1+\beta(1-y)}=0
\]

с граничными условиями $(6.1 .16 \mathrm{a}, \mathrm{c})$. Рассмотрим далее для простоты предельный случай $\mathrm{Nu} \rightarrow \infty, \mathrm{Sh} \rightarrow \infty$ (см. гл. 4, задача 16). Выберем сетку узловых точек следующим образом: $r_{i}=i h, i=0,1, \ldots, n, h=1 / n$. Тогда разностный аналог уравнения (6.1.21) принимает вид
\[
\begin{array}{c}
\frac{y_{i-1}-2 y_{i}+y_{i+1}}{h^{2}}+\frac{a}{i h} \frac{y_{i+1}-y_{i-1}}{2 h}-\Phi^{2} y_{i} \exp \frac{\gamma \beta\left(1-y_{i}\right)}{1+\beta\left(1-y_{i}\right)}=0, \\
i=1,2, \ldots, n-1 .
\end{array}
\]

Разностный аналог граничного условия (6.1.16a) получим, введя виртуальную точку $r_{-1}$. В случае $r=0$, однако, в уравнении (6.1.21) имеется неопределенное выражение $(a / r) y^{\prime}$ типа $0 / 0$. При $r \rightarrow 0$ находим
\[
y^{\prime \prime}+(a / r) y^{\prime} \rightarrow(a+1) y^{\prime \prime} .
\]

Далее, используя разностную замену для уравнения (6.1.21) при $r=0$, с учетом требования $y_{-1}=y_{1}$ (которое вытекает из условия (6.1.16a)) получаем окончательно
\[
\frac{a+1}{h^{2}}\left(2 y_{1}-2 y_{0}\right)-\Phi^{2} y_{0} \exp \frac{\gamma \beta\left(1-y_{0}\right)}{1+\boldsymbol{\beta}\left(1-y_{0}\right)}=0 .
\]

Наконец, граничное условие (6.1.16c) при $\mathrm{Sh} \rightarrow \infty$ аппроксимируется с помощью соотношения
\[
y_{n}-1=0 .
\]

Читатель может легко построить матрицу Якоби $\mathbf{G}=\left[g_{i j}\right]$, $i, j=0,1, \ldots, n$, для решения системы методом Ньютона. Эта матрица будет трехдиагональной, а ее элементы (производные левой части соотношения (6.1.22a)) имеют вид
\[
\begin{array}{c}
g_{i, i-1}=\frac{1}{h^{2}}-\frac{a}{2 i h^{2}}, \quad i=1,2, \ldots, n-1, \\
g_{i, i+1}=\frac{1}{h^{2}}+\frac{a}{2 i h^{2}}, \quad i=1,2, \ldots, n-1, \\
g_{i, i}=-\frac{2}{h^{2}}-\Phi^{2}\left[1+\frac{\gamma \beta y_{i}}{\left(1+\beta\left(1-y_{i}\right)\right)^{2}}\right] \exp \frac{\gamma \beta\left(1-y_{i}\right)}{1+\beta\left(1-y_{i}\right)}, \\
i=1,2, \ldots, n-1 .
\end{array}
\]

Соответствующие выражения для $g_{00}, g_{01}, g_{n n}$ нетрудно получить дифференцированием уравнений $(6.1 .22 \mathrm{~b}, \mathrm{c})$. Ход итерационного процесса в случае применения метода Ньютона представлен в табл. 6.1. Точное значение $y(0)$ в этой задаче равно 0,5521 , так что погрешность аппроксимации при $h=0,1$ сказывается лишь в четвертом знаке. Заметим, что приведенный выше пример оказывается настолько простым, что его вполне можно анализировать с помощью небольшой персональной ЭВМ.

В качестве более сложного примера рассмотрим решение з $a$ дачи 17. Соответствующие разностные аналоги уравнений $(\mathrm{P} 17-16)-(\mathrm{P} 17-20)$ имеют вид
\[
H_{i}-H_{i-1}+h \sqrt{\operatorname{Re}}\left(F_{i}+F_{i-1}\right)=0,
\]
\[
\begin{array}{c}
F_{i-1}-2 F_{i}+F_{i+1}- \\
-\frac{1}{2} h \sqrt{\operatorname{Re}} H_{i}\left(F_{i+1}-F_{i-1}\right)- \\
\quad-h^{2} \operatorname{Re}\left(F_{i}^{2}-G_{i}^{2}+k\right)=0, \\
G_{i-1}-2 G_{i}+G_{i+1}-2 h^{2} \operatorname{Re} F_{i} G_{i}- \\
-\frac{1}{2} h \sqrt{\operatorname{Re}}\left(G_{i+1}-G_{i-1}\right) H_{i}=0, \\
i=1,2, \ldots, n-1, \quad(6.1 .23 \mathrm{~b})
\end{array}
\]
\[
H_{n}-H_{n-1}+h \sqrt{\operatorname{Re}}\left(F_{n}+F_{n-1}\right)=0 \text {, }
\]
\[
\begin{array}{l}
H_{n}=0, \\
F_{n}=0, \\
G_{n}=: S .
\end{array}
\]

Допустим, что значения параметров $\operatorname{Re}$ и $S$ заданы. Упорядочим неизвестные следующим образом: $H_{0}, F_{0}, G_{0}, H_{1}, F_{1}$, $G_{1}, \ldots, H_{n-1}, F_{n-1}, G_{n-1}, H_{n}, F_{n}, G_{n}, k$. Мы получим тогда семидиагональную матрицу размещения с полностью заполненным последним столбцом (что соответствует появлению неизвестной $k$ во всех уравнениях (6.1.23b)). В этом случае общее число уравнений оказывается равным $3(n-1)+7$. При интересных для этой задачи числах Рейнольдса ( $\mathrm{Re} \sim 500 \div 1000$ ) нужно довольно много узлов. Например, для $\mathrm{Re}=625$ необходимо $n=100 \div 200$ [6.11]. Подобную задачу следует решать уже на достаточно мощных ЭВМ, даже если использовать специальную программу для решения систем линейных алгебраических уравнений с почти семидиагональной матрицей. Кроме того, в данной задаче при достаточно больших значениях параметра $\mathrm{Re}$ ( $\mathrm{Re}>700$ ) появляются «паразитные» решения. Так, при $n=100$ и $n=200$ существуют решения, которые в случае более мелкого разбиения, например при $n=800$, исчезают (см. [6.11]). Учитывая это обстоятельство, при использовании разностных методов представляется необходимым результаты, полученные при достаточно грубом разбиении, пересчитывать на более мелкой сетке узловых точек. При этом только хорошее совпадение результатов для нескольких последовательных дроблений может служить критерием правильности решения.

6.1.2. Метод стрельбы

Основная идея метода стрельбы заключается в сведении краевой задачи к задаче Қоши, для решения которой можно использовать стандартные программные средства, имеющиеся практически на каждой ЭВМ (см. § 5.7). Ниже мы рассмотрим простейшие варианты метода стрельбы, вполне достаточные для понимания его концептуальной стороны. Более подробное изложение этого метода читатель может найти в учебной литературе $[6.7,6.8]$.

Опишем подробнее метод стрельбы на примере задачи (6.1.1), (6.1.2) с граничными условиями (6.1.4). Для того чтобы данную задачу четвертого порядка можно было решать как задачу Коши, мы должны задать в некоторой точке четыре начальных условия. Воспользуемся тем, что в точке $z=0$ уже имеются два заданных условия, а именно условия (6.1.4a), и выберем два дополнительных условия вида
\[
x(0)=\eta_{1}, \quad y(0)=\eta_{2} .
\]

Полученную задачу Коши на промежутке от $z=0$ до $z=1$ можно решать с помощью какого-либо из методов, рассмотрен-

ных в §5.7. (Эти методы предназначались для решения систем дифференциальных уравнений первого порядка. Уравнения (6.1.1), (6.1.2) легко преобразуются в систему четырех уравнений первого порядка: $x^{\prime}=u, \quad u^{\prime}=-L^{2} f(x, y) / D_{x}, y^{\prime}=v$, $v^{\prime}=-L^{2} g(x, y) / D_{\text {y }}$. Начальные условия (6.1.4а) и (6.1.24) преобразуются к виду $u(0)=v(0)=0, x(0)=\eta_{1}, y(0)=\eta_{2}$. Мы будем, далее придерживаться исходных обозначений).

По окончании интегрирования, т. е. в точке $z=1$, мы получаем значения решения, зависящие от выбора условий (6.1.24): $x(1, \boldsymbol{\eta}), x^{\prime}(1, \boldsymbol{\eta}), y(1, \boldsymbol{\eta}), y^{\prime}(1, \boldsymbol{\eta})$. Для того чтобы найденное решение задачи Коши было одновременно и решением исходной краевой задачи, нам необходимо удовлетворить условиям (6.1.4b), т. е.
\[
F_{1}(\boldsymbol{\eta})=x^{\prime}(1, \boldsymbol{\eta})=0, \quad F_{2}(\boldsymbol{\eta})=y^{\prime}(1, \boldsymbol{\eta})=0 .
\]

Соотношения (6.1.25) представляют собой два нелинейных уравнения относительно неизвестных $\eta_{1}$ и $\eta_{2}$, которое можно решить, используя любой из методов решения нелинейных уравнений, например, метод Ньютона (см. § 5.1). Для вычисления элементов матрицы Якоби, т. е. производных $\partial F_{i} / \partial \eta_{j}$, можно воспользоваться соответствующими разностными заменами, вычислив функщии $F_{1}$ и $F_{2}$ в точках $\left(\eta_{1}, \eta_{2}\right),\left(\eta_{1}+\Delta \eta, \eta_{2}\right)$, $\left(\eta_{1}, \eta_{2}+\Delta \eta\right)$. Каждое такое вычисление требует решения задачи Коши (6.1.1), (6.1.2), (6.1.4a), (6.1.24) при соответствующем задании начальных условий. Погрешности, возникающие в процессе интегрирования, могут при этом перекрываться погрешностями аппроксимации использованных разностных формул, и элементы матрицы Якоби могут вычисляться весьма неточно. Поэтому часто, когда это возможно, используются вариационные дифференциальные уравнения (уравнения в вариациях) относительно переменных
\[
p_{x i}=\partial x / \partial \eta_{i}, \quad p_{y i}=\partial y / \partial \eta_{i}, \quad i=1,2 .
\]

Вариационные уравнения получаются посредством дифференцирования исходных уравнений ${ }^{1)}(6.1 .1)$, (6.1.2) по $\eta_{1}$ и $\eta_{2}$ и перестановки дифференцирования по $z$ и $\eta_{i}$ :
\[
\begin{array}{ll}
p_{x i}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{x}}}\left(\frac{\partial f}{\partial x} p_{x i}+\frac{\partial f}{\partial y} p_{y i}\right)=0, & i=1,2, \\
p_{y i}^{\prime \prime}+\frac{L^{2}}{D_{\mathrm{y}}}\left(\frac{\partial g}{\partial x} p_{x i}+\frac{\partial g}{\partial y} p_{y i}\right)=0, & i=1,2 .
\end{array}
\]
1) Более точно: дифференцированием тождеств, возникающих при подстановке решения $x(z, \eta), y(z, \eta)$ в уравнения (6.1.1), (6.1.2).

Начальные условия для варьируемых переменных получаются в результате дифференцирования исходных начальных условий (6.1.4а), (6.1.24) по $\eta_{1}$ и $\eta_{2}$
\[
\begin{array}{c}
p_{x 1}(0)=1, \quad p_{x 2}(0)=0, \quad p_{y 1}(0)=0, \quad p_{y 2}(0)=1, \\
p_{x i}^{\prime}(0)=p_{y i}^{\prime}(0)=0, \quad i=1,2 .
\end{array}
\]

Если проинтегрировать дифференциальное уравнение (6.1.27) вместе с уравнениями (6.1.1) и (6.1.2) и начальными условиями (6.1.4a), (6.1.24), (6.1.28), что в точке $z=1$ мы получим
\[
\begin{array}{ll}
\frac{\partial F_{1}}{\partial \eta_{1}}=p_{x 1}^{\prime}(1), & \frac{\partial F_{1}}{\partial \eta_{2}}=p_{x 2}^{\prime}(1), \\
\frac{\partial F_{2}}{\partial \eta_{1}}=p_{y 1}^{\prime}(1), & \frac{\partial F_{2}}{\partial \eta_{2}}=p_{y 2}^{\prime}(1) .
\end{array}
\]

Тем самым, мы имеем матрицу Якоби для применения метода Ньютона, вычисленную на основе интегрирования вариационных уравнений. Заметим, что в данном случае для одной итерации метода Ньютона нам приходится интегрировать систему 12 дифференциальных уравнений первого порядка. При аппроксимации элементов матрицы Якоби соответствующими разностными формулами нужно трижды интегрировать систему 4-х дифференциальных уравнений первого порядка. При этом объем вычислений для большой группы различных функций $f$ и $g$ оказывается приблизительно одинаковым. Иногда, правда, удается существенно сократить затраты машинного времени при переходе к вариационным уравнениям, например, при появлении функций exp, sin, cos, дифференцирование которых дает те же функции.

В случае ГУ1 дополнительные начальные условия к (6.1.3а) выбираются следующим образом:
\[
x^{\prime}(0)=\eta_{1}, \quad y^{\prime}(0)=\eta_{2} .
\]

При этом условия, которым нужно удовлетворить после интегрирования от $z=0$ до $z=1$, имеют вид (см. (6.1.3b)
\[
F_{1}(\boldsymbol{\eta})=x(1)-\bar{x}=0, \quad F_{2}(\boldsymbol{\eta})=y(1)-\bar{y}=0 .
\]

Вариационные уравнения (6.1.27) остаются теми же самыми, а начальные условия для варьируемых переменных заменяются на условия
\[
\begin{array}{c}
p_{x i}(0)=p_{y i}(0)=0, \quad i=1,2, \\
p_{x 1}^{\prime}(0)=1, \quad p_{x 2}^{\prime}(0)=0, \quad p_{y 1}^{\prime}(0)=0, \quad p_{y 2}^{\prime}(0)=1 .
\end{array}
\]

Для ГУЗ в форме (6.1.5) с ненулевыми коэффициентами (с тем чтобы в точке $z=0$ эти условия не совпали в ГУ1 или ГУ2)

дополнительные начальные условия в этой точке можно выбрать в виде (6.1.24) или (6.1.30). В случае выбора условий (6.1.24) из (6.1.5a) следует
\[
x^{\prime}(0)=\left(\gamma_{x 0}-\beta_{x 0} \eta_{1}\right) / \alpha_{x 0}, \quad y^{\prime}(0)=\left(\gamma_{y 0}-\beta_{y 0} \eta_{2}\right) / \alpha_{y 0} .
\]

При этом вариационные дифференциальные уравнения остаются теми же, а начальные условия заменяются следующими:
\[
\begin{array}{l}
p_{x 1}(0)=1, \quad p_{x 2}(0)=0, \quad p_{y 1^{1}}^{r}(0)=0, \quad p_{y 2}(0)=1, \\
\left.\left.p_{x 1}^{\prime}(0)=-\beta_{x 0}\right) / \alpha_{x 0}, \quad p_{x 2}^{\prime}(0)=0, \quad p_{y 1}^{\prime}(0)=0, \quad p_{y 2}^{\prime}(0)=-\beta_{y 0}\right) / \alpha_{y 0} . \\
\end{array}
\]

Конкретный вид функций $F_{i}$ и производных $\partial F_{i} / \partial \eta_{i}$ читатель может легко получить с помощью формул (6.1.5b).

Конечно, мы могли выбирать недостающие начальные условия в точке $z=1$ взамен точки $z=0$ и находить решения соответствующих задач Коши на промежутке от $z=1$ до $z=0$. В некоторых задачах выбор направления интегрирования может играть существенную роль, поскольку иногда численное интегрирование в одном из направлений оказывается труднореализуемым (соответствующая задача Коши неустойчива по отношению к начальным условиям). Это имеет место, например, в задачах 14 и 15 при больших значениях критерия Пекле, а также в задаче 16. Если задача Коши неустойчива в обоих направлениях, метод стрельбы применять нельзя и следует использовать, например, разностные методы решения, описанные в подпункте 6.1.1.

Рассмотрим теперь примеры расчетов методом стрельбы, иллюстрирующие процесс построения стационарных решений некоторых задач из гл. 4. В качестве первой из этих задач исследуем задачу 11, т. е. систему типа «реакция-диффузия» (6.1.1), (6.1.2), функции $f$ и $g$ для которой имеют вид (Р11-1). В табл. 6.2 приведены некоторые результаты для ГУ 2 при выборе начальных условий типа (6.1.24). Отметим, что второе из представленных в таблице решений однородно по пространству. Первое и третье решения зависят от $z$ и получаются друг из друга симметрией: $x^{\mathrm{I}}(z)=x^{\mathrm{III}}(1-z), \quad y^{\mathrm{I}}(z)=y^{\mathrm{III}}(1-z)$. В табл. 6.3 представлены соответствующие результаты для ГУ1, при выборе начальных условий типа (6.1.30). Здесь третье решение пространственно однородно: $x(z) \equiv 2 ; y(z) \equiv 2,3$. В обеих этих ситуациях применялся метод Ньютона, причем матрица Якоби подсчитывалась с помощью соответствующих уравнений в вариациях.

Таблица 6.2. Метод стрельбы для задачи 11 в случае ГУ2 $(A=2, B=4,6$, $D_{x}=0,0016, D_{y}=0,008, L=0,12$ ).

Таблица 6.3. Метод стрельбы для задачи 11 в случае ГУ 1 ( $A=2 ; B=4,6$; $\left.D_{x}=0,0016 ; D_{y}=0,008 ; L=0,12 ; x(0)=\bar{x}=2 ; y(0)=\bar{y}=2,3\right)$.

Используем теперь метод стрельбы для нахождения стационарных решений задачи 14 , т. е. для решения уравнений
\[
\begin{array}{c}
\frac{1}{\mathrm{Pe}_{\mathrm{M}}} y^{\prime \prime}-y^{\prime}+\mathrm{Da}(1-y) \exp \frac{\Theta}{1+\Theta / \gamma}=0, \\
\frac{1}{\mathrm{Pe}_{\mathrm{H}}} \Theta^{\prime \prime}-\Theta^{\prime}+B \mathrm{Da}(1-y) \exp \frac{\Theta}{1+\Theta / \gamma}-\beta\left(\Theta-\Theta_{c}\right)=0,
\end{array}
\]

соответствующих уравнениям (P14-7), (P14-8). Граничные условия при этом имеют вид (Р14-9), (P14-10), т. е.
\[
\begin{aligned}
\mathrm{Pe}_{\mathrm{M}} y(0)-y^{\prime}(0) & =0, \quad \mathrm{Pe}_{\mathrm{H}} \Theta(0)-\Theta^{\prime}(0)=0, \\
y^{\prime}(1) & =0, \quad \Theta^{\prime}(1)=0 .
\end{aligned}
\]

Выберем два недостающих начальных условия в точке $z=1$ :
\[
y(1)=\eta_{1}, \quad \Theta(1)=\eta_{2} .
\]

Тогда после интегрирования уравнений (6.1.35) с начальными условиями (6.1.36b) и (6.1.37) на промежутке от $z=1$ до $z=0$ мы имеем
\[
F_{1}(\boldsymbol{\eta})=\mathrm{Pe}_{\mathrm{M}} y(0)-y^{\prime}(0)=0, \quad F_{2}(\boldsymbol{\eta})=\mathrm{Pe}_{\mathrm{H}} \Theta(0)-\boldsymbol{\theta}^{\prime}(0)=0 .
\]

Результаты решения уравнений (6.1.38) методом Ньютона с использованием соответствующих вариационных уравнений приведены в табл. 6.4. В ней представлен случай, когда существует пять решений данной задачи. Отметим, что область сходимости метода Ньютона для некоторых решений мала. В частности, так обстоит дело для пятого решения. Сходимость к этому решению даже из близкого к нему начального приближения может быть медленной (см. вторую половину таблицы). Наконец, в табл. 6.5 приведены окончательные решения задачи для нескольких различных значений числа Дамкёлера.

Замечание. Во многих задачах концы интервала равноправны и можно использовать метод стрельбы в любом направлении – как «справа налево», так и «слева направо». В этой задаче есть выделенное направление: решение задачи Коши от $z=0$ к $z=1$ при больших числах Пекле сильно неустойчиво, и такая реализация метода стрельбы здесь не годится.

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

Будем искать решение $y(z), \Theta(z)$ дифференциальных уравнений (Р15-6) и (Р15-7) с граничными условиями вида (6.1.36).

Таблица 6.4. Метод стрельбы для задачи 14, $\mathrm{Pe}_{\mathrm{H}}=5, \mathrm{Pe}_{\mathrm{M}}=10, \gamma=20$, $B=15, \beta=2, \boldsymbol{\theta}_{c}=0, \mathrm{Da}=0,07$.

Таблица 6.5. Метод стрельбы для задачи 14 ( $\mathrm{Pe}_{\mathrm{H}}=5, \mathrm{Pe}_{\mathrm{M}}=10, \gamma=20$, $B=15, \beta=2, \Theta_{c}=0$ ).
Решения для последовательных значений Da

При этом в уравнения входят также неизвестные функции $\omega(z)$ и $\theta(z)$, которые должны удовлетворять соотношениям (Р15-8) и (Р15-9). Для практических вычислений вместо этих нелинейных (алгебраических) соотношений более удобно использовать формулы (P15-10) и (P15-11). При фиксированных значениях параметров. $J_{\mathrm{M}}, J_{\mathrm{H}}$, Da, $B, \gamma$ и заданных значениях функций $y(z)$ и $\Theta(z)$ уравнение (P15-11) представляет собой при каждом $z$ нелинейное уравнение относительно неизвестной $\omega(z)$. Решая уравнение (Р15-11), мы находим значения $\omega(z)$, после чего подставляя их в формулу (Р15-10), вычисляем значения функции $\theta(z)$.

Используем для решения этой задачи метод стрельбы. Два недостающих начальных условия опять выберем в точке $z=1$ в форме (6.1.37). Затем проинтегрируем уравнения (P15-6), (P15-7) с начальными условиями (6.1.36b) и (6.1.37) на промежутке от $z=1$ до $z=0$. На каждом шаге интегрирования нам необходимо вычислять правые части дифференциальных уравнений, в которые наряду с параметрами и значениями функций $y$ и $\Theta$, входят также значения функций $\omega$ и $\theta$. Эти значения определяются из формул (Р15-11) и (Р15-10). После проведенного таким образом интегрирования, когда на каждом шаге решается одно нелинейное уравнение, мы получаем в точке $z=0$ систему двух уравнений (6.1.38) относительно двух неизвестных $\eta_{1}$ и $\eta_{2}$. Эта система может быть решена методом Ньютона с использованием вариационных уравнений для переменных
\[
p_{y i}=\partial y / \partial \eta_{i}, \quad p_{\Theta i}=\partial \Theta / \partial \eta_{i}, \quad i=1,2 .
\]

Первое из этих уравнений имеет вид
\[
p_{y 1}^{\prime \prime}-\mathrm{Pe}_{M} p_{y 1}^{\prime}+\mathrm{Pe}_{M} J_{M}\left(\frac{\partial \omega}{\partial y} p_{y 1}+\frac{\partial \omega}{\partial \Theta} p_{\Theta 1}-p_{y 1}\right)=0 .
\]

Из этого примера видно, что наряду с вариационными переменными, определяемыми формулами (6.1.39), у нас возникают еще и другие переменные
\[
\frac{\partial \omega}{\partial y}, \frac{\partial \omega}{\partial \Theta}, \frac{\partial \theta}{\partial y}, \frac{\partial \theta}{\partial \Theta} .
\]

Значения этих переменных на каждом шаге интегрирования вычисляются из соотношений (Р15-8), (Р15-9), определяющих $\omega=\omega(y, \Theta)$ и $\theta=\theta(y, \Theta)$. Запишем их в виде
\[
G_{i}(y, \Theta, \omega, \theta)=0, \quad i=1,2 .
\]

Предположим, что функции $G_{i}$ удовлетворяют условиям теоремы о неявных функциях, и пусть значения $\omega(y, \Theta)$ и $\theta(y, \Theta)$

для заданных у и $\Theta$ найдены. Тогда производные (6.1.41) можно найти, решая две системы линейных алгебраических уравнений
\[
\begin{array}{c}
\boldsymbol{\Gamma} \cdot\left[\begin{array}{l}
\frac{\partial \omega}{\partial y} \\
\frac{\partial \theta}{\partial y}
\end{array}\right]=-\left[\begin{array}{c}
\frac{\partial G_{1}}{\partial y} \\
\frac{\partial G_{2}}{\partial y}
\end{array}\right], \\
\boldsymbol{\Gamma} \cdot\left[\begin{array}{c}
\frac{\partial \omega}{\partial \Theta} \\
\frac{\partial \theta}{\partial \Theta}
\end{array}\right]=-\left[\begin{array}{c}
\frac{\partial G_{1}}{\partial \Theta} \\
\frac{\partial G_{2}}{\partial \Theta}
\end{array}\right],
\end{array}
\]

где $\partial G_{1} / \partial y=-J_{\mathrm{M}}, \quad \partial G_{2} / \partial y=0, \quad \partial G_{1} / \partial \Theta=0, \quad \partial G_{2} / \partial \Theta=-J_{\mathrm{H}}$, а матрица $\Gamma$ определяется как
\[
\boldsymbol{\Gamma}=\left[\begin{array}{ll}
\frac{\partial G_{1}}{\partial \omega} & \frac{\partial G_{1}}{\partial \theta} \\
\frac{\partial G_{2}}{\partial \omega} & \frac{\partial G_{2}}{\partial \theta}
\end{array}\right] .
\]

Описанный подход несложно реализовать в случае, когда нелинейное уравнение (P15-11) имеет одно и только одно решение $\omega$, которое должно удовлетворять естественному с физической точки зрения требованию $y<\omega<1$ (см. формулу (Р15-5)).

Значительные сложности возникают тогда, когда уравнение (P15-11) при определенной комбинации параметров и переменных $y$ и $\Theta$ будет иметь несколько допустимых решений (как правило, три различных решения $\omega$ ). Такая ситуация имеет место, в частности, при больших значениях параметра $B$. В этом случае метод Ньютона, используемый для решения уравнения (P15-11), может оказаться расходящимся или же будет сходиться к какому-либо другому решению, а не к тому, которое ожидалось. В зависимости от того, какой из корней уравнения (P15-11) выбрать, получаются различные функции $y(z)$ и $\Theta(z)$.

Из физических соображений ${ }^{1)}$ (в случае трех решений) интерес представляют прежде всего два крайних решения $\omega$ уравнения (Р15-11), а именно значение $\omega$, лежащее в окрестности значения $y, \omega>y$ (так называемое нижнее решение), и корень, располагающийся в окрестности $1, \omega<1$ (так называемое верхнее решение). При решении дифференциальных уравнений (P15-6), (Р15-7) мы поступаем следующим образом: в областях параметров, где существует несколько решений уравнения (P15-11), мы всегда рассматриваем (если это возможно) только нижнее решение $\omega$ (или соответственно верхнее решение $\omega$ ) и

тем самым получаем два различных (основных) решения дифференциальных уравнений. Если же комбинировать на разных подынтервалах $z \in[0,1]$ верхнее и нижнее решения $\omega$, то мы получим разрывные профили $\omega(z), \theta(z)$.

На рис. 6.1а приведены профили $\Theta(z), y(z), \omega(z), \theta(z)$, найденные для случая, когда при всех $z$ выбираліось нижнее решение $\omega$ уравнения (P15-11). На рис. 6.1b представлен пример профилей решений (для тех же значений параметров), когда

Рис. 6.1. Аксиальные профили конверсии и температуры, задача $15 ; \mathrm{Da}=$ $=0,065, P e_{M}=20, P e_{H}=10, \beta=1, B=15, \gamma=20, \boldsymbol{\theta}_{c}=0, J_{M}=J_{H}=25$; a) непрерывные профили $\omega$ и $\theta, b$ ) разрывные профили $\omega$ и $\theta$.
$\omega(z)$ и $\theta(z)$ не являются непрерывными функциями (на промежутке $z \in[0,4995 ; 0,5015]$ мы рассматривали верхнее решение $\omega$ ). Из сравнения рисунков видно, что профили решений $y(z)$ и $\Theta(z)$ изменились.

6.1.3. Метод многократной стрельбы

Метод стрельбы, описанный в предыдущем пункте, иногда не позволяет получить удовлетворительные результаты. Так, решение соответствующих задач Коши (включая дифференциальные уравнения в вариациях) может оказаться практически невозможным при наличии сильной чувствительности к начальным условиям. В таких случаях часто оказывается удобным использовать метод многократной стрельбы (см., например, [6.34]). Опишем кратко идею этого метода на примере задачи (6.1.1), (6.1.2) с граничными условиями ГУ2 вида (6.1.4). По аналогии с п. 6.1.1 выберем на промежутке $z \in[0,1]$ сетку узловых точек $z_{0}=0, z_{1}, \ldots, z_{n}=1, z_{i+1}>z_{i}$, которая, вообще говоря,

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

На каждом подынтервале $\left[z_{i}, z_{i+1}\right], i=0, \ldots, n-1$ уравнения (6.1.1), (6.1.2) интегрируются независимо. Для этого нам необходимо задать начальные условия в точках $z_{i}, i=0, \ldots$ $\ldots, n-1$; эти начальные значения мы обозначим как $\boldsymbol{\eta}_{i}$. В точке $z_{0}=0$ уже заданы два условия (6.1.4a), и поэтому здесь, аналогично тому, как это делалось в п. 6.1.2, мы выбираем два дополнительных условия вида (6.1.24). При этом вектор $\eta_{0}$ будет иметь только две составляющих $\eta_{01}=x(0)$ и $\eta_{02}=$ $=y(0)$. Остальные векторы $\eta_{i}, i=1, \ldots, n-1$ будут иметь по четыре составляющих
\[
\boldsymbol{\eta}_{i}=\left(x\left(z_{i}\right), x^{\prime}\left(z_{i}\right), y\left(z_{i}\right), y^{\prime}\left(z_{i}\right)\right) .
\]

Обозначим решение соответствующей задачи Коши на каждом подынтервале $\left[z_{i}, z_{i+1}\right]$ через $x(z)=x\left(z ; z_{i}, \eta_{i}\right)$ (и аналогично для $\left.x^{\prime}(z), y(z), y^{\prime}(z)\right)$. Метод многократной стрельбы состоит в следующем: так подобрать векторы $\boldsymbol{\eta}_{i}, i=0, \ldots$ $\ldots, n-1$, чтобы функции $x(z), x^{\prime}(z), y(z), y^{\prime}(z)$, полученные «частями» на отдельных подынтервалах $\left[z_{i}, z_{i+1}\right]$, оказались непрерывными и чтобы при этом были удовлетворены граничные условия $(6.1 .4 \mathrm{~b})$.

Для этого в точках $z_{i}, i=1, \ldots, n-1$, должны выполняться условия
\[
\begin{array}{c}
\left(x\left(z_{i} ; z_{i-1}, \boldsymbol{\eta}_{i-1}\right), x^{\prime}\left(z_{i} ; z_{i-1}, \boldsymbol{\eta}_{i-1}\right), y\left(z_{i} ; z_{i-1}, \boldsymbol{\eta}_{i-1}\right),\right. \\
\left.y^{\prime}\left(z_{i} ; z_{i-1}, \boldsymbol{\eta}_{i-1}\right)\right)=\boldsymbol{\eta}_{i},
\end{array}
\]

а в точке $z_{n}=1$ (в соответствии с условиями (6.1.4b)) должны быть выполнены соотношения
\[
x^{\prime}\left(1 ; z_{n-1}, \boldsymbol{\eta}_{n-1}\right)=0, \quad y^{\prime}\left(1 ; z_{n-1}, \boldsymbol{\eta}_{n-1}\right)=0 .
\]

Уравнения (6.1.46) и (6.1.47) представляют собой систему из $4(n-1)+2$ нелинейных уравнений относительно $4(n-1)+2$ неизвестных – составляющих векторов $\boldsymbol{\eta}_{0}, \boldsymbol{\eta}_{1}, \ldots, \boldsymbol{\eta}_{n-1}$. Эти уравнения можно решать с помощью любого подходящего метода (например, с помощью метода Ньютона). Для вычисления матрицы Якоби можно опять использовать вариационные дифференциальные уравнения (см. п. 6.1.2). Указанный алгоритм легко модифицируется на случай ГУ1 и ГУЗ; при этом в случае ГУЗ мы выбираем две неизвестных, например $x(0)$, $y(0)$, а затем из уравнений (6.1.5a) находим остальные неизвестные, т. е. $x^{\prime}(0), y^{\prime}(0)$. Далее, с помощью условий (6.1.5b) мы легко получаем соотношения, аналогичные (6.1.47).