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

В полуполосе $0 \leqslant z \leqslant 1,0 \leqslant t<\infty$ выберем равномерную сетку узловых точек $\left(z_{i}, t_{j}\right)$;
\[
\begin{array}{c}
z_{i}=i h, \quad i=0,1, \ldots, n, \quad h=1 / n, \\
t_{j}=j \tau, \quad j=0,1,2, \ldots, \quad \tau>0 .
\end{array}
\]
1) Имеется в виду численное моделирование динамики, т. е. численное решение задачи Коши (точнее, задачи с начальными и краевыми условиями).

Значения приближенного решения в этих узловых точках будем обозначать как
\[
x_{i}^{i} \sim x\left(z_{i}, t_{j}\right), \quad y_{i}^{i} \sim y\left(z_{i}, t_{i}\right) .
\]

Наиболее часто используется шеститочечная разностная схема типа Кранка – Николсона. Для системы (4.3.7) она дает:
\[
\begin{aligned}
\frac{1}{\tau}\left(x_{i}^{j+1}-x_{i}^{j}\right)= & \frac{D_{\mathrm{x}} w}{L^{2} h^{2}}\left(x_{i-1}^{i+1}-2 x_{i}^{j+1}+x_{i+1}^{j+1}\right)+ \\
& +\frac{D_{\mathrm{x}}(1-w)}{L^{2} h^{2}}\left(x_{i-1}^{i}-2 x_{i}^{j}+x_{i+1}^{i}\right)+f_{i}^{j+w}, \\
\frac{1}{\tau}\left(y_{i}^{j+1}-y_{i}^{j}\right)= & \frac{D_{\mathrm{y}} w}{L^{2} h^{2}}\left(y_{i-1}^{j+1}-2 y_{i}^{j+1}+y_{i+1}^{j+1}\right)+ \\
& +\frac{D_{\mathrm{y}}(1-w)}{L^{2} h^{2}}\left(y_{i-1}^{i}-2 y_{i}^{j}+y_{i+1}^{j}\right)+g_{i}^{j+w} .
\end{aligned}
\]

Формулы (6.4.3) включают в себя явную схему ( $=0$ ), чисто неявную схему ( $=1$ ) (обе с погрешностью аппроксимации порядка $O\left(\tau+h^{2}\right)$ ), а также классическую схему Кранка Николсона с погрешностью аппроксимации порядка $O\left(\tau^{2}+\right.$ $\left.+h^{2}\right)(w=1 / 2)^{1)}$. В случае граничных условий 1 рода естественно положить при всех $j$
\[
x_{0}^{j+1}=\bar{x}, \quad x_{n}^{j+1}=\bar{x} ; \quad y_{0}^{j+1}=\bar{y}, \quad y_{n}^{j+1}=\bar{y} .
\]

При этом формулы (6.4.3) записываются для $i=1,2, \ldots, n-1$. В случае граничных условий второго рода (ГУ2) уравнения (6.4.3) рассматриваются для $i=0,1, \ldots, n$, а граничные условия (4.3.12) заменяются соотношениями
\[
x_{-1}^{j+1}=x_{1}^{j+1}, \quad x_{n+1}^{j+1}=x_{n-1}^{j+1} ; \quad y_{-1}^{j+1}=y_{1}^{j+1}, \quad y_{n+1}^{j+1}=y_{n-1}^{j+1} .
\]

Схема (6.4.3) является одношаговой по переменной $t$, а это значит, что новый слой $j+1$ можно вычислить исходя из старого слоя $j$. Это значит, что значения $f^{j+w}$ и $g^{j+w}$ зависят только от $x^{j}, x^{j+1}, y^{j}$ и $y^{j+1}$. Явная схема $(w=0)$ представляет собой наиболее простой случай, поскольку значения неизвестных функций на новом слое (с верхним индексом $j+1$ ) непосредственно выражаются через значения на старом слое (с верхним индексом $j$ ). При этом нелинейные члены $f_{i}^{i}$ и $g_{i}^{j}$ легко аппроксимируются подстановкой значений $x_{i}^{i}$ и $y_{i}^{i}$ в функции $f$ и $g$. Вместе с тем, в случае явной схемы должно выполняться условие численной устойчивости, накладывающее определенные

