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

В полуполосе $0⩽z⩽1,0⩽t<\mathrm{\infty }$ выберем равномерную сетку узловых точек $(z_i,t_j)$;
$$\begin{array}{c}{z}_i=ih,i=0,1,\dots ,n,h=1/n,\\ t_j=j\tau ,j=0,1,2,\dots ,\tau >0.\end{array}$$
1) Имеется в виду численное моделирование динамики, т. е. численное решение задачи Коши (точнее, задачи с начальными и краевыми условиями).

Значения приближенного решения в этих узловых точках будем обозначать как
$$x_i^i\sim x(z_i,t_j),y_i^i\sim y(z_i,t_i).$$

Наиболее часто используется шеститочечная разностная схема типа Кранка — Николсона. Для системы (4.3.7) она дает:
$$\begin{array}{rl}\frac{1}{\tau }(x_i^{j+1}-x_i^j)=& \frac{D_{\mathrm{x}}w}{L^2{h}^2}(x_{i-1}^{i+1}-2x_i^{j+1}+x_{i+1}^{j+1})+\\ & +\frac{D_{\mathrm{x}}(1-w)}{L^2{h}^2}(x_{i-1}^i-2x_i^j+x_{i+1}^i)+f_i^{j+w},\\ \frac{1}{\tau }(y_i^{j+1}-y_i^j)=& \frac{D_{\mathrm{y}}w}{L^2{h}^2}(y_{i-1}^{j+1}-2y_i^{j+1}+y_{i+1}^{j+1})+\\ & +\frac{D_{\mathrm{y}}(1-w)}{L^2{h}^2}(y_{i-1}^i-2y_i^j+y_{i+1}^j)+g_i^{j+w}.\end{array}$$

Формулы (6.4.3) включают в себя явную схему ( $=0$ ), чисто неявную схему ( $=1$ ) (обе с погрешностью аппроксимации порядка $O(\tau +h^2)$ ), а также классическую схему Кранка Николсона с погрешностью аппроксимации порядка $O({\tau }^2+$ $+h^2)(w=1/2{)}^{1)}$. В случае граничных условий 1 рода естественно положить при всех $j$
$$x_0^{j+1}=\overline{x},x_n^{j+1}=\overline{x};y_0^{j+1}=\overline{y},y_n^{j+1}=\overline{y}.$$

При этом формулы (6.4.3) записываются для $i=1,2,\dots ,n-1$. В случае граничных условий второго рода (ГУ2) уравнения (6.4.3) рассматриваются для $i=0,1,\dots ,n$, а граничные условия (4.3.12) заменяются соотношениями
$$x_{-1}^{j+1}=x_1^{j+1},x_{n+1}^{j+1}=x_{n-1}^{j+1};y_{-1}^{j+1}=y_1^{j+1},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 ⩽min(\frac{L^2{h}^2}{2D_{\mathrm{x}}},\frac{L^2{h}^2}{2D_{\mathrm{y}}}).$$

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

При $weq0$ соотношение (6.4.3) представляет собой уже систему уравнений относительно неизвестных с верхним индексом $j+1$. Такой метод называется неявным. При $w⩾1/2$ метод оказывается численно устойчивым, причем длина шага $\tau$ ничем не ограничивается. Записав нелинейные члены $f_i^{j+w},g_i^{j+w}$ в уравнениях (6.4.3) в виде
$$f_i^{j+w}\sim f(w{x}_i^{j+1}+(1-w)x_i^j,w{y}_i^{j+1}+(1-w)y_i^j)$$
(член $g_i^{j+w}$ записывается аналогично), мы получаем систему $2n+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 ai+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)}$ произведением $(x_i^{i+1}\cdot x_i^j)$; точно так же ${\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[(x_i^j+x_i^{i+1})/2\right]}^2\sim \frac{1}{4}[{\left(x_i^i\right)}^2+2x_i^i{x}_i^{j+1}+x_i^i{x}_i^{j+1}]$. Довольно часто используется еще более простая замена
$$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+{\frac{\partial f}{\partial x}|}_i^i(x_i^{j+w}-x_i^i)+{\frac{\partial f}{\partial y}|}_i^i(y_i^{i+w}-y_i^i){}^{3)}$$

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

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

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

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(x_i^{j+w},y_i^{j+w})$

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}}=\Vert x^{j+1(\tau )}-x^{j+1(\tau /2)}\Vert$$
(аналогичная формула записывается для $E_{\text{y }}$ ). Для более тонкой оценки $E$ можно также использовать экстраполяцию по Ричардсону [6.10]. Если действовать так, то для одного удачного и контролируемого шага $\tau$ нам потребуется сделать три «пробных» шага. По этой причине указанный подход используется сравнительно редко и в большинстве случаев после нескольких постоянных (т. е. не регулируемых) шагов по времени.
b) Оценка различия между экстраполированным значением $x^{j+w}$, найденным с помощью формул (6.4.11) или (6.4.12), и вычисленным значением $x^{j+w}$ :
$$E_{\mathrm{x}}=\underset{i}{max}|x_i^{j+w}-w{x}_i^{j+1}-(1-w)x_i^j|.$$