ограннчения на величину шага $\tau$, а именно
\[
\tau \leqslant \min \left(\frac{L^{2} h^{2}}{2 D_{\mathrm{x}}}, \frac{L^{2} h^{2}}{2 D_{\mathrm{y}}}\right) .
\]

Эффективность указанной явной схемы может оказаться довольно низкой, если из условия (6.4.6) приходится выбирать очень малый шаг $\tau$.

При $w
eq 0$ соотношение (6.4.3) представляет собой уже систему уравнений относительно неизвестных с верхним индексом $j+1$. Такой метод называется неявным. При $w \geqslant 1 / 2$ метод оказывается численно устойчивым, причем длина шага $\tau$ ничем не ограничивается. Записав нелинейные члены $f_{i}^{j+w}, g_{i}^{j+w}$ в уравнениях (6.4.3) в виде
\[
f_{i}^{j+w} \sim f\left(w x_{i}^{j+1}+(1-w) x_{i}^{j}, \quad w y_{i}^{j+1}+(1-w) y_{i}^{j}\right)
\]
(член $g_{i}^{j+w}$ записывается аналогично), мы получаем систему $2 n+2$ нелинейных уравнений, которую можно решать итерациями (например, методом Ньютона). Поскольку у нас имеется хорошее начальное приближение (а именно значения на старом слое $j$ ), обыкновенно оказывается достаточным одной-двух итераций. Метод Ньютона предпочтительнее использовать при анализе быстрых химических реакций (т. е. тогда, когда производные $f$ и $g$ велики. – Peд.). При этом число итераций может оказаться более высоким. Для вычисления нового слоя используются также неитерационные подходы, в которых нелинейные члены аппроксимируются иначе – исходя из уже найденных значений неизвестных. По существу в таких подходах мы имеем дело с одной итерацией, проводимой каким-либо методом, аналогичным методу Ньютона.
Остановимся на этом вопросе несколько подробнее.

6.4.1. Замена нелинейных членов

Простейшая разностная замена нелинейных членов в схеме (6.4.3) заключается в использовании линеаризации по $x^{j+1}$ и $y^{j+1}$, т. е. для члена $f_{i}^{i+w}$ мы полагаем
\[
f_{i}^{i+w} \sim a i+b_{i}^{i} x_{i}^{i+1}+c_{i}^{i} y_{i}^{i+1},
\]

тде коэффициенты $a_{i}^{l}, b_{i}^{i}, c_{i}^{i}$ зависят от известных значений решения на старом ( $j$-м) слое. Так, например, при $w=1$ член

$\left(x_{i}^{i+1}\right)^{2}$ можно заменить ${ }^{1)}$ произведением $\left(x_{i}^{i+1} \cdot x_{i}^{j}\right)$; точно так же $\left(x_{i}^{i+1}\right)^{3} \sim\left(x_{i}^{i}\right)^{2} x_{i}^{j+1}$. При $=1 / 2$ член $\left(x_{i}^{i+1 / 2}\right)^{2}$ можно заменить выражением $\left[\left(x_{i}^{j}+x_{i}^{i+1}\right) / 2\right]^{2} \sim \frac{1}{4}\left[\left(x_{i}^{i}\right)^{2}+2 x_{i}^{i} x_{i}^{j+1}+x_{i}^{i} x_{i}^{j+1}\right]$. Довольно часто используется еще более простая замена
\[
f_{i}^{i+w} \sim f_{i}^{!} .
\]

Если функции $f$ и $g$ легко продифференцировать аналитически, то можно воспользоваться формулой Тейлора, т. е. (например, для $f_{i}^{j+w}$ ) положить ${ }^{2)}$ (ср. с (6.4.7))
\[
f_{i}^{i+w} \sim f_{i}^{j}+\left.\frac{\partial f}{\partial x}\right|_{i} ^{i}\left(x_{i}^{j+w}-x_{i}^{i}\right)+\left.\frac{\partial f}{\partial y}\right|_{i} ^{i}\left(y_{i}^{i+w}-y_{i}^{i}\right){ }^{3)}
\]

Здесь $x_{i}^{i+w} \sim w x_{i}^{j+1}+(1-w) x_{i}^{i}$ (аналогично определяется $\left.y_{i}^{i+w}\right)$.
Если аналитическое дифференцирование $f$ и $g$ оказывается сложным, используются процедуры экстраполяции, при которых значение $x_{i}^{j+w}$ для вычисления $f_{i}^{j+w}$ получают экстраполяцией, исходя из значений на нескольких предыдущих слоях (с верхними индексами $j, j-1, \ldots)$. Например, полагают
\[
\begin{array}{c}
x_{i}^{i+w} \sim(1+w) x_{i}^{i}-w x_{i}^{j-1}, \\
x_{i}^{j+w} \sim\left[1+1,5 w+0,5 w^{2}\right] x_{i}^{j}-\left(w^{2}+2 w\right) x_{i}^{j-1}+0,5\left(w^{2}+w\right) x_{i}^{j-2} .
\end{array}
\]

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

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

1) Здесь имеется в виду, что $f(x, y)$ (или $g(x, y)$ ) содержит слагаемое $x^{2}$, и обсуждается, какое слагаемое можно ввести в определение величины. $f^{i+w}$.
2) Разумно сначала определить $f^{f+w}$ по (6.4.7) и потом оставить только первые члены тейлоровского разложения. По существу об этом уже шла речь (см. абзац после формулы (6.4.7)).
${ }^{3}$ ) Здесь имеется в виду, что (по определению) $f_{i}^{j+w}=f\left(x_{i}^{j+w}, y_{i}^{j+w}\right)$

6.4.2. Автоматическое изменение шага по времени

Нахождение значений на одном слое при больших $n$ может потребовать много вычислений, в связи с чем часто оказывается выгодным регулировать длину шага по времени $\tau$ в зависимости от требуемой точности решения. Характер самого решения и его вариаций во времени может изменяться, а вместе с ними будет существенно изменяться и шаг $\tau$, который обеспечивает заданную точность. В связи с этим включение в алгоритм изменения шага по времени может значительно повысить эффективность. используемой схемы. Основная проблема при этом – оценка погрешности решения, вызванной дискретизацией по переменной $t$. Отметим, что используемые подходы носят эвристический характер. Опишем здесь три подхода.
a) Сравнение значений решения на новом слое, полученных с шагом $\tau\left(x^{i+1(\tau)}\right)$, со значениями, полученными с помощью двух шагов длиной $\tau / 2\left(x^{j+1(\tau / 2)}\right)$. В этом случае оценкой погрешности с целью регулирования шага может служить величина ${ }^{1)}$
\[
E_{\mathrm{x}}=\left\|x^{j+1(\tau)}-x^{j+1(\tau / 2)}\right\|
\]
(аналогичная формула записывается для $E_{\text {y }}$ ). Для более тонкой оценки $E$ можно также использовать экстраполяцию по Ричардсону [6.10]. Если действовать так, то для одного удачного и контролируемого шага $\tau$ нам потребуется сделать три «пробных» шага. По этой причине указанный подход используется сравнительно редко и в большинстве случаев после нескольких постоянных (т. е. не регулируемых) шагов по времени.
b) Оценка различия между экстраполированным значением $x^{j+w}$, найденным с помощью формул (6.4.11) или (6.4.12), и вычисленным значением $x^{j+w}$ :
\[
E_{\mathrm{x}}=\max _{i}\left|x_{i}^{j+w}-w x_{i}^{j+1}-(1-w) x_{i}^{j}\right| .
\]