Аналогично определяется $E_{\mathrm{y}}$. Такую же оценку можно получить и для метода «предиктор-корректор».
c) В случае использования какой-то схемы линеаризации (6.4.8) мы можем положить
$$\begin{array}{rl}{E}_{\mathrm{x}}=\underset{i}{max}\mid f(w{x}_i^{i+1}+(1-w)x_i^j,w_{y_i^i}+1& +(1-w)y_i^j)-\\ & -a_i^j-b_i^i{x}_i^{j+1}-c_i^j{y}_i^{i+1}\mid .\end{array}$$

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

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

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

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

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

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

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

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

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

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

Тогда для изменения (регулировки) $n$ можно использовать следующую стратегию:
1. $E>\epsilon$; шаг $h$ слишком велик, и его следует уменьшить вдвое. Значения решения в узловых точках, которые появятся между исходными узловыми точками, после уменьшения шага, можно определить с помощью интерполяции, например, положить
$$\begin{array}{rl}{x}_{i+1/2}^j& =\frac{1}{16}[-x_{i-1}^j+9x_i^j+9x_{i+1}^i-x_{i+2}^j],i=1,2,\dots ,n-2,\\ x_{1/2}^j& =\frac{1}{8}[3x_0^j+6x_1^j-x_2^j],\\ x_{n-1/2}^j& =\frac{1\ast }{8}[3x_n^j+6x_{n-1}^j-x_{n-2}^j].\end{array}$$

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

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

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

Пусть у нас имеется некоторая, вообще говоря, неравномерная сетка узловых точек $z_0=0,z_1,\dots ,z_n=1$. Тогда вторую производную в направлении $z$ в точке $(z_i,t_j)$ можно заменить трехточечной разностной формулой
$$\begin{array}{c}\frac{{\partial }^{2}x}{\partial z^2}(z_i,t_j)\sim \frac{2}{(z_i-z_{i-1})(z_{i+1}-z_{i-1})}{x}_{i-1}^i-\\ -\frac{2}{(z_i-z_{i-1})(z_{i+1}-z_i)}{x}_i^i+\frac{2}{(z_{i+1}-z_{i-1})(z_{i+1}-z_i)}{x}_{i+1}^i.\end{array}$$

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

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

во втором случае — интерполированием по узловым точкам $z_{i-1},z_i,z_{i+1}$,
$$\begin{array}{c}{x}_{i-1/2}^{(2)}=\frac{2z_{i+1}-z_i-z_{i-1}}{4(z_{i+1}-z_{i-1})}{x}_{i-1}+\frac{2z_{i+1}-z_i-z_{i-1}}{4(z_{i+1}-z_i)}{x}_i-\\ -\frac{{(z_i-z_{i-1})}^2}{4(z_{i+1}-z_{i-1})(z_{i+1}-z_i)}{x}_{i+1}.\end{array}$$

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