Аналогично определяется $E_{\mathrm{y}}$. Такую же оценку можно получить и для метода «предиктор-корректор».
c) В случае использования какой-то схемы линеаризации (6.4.8) мы можем положить
\[
\begin{aligned}
E_{\mathrm{x}}=\max _{i} \mid f\left(w x_{i}^{i+1}+(1-w) x_{i}^{j}, w_{y_{i}^{i}}+1\right. & \left.+(1-w) y_{i}^{j}\right)- \\
& -a_{i}^{j}-b_{i}^{i} x_{i}^{j+1}-c_{i}^{j} y_{i}^{i+1} \mid .
\end{aligned}
\]

Аналогичным образом записывается и оценка $E_{\text {у. }}$. Отметим, что в данном случае мы находим отклонение в значениях функций $f$ и $g$, но не в переменных состояния $x$ и $y$.

1) Обычно норму \|\| определяют так: $\|u\|=\max \left|u_{i}\right|$.

Каждую из величин $E=\max \left(E_{\mathrm{x}}, E_{\mathrm{y}}\right)$, приведенных выше, можно использовать для регулировки шага по времени $\tau$. С этой целью обычно выбирается некоторое характерное значение $\varepsilon$, причем шаг $\tau$ уменьшается, если $E>\varepsilon$. Если же $E$ существенно меньше, чем $\varepsilon$, то шаг $\tau$ увеличивается.

Мы ограничились здесь лишь идейной стороной алгоритма изменения шага $\tau$. Очень вероятно, что действительная величина погрешности решения будет заметно отличаться от $\varepsilon$. Поэтому если мы хотим удостовериться в том, что заданная точность обеспечивается, следует провести вычисления при двух различных значениях $\varepsilon$ (отличающихся, например, на порядок) и сравнить полученные результаты. Таким образом, однако, можно оценить лишь влияние погрешностей аппроксимации «в направлении оси $t$ ». Влияние погрешности аппроксимации «в направлении координаты $z$ » и адаптационные алгоритмы для выбора шага $h$ будут обсуждаться в следующем пункте.

6.4.3. Автоматическое изменение шага $h$

В этом пункте мы исследуем проблему автоматического изменения шага $h$ на равномерной сетке. Для решения вопроса о том, достаточно ли велико $n$ (т. е. достаточно ли мал шаг $h$ ), при заданной «характерной величине погрешности» $\varepsilon$, будем искать оценку погрешности аппроксимации в направлении оси $z$. В принципе мы могли бы провести вычисления для момента $t_{j}$ дважды – ш шагом $h$ и $2 h$ – и результаты сравнить. Такой численный подход не был бы, однако, достаточно простым и эффективным.

Другой, более подходящий путь заключается в том, чтобы использовать выражения для остаточного члена в соответствующих разностных формулах. При замене второй производной в (6.4.3) имеем
\[
\frac{\partial^{2} x}{\partial z^{2}}\left(z_{i}, t_{j}\right)=\frac{1}{h^{2}}\left[x_{i-1}^{i}-2 x_{i}^{j}+x_{i+1}^{l}\right]-\frac{1}{12} h^{2} \frac{\partial^{4} x}{\partial z^{4}}\left(\xi, t_{j}\right),
\]

тде $\xi \in\left(z_{i-1}, z_{t+1}\right)$. Для аппроксимации четвертой производной можно воспользоваться разностной формулой вида
\[
\frac{\partial^{4} x}{\partial z^{4}}\left(\xi, t_{j}\right) \sim \frac{1}{h^{4}}\left[x_{i-2}^{j}-4 x_{-1}^{i}+6 x_{i}^{j}-4 x_{i+1}^{j}+x_{i+2}^{j}\right]
\]

при $i=2,3, \ldots, n-2$. Тогда оценку погрешности «в направлении оси $z$ » можно записать в виде
\[
E_{\mathrm{x}}=\max _{i=2 \ldots \ldots n-2} \frac{1}{12 h^{2}}\left|x_{i-2}^{j}-4 x_{i-1}^{j}+6 x_{i}^{j}-4 x_{i+1}^{j}+x_{i+2}^{j}\right| .
\]

Аналогичное выражение получается и для $E_{\text {у }}$ (при $i=1$ и $i=$ $=n-1$ для четвертой производной можно использовать асимметричные формулы). Обозначим $E=\max \left(E_{\mathrm{x}}, E_{\mathrm{y}}\right)$.

Тогда для изменения (регулировки) $n$ можно использовать следующую стратегию:
1. $E>\varepsilon$; шаг $h$ слишком велик, и его следует уменьшить вдвое. Значения решения в узловых точках, которые появятся между исходными узловыми точками, после уменьшения шага, можно определить с помощью интерполяции, например, положить
\[
\begin{aligned}
x_{i+1 / 2}^{j} & =\frac{1}{16}\left[-x_{i-1}^{j}+9 x_{i}^{j}+9 x_{i+1}^{i}-x_{i+2}^{j}\right], i=1,2, \ldots, n-2, \\
x_{1 / 2}^{j} & =\frac{1}{8}\left[3 x_{0}^{j}+6 x_{1}^{j}-x_{2}^{j}\right], \\
x_{n-1 / 2}^{j} & =\frac{1 *}{8}\left[3 x_{n}^{j}+6 x_{n-1}^{j}-x_{n-2}^{j}\right] .
\end{aligned}
\]

После проведения интерполяции перенумеруем заново узловые точки, удвоим число $n$ и вдвое уменьшим шаг $h$; затем продолжим вычисления.
2. $E<\varepsilon / \omega$, где значение $\omega$, как правило, выбирается между 4 и 8. Шаг $h$ излишне мал и может быть удвоен. При этом узловые точки с нечетными индексами $1,3,5, \ldots, n-1$ выбрасываются, оставшиеся точки перенумеровываются, шаг $h$ удваивается, а $n$ уменьшается вдвое.
3. В остальных случаях шаг $h$ остается неизменным и проводятся дальнейшие вычисления.

При формулировке алгоритма удобно задавать максимальный и минимальный допустимый шаг $h$.

6.4.4. Адаптивная неравномерная сетка

Пусть у нас имеется некоторая, вообще говоря, неравномерная сетка узловых точек $z_{0}=0, z_{1}, \ldots, z_{n}=1$. Тогда вторую производную в направлении $z$ в точке $\left(z_{i}, t_{j}\right)$ можно заменить трехточечной разностной формулой
\[
\begin{array}{c}
\frac{\partial^{2} x}{\partial z^{2}}\left(z_{i}, t_{j}\right) \sim \frac{2}{\left(z_{i}-z_{i-1}\right)\left(z_{i+1}-z_{i-1}\right)} x_{i-1}^{i}- \\
-\frac{2}{\left(z_{i}-z_{i-1}\right)\left(z_{i+1}-z_{i}\right)} x_{i}^{i}+\frac{2}{\left(z_{i+1}-z_{i-1}\right)\left(z_{i+1}-z_{i}\right)} x_{i+1}^{i} .
\end{array}
\]

При этом разностные уравнения (6.4.3) изменяются только за счет подстановки формулы (6.4.20) вместо простейшей фор-

мулы ${ }^{1)}$ для $\frac{\partial^{2} x}{\partial z^{2}}$. Все остальное остается неизменным, включая вычисления нового слоя $j+1$. Погрешность аппроксимации в направлении оси $z$ оценивается на каждом подынтервале $\left(z_{i-1}, z_{i}\right)$ отдельно. Айгенбергер и Батт [6.29] предложили следующий подход. Приближенное значение решения при $t=t_{j}$ и $z=\left(z_{i-1}+z_{i}\right) / 2$ отыскивается двумя способами (ниже индекс $j$ опущен). В первом случае мы получаем его интерполированием по узловым точкам $z_{i-2}, z_{i-1}, z_{i}$,
\[
\begin{array}{c}
x_{i-1 / 2}^{(1)}=\frac{-\left(z_{i}-z_{i-1}\right)^{2}}{4\left(z_{i}-z_{i-2}\right)\left(z_{i-1}-z_{i-2}\right)} x_{i-2}+\frac{z_{i}+z_{i-1}-2 z_{i-2}}{4\left(z_{i-1}-z_{i-2}\right)} x_{i-1}+ \\
+\frac{z_{i}+z_{i-1}-2 z_{i-2}}{4\left(z_{i}-z_{i-2}\right)} x_{i},
\end{array}
\]