1. $E^i>\epsilon$. Тогда между $z_i$ и $z_{i-1}$ мы добавляем еще одну узловую точку $(z_{i-1}+z_i)/2$. Приближенное значение решения в этой точке находим посредством интерполяции по четырем соседним узловым точкам $z_{i-2},z_{i-1},z_i,z_{i+1}$. Проделаем эту операцию для $i=1,2,\dots ,n$. Далее перенумеруем узловые точки (добавим новые точки с сохранением порядка), после чего продолжим вычисления на построенной более густой сетке.
2. $E^i<\epsilon /\omega$ для $i=k$ и одновременно для $i=k+1$. В этом случае мы отбрасываем узловую точку $z_k$ и вновь проводим
1) $\frac{{\partial }^{2}x}{\partial z^2}\sim \frac{1}{h^2}(x_{i-1}-2x_i+x_{i+1})$. Формула (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(z_i,t),y_i(t)\approx y(z_i,t)$ и рассмотрим равномерную сетку $z_i=$ $=ih,i=0,1,\dots ,n;h=1/n$. Уравнение (4.3.7a) заменяется системой уравнений
$$\begin{array}{l}\frac{d{x}_i}{dt}=\frac{D_{\mathrm{x}}}{12L^2{h}^2}(-x_{i-2}+16x_{i-1}-30x_i+16x_{i+1}-x_{i+2})+f(x_i,y_i)\text{, }\\ i=2,3,\dots ,n-2\text{, }\\ \frac{d{x}_1}{dt}=\frac{D_{\mathrm{x}}}{12L^2{h}^2}(11x_0-20x_1+6x_2+4x_3-x_4)+f(x_1,y_1)\text{, }\\ \frac{d{x}_{n-1}}{dt}=\frac{D_{\mathrm{x}}}{12L^2{h}^2}(11x_n-20x_{n-1}+6x_{n-2}+\\ +4x_{n-3}-x_{n-4})+f(x_{n-1},y_{n-1})\text{. }\end{array}$$

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

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

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

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

Продемонстрируем указанный подход, построив аппроксимацию решения задачи 16 с помощью метода прямых [6.10]. Зададим граничные условия (P16-13) при $\mathrm{Nu}\to \mathrm{\infty },\mathrm{Sh}\to \mathrm{\infty }$ и выберем $n$ узловых точек $r_1,r_2,\dots ,r_n=1$. Предположим теперь, что пространственный дифференциальный оператор в уравнении (Р16-11) заменяется в узловой точке $r_j$ некоторой линейной комбинацией значений решения в узловых точках, а именно
$$\frac{{\partial }^{2}y}{\partial r^2}+{\frac{a}{r}\frac{\partial y}{\partial r}|}_{r=r_j}\sim \sum _{i=1}^n{A}_{ij}y(r_i,t).$$

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

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

Используя аппроксимацию (6.4.25) для $у$ и $\mathrm{\Theta }$, мы получаем из исходных уравнений (P16-11), (P16-12) систему обыкновенных дифференциальных уравнений для функций $y_j(t)=y(r_j,t)$, ${\mathrm{\Theta }}_j=\mathrm{\Theta }(r_j,t),j=1,2\dots ,n-1$ :
$$\begin{array}{c}\mathrm{Lw}\frac{d{y}_j}{dt}=\sum _{i=1}^n{A}_{ij}{y}_i-{\mathrm{\Phi }}^2{y}_j\exp \frac{{\mathrm{\Theta }}_j}{1+{\mathrm{\Theta }}_j/\gamma },\\ \frac{d{\mathrm{\Theta }}_j}{dt}=\sum _{i=1}^n{A}_{ij}{\mathrm{\Theta }}_i+\gamma \beta {\mathrm{\Phi }}^2{y}_j\exp \frac{{\mathrm{\Theta }}_j}{1+{\mathrm{\Theta }}_j/\gamma }\end{array}$$
(напомним, что $y_n(t)\equiv 1$ и ${\mathrm{\Theta }}_n(t)\equiv 0$ ).
Положим теперь $n=2$, т. е. выберем лишь одну внутреннюю точку $r_1$, а также точку $r_2=1$. Система уравнений (6.4.26) для $j=1$ приобретает вид
$$A_{21}+A_{11}=0,A_{21}+r_1^2{A}_{11}=2(1+a)$$

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

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

Для стационарного решения системы (6.4.30) имеет место соотношение
$${\mathrm{\Theta }}_1=\gamma \beta (1-y_1),$$

аналогичное формуле (P16-15). Подставляя его во второе из уравнений (6.4.30) (при ${\mathrm{\Theta }}_1=0$ ), получаем
$$\rho (y_1-1)-{\mathrm{\Phi }}^2{y}_1\exp \frac{\gamma \beta (1-y_1)}{1+\beta (1-y_1)}=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].