во втором случае – интерполированием по узловым точкам $z_{i-1}, z_{i}, z_{i+1}$,
\[
\begin{array}{c}
x_{i-1 / 2}^{(2)}=\frac{2 z_{i+1}-z_{i}-z_{i-1}}{4\left(z_{i+1}-z_{i-1}\right)} x_{i-1}+\frac{2 z_{i+1}-z_{i}-z_{i-1}}{4\left(z_{i+1}-z_{i}\right)} x_{i}- \\
-\frac{\left(z_{i}-z_{i-1}\right)^{2}}{4\left(z_{i+1}-z_{i-1}\right)\left(z_{i+1}-z_{i}\right)} x_{i+1} .
\end{array}
\]

Для оценки погрешности аппроксимации на интервале $\left(z_{i-1}, z_{i}\right.$ ) используется величина
\[
E_{\mathrm{x}}^{i}=\left|x_{i-1 / 2}^{(1)}-x_{i-1 / 2}^{(2)}\right| ;
\]
$E_{\mathrm{y}}^{i}$ определяется аналогично. Тогда $E^{i}=\max \left(E_{x}^{i}, E_{y}^{i}\right)$. Далее сетка узловых точек видоизменяется с помощью следующего алгоритма (здесь снова $\varepsilon$ – заданная характерная величина погрешности).

1. $E^{i}>\varepsilon$. Тогда между $z_{i}$ и $z_{i-1}$ мы добавляем еще одну узловую точку $\left(z_{i-1}+z_{i}\right) / 2$. Приближенное значение решения в этой точке находим посредством интерполяции по четырем соседним узловым точкам $z_{i-2}, z_{i-1}, z_{i}, z_{i+1}$. Проделаем эту операцию для $i=1,2, \ldots, n$. Далее перенумеруем узловые точки (добавим новые точки с сохранением порядка), после чего продолжим вычисления на построенной более густой сетке.
2. $E^{i}<\varepsilon / \omega$ для $i=k$ и одновременно для $i=k+1$. В этом случае мы отбрасываем узловую точку $z_{k}$ и вновь проводим
1) $\frac{\partial^{2} x}{\partial z^{2}} \sim \frac{1}{h^{2}}\left(x_{i-1}-2 x_{i}+x_{i+1}\right)$. Формула (6.4.20), очевидно, переходит в простейшую при $z_{i}-z_{i-1}=z_{i+1}-z_{i}=h$.

соответствующие вычисления $E^{i}$ для соседних узловых точек. Затем перенумеровываем узловые точки и продолжаем вычисления. При этом значение $\omega$ выбирается обычно в диапазоне $10-20$.
3. В остальных случаях сетка не изменяется. При этом бывает удобным задавать максимальное и минимальное расстояние между двумя соседними точками.

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

6.4.5. Метод прямых

С методом прямых мы уже встречались в п. 6.3 .3 (см. уравнения (6.3.31)). Его применение оказывается эффективным в тех случаях, когда число узловых точек сравнительно невелико. Последнего можно добиться увеличением порядка аппроксимации. В данном пункте на примере системы типа «реакция-диффузия» (4.3.7) будет продемонстрировано использование этого метода для случая, когда пространственные производные заменяются разностными отношениями с пятью узловыми точками. При этом порядок аппроксимации по сравнению с трехточечной схемой (6.3.31) возрастает с $O\left(h^{2}\right)$ до $O\left(h^{4}\right)$. Обозначим $x_{i}(t) \approx$ $\approx x\left(z_{i}, t\right), y_{i}(t) \approx y\left(z_{i}, t\right)$ и рассмотрим равномерную сетку $z_{i}=$ $=i h, i=0,1, \ldots, n ; h=1 / n$. Уравнение (4.3.7a) заменяется системой уравнений
\[
\begin{array}{l}
\frac{d x_{i}}{d t}=\frac{D_{\mathrm{x}}}{12 L^{2} h^{2}}\left(-x_{i-2}+16 x_{i-1}-30 x_{i}+16 x_{i+1}-x_{i+2}\right)+f\left(x_{i}, y_{i}\right) \text {, } \\
i=2,3, \ldots, n-2 \text {, } \\
\frac{d x_{1}}{d t}=\frac{D_{\mathrm{x}}}{12 L^{2} h^{2}}\left(11 x_{0}-20 x_{1}+6 x_{2}+4 x_{3}-x_{4}\right)+f\left(x_{1}, y_{1}\right) \text {, } \\
\frac{d x_{n-1}}{d t}=\frac{D_{\mathrm{x}}}{12 L^{2} h^{2}}\left(11 x_{n}-20 x_{n-1}+6 x_{n-2}+\right. \\
\left.+4 x_{n-3}-x_{n-4}\right)+f\left(x_{n-1}, y_{n-1}\right) \text {. } \\
\end{array}
\]

В случае граничных условий 1-го рода вида (4.3.9), мы имеем $x_{0}(t) \equiv \bar{x}, x_{n}(t) \equiv \bar{x}$. Для уравнения (4.3.7b) соответствующая замена строится совершенно аналогично. Таким образом, мы получаем систему обыкновенных дифференциальных уравнений, которую можно решить каким-либо из методов, описанных в $\$ 5.7$, например методом Рунге-Кутты.

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

ограничение, которое определяется требованием численной устойчивости. Если же мы применяем схему интегрирования с автоматическим изменением шага (например, схему Мерсона), то будет выбран относительно короткий шаг интегрирования ( $\sim h^{2}$ ), что, как правило, приводит к большим затратам машинного времени.

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

Продемонстрируем указанный подход, построив аппроксимацию решения задачи 16 с помощью метода прямых [6.10]. Зададим граничные условия (P16-13) при $\mathrm{Nu} \rightarrow \infty, \mathrm{Sh} \rightarrow \infty$ и выберем $n$ узловых точек $r_{1}, r_{2}, \ldots, r_{n}=1$. Предположим теперь, что пространственный дифференциальный оператор в уравнении (Р16-11) заменяется в узловой точке $r_{j}$ некоторой линейной комбинацией значений решения в узловых точках, а именно
\[
\frac{\partial^{2} y}{\partial r^{2}}+\left.\frac{a}{r} \frac{\partial y}{\partial r}\right|_{r=r_{j}} \sim \sum_{i=1}^{n} A_{i j} y\left(r_{i}, t\right) .
\]

Коэффициенты $A_{i j}$ находятся из требования, чтобы левая часть соотношения (6.4.25) точно равнялась правой для функций $y=1, \quad y=r^{2}, \quad y=r^{4}, \ldots, \quad y=r^{2 n-2}$. Последовательно подставляя эти функции в формулу (6.4.25), получим
\[
\begin{array}{c}
0=\sum_{i=1}^{n} A_{i j}, \\
2+\frac{a}{r_{j}} 2 r_{j}=\sum_{i=1}^{n} A_{i j} r_{i}^{2} \\
(2 n-2)(2 n-3) r_{j}^{2 n-4}+\frac{a}{r_{j}}(2 n-2) r_{j}^{2 n-3}=\sum_{i=1}^{n} A_{i j} r_{i}^{2 n-2} .
\end{array}
\]

При каждом фиксированном $j$ соотношения (6.4.26) представляют собой систему $n$ линейных алгебраических уравнений
относительно неизвестных коэффициентов $A_{i j}$. Найдем решения этих систем при $j=1,2, \ldots, n-1$ (в случае $j=n$ это оказывается излишним, поскольку в точке $r=1$ у нас задано $y=1$ ).

Используя аппроксимацию (6.4.25) для $у$ и $\Theta$, мы получаем из исходных уравнений (P16-11), (P16-12) систему обыкновенных дифференциальных уравнений для функций $y_{j}(t)=y\left(r_{j}, t\right)$, $\Theta_{j}=\Theta\left(r_{j}, t\right), j=1,2 \ldots, n-1$ :
\[
\begin{array}{c}
\mathrm{Lw} \frac{d y_{j}}{d t}=\sum_{i=1}^{n} A_{i j} y_{i}-\Phi^{2} y_{j} \exp \frac{\Theta_{j}}{1+\Theta_{j} / \gamma}, \\
\frac{d \Theta_{j}}{d t}=\sum_{i=1}^{n} A_{i j} \Theta_{i}+\gamma \beta \Phi^{2} y_{j} \exp \frac{\Theta_{j}}{1+\Theta_{j} / \gamma}
\end{array}
\]
(напомним, что $y_{n}(t) \equiv 1$ и $\Theta_{n}(t) \equiv 0$ ).
Положим теперь $n=2$, т. е. выберем лишь одну внутреннюю точку $r_{1}$, а также точку $r_{2}=1$. Система уравнений (6.4.26) для $j=1$ приобретает вид
\[
A_{21}+A_{11}=0, \quad A_{21}+r_{1}^{2} A_{11}=2(1+a)
\]

откуда
\[
A_{11}=\frac{2(1+a)}{r_{1}^{2}-1}=\rho, \quad A_{21}=-\frac{2(1+a)}{r_{1}^{2}-1}=-\rho .
\]

Дифференциальные уравнения (6.4.27), (6.4.28) принимают в этом случае вид
\[
\begin{aligned}
\text { Lw } \frac{d y_{1}}{d t} & =\rho\left(y_{1}-1\right)-\Phi^{2} y_{1} \exp \frac{\Theta_{1}}{1+\Theta_{1} / \gamma}, \\
\frac{d \Theta_{1}}{d t} & =\rho \Theta_{1}+\gamma \beta \Phi^{2} y_{1} \exp \frac{\Theta_{1}}{1+\Theta_{1} / \gamma}
\end{aligned}
\]
(здесь уже использованы условия $y_{2}(t) \equiv 1, \Theta_{2}(t) \equiv 0$ ).
Весьма существенным является в данном случае выбор координаты $r_{1}$. Вилладсен и Стюарт [6.30] рекомендуют выбирать в качестве $r_{1}$ нуль подходящего ортогонального полинома. Рекомендуемые ими значения $r_{1}$ (нули неких полиномов Якоби) даются следующей таблицей:

Для стационарного решения системы (6.4.30) имеет место соотношение
\[
\Theta_{1}=\gamma \beta\left(1-y_{1}\right),
\]

аналогичное формуле (P16-15). Подставляя его во второе из уравнений (6.4.30) (при $\Theta_{1}=0$ ), получаем
\[
\rho\left(y_{1}-1\right)-\Phi^{2} y_{1} \exp \frac{\gamma \beta\left(1-y_{1}\right)}{1+\beta\left(1-y_{1}\right)}=0 .
\]

В табл. 6.12 приведена зависимость $y_{1}$ от Ф, подсчитанная с помощью метода отображения параметра Ф для последовательности значений $y_{1}$. Сравнивая ее с табл. 6.9, можно заметить качественное совпадение зависимостей решений от параметра ( $y_{1}$ соответствует значению $r=0,4472 !$ ). При $\beta=0,05$ зависимость однозначна, при $\beta \approx 0,25$ у нас возникают кратные решения, а при $\beta=0,4$ в некотором диапазоне изменения параметра Ф существуют три стационарных решения данной задачи.

Таблища 6.12. Зависимость решений уравнения (6.4.32) от параметра Ф для задачи 16, $\gamma=20, a=0, \rho=-2,5$.

Динамическое поведение системы (6.4.30) может подсказать нам, как будет вести себя исходная система двух дифференциальных уравнений с частными производными, динамическое моделирование которой требует больших затрат машинного времени. Выбор $n=3$ может дать нам более точные результаты. Более подробно результаты исследования задачи 16 приведены в работе [6.19], а для других задач – в работе [6.4].