Динамическая модель с сосредоточенными параметрами может быть представлена в виде системы обыкновенных дифференциальных уравнений
\[
\frac{d x_{i}}{d t}=f_{i}\left(x_{1}, x_{2}, \ldots, x_{n}, \alpha_{1}, \ldots, \alpha_{m}\right), \quad i=1,2, \ldots, n .
\]

Систему (5.1.1) можно записать также в векторной форме
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}(\mathbf{x}, \boldsymbol{\alpha}),
\]

где использованы обозначения
\[
\mathbf{x}=\left(x_{1}, x_{2}, \ldots, x_{n}\right), \quad \mathbf{f}=\left(f_{1}, f_{2}, \ldots, f_{n}\right), \quad \boldsymbol{\alpha}=\left(\alpha_{1}, \alpha_{2}, \ldots, \alpha_{m}\right) .
\]

Здесь $x_{i}$ – это переменные состояния, $t$-время, а $\boldsymbol{\alpha}$-вектор: параметров системы. В большинстве случаев мы будем рассматривать свойства системы в зависимости от одного скалярного параметра ( $m=1$ ). Это означает, что значения оставшихся параметров фиксированы и, следовательно, такие параметры являются «составной частью» функции f. Кроме того, мы будем считать систему автономной, т. е. будем считать, что время $t$ не входит явным образом в правые части уравнений (5.1.1). Неавтономным задачам будет посвящен § 5.11. Мы предполагаем, что правые части системы (5.1.1) представляют собой непрерывно дифференцируемые функции своих аргументов.

Стационарное решение уравнений (5.1.1) удовлетворяет системе нелинейных уравнений
\[
f_{i}\left(x_{1}, x_{2}, \ldots, x_{n}, \boldsymbol{\alpha}\right)=0, \quad i=1,2, \ldots, n .
\]

Для решения этой системы (при фиксированном $\alpha$ ) мы можем в принципе использовать любые методы решения систем нелинейных уравнений. Читатель, который более глубоко заинтересуется этой проблемой, может найти обзор соответствующих методов в книге [5.1]. Если функции $f_{i}$ можно продифференцировать аналитически, то удобнее всего воспользоваться методом

типа метода Ньютона, представляющим собой один из вариантов метода итераций. Для этого построим последовательность $\left\{\mathbf{x}^{k}\right\}$ в соответствии с рекуррентной формулой вида
\[
\mathbf{x}^{k+1}=\mathbf{x}^{k}+\omega \Delta \mathbf{x}^{k}, \quad k=0,1,2, \ldots,
\]

где $\Delta \mathbf{x}^{k}$ для данного $k$ находится из решения системы линейных алгебраических уравнений
\[
\mathbf{J}\left(\mathbf{x}^{k}, \boldsymbol{\alpha}\right) \Delta \mathbf{x}^{k}=-\mathbf{f}\left(\mathbf{x}^{k}\right) .
\]

Здесь через $\mathbf{J}$ обозначена матрица Якоби для системы (5.1.3), т. е. матрица $\mathbf{J}=\left[\partial f_{i} / \partial x_{j}\right]$. Параметр $\omega \in(0,1]$ выбирается обычно так, чтобы
\[
\left\|\mathbf{f}\left(\mathbf{x}^{k+1}\right)\right\|<\left\|\mathbf{f}\left(\mathbf{x}^{k}\right)\right\| .
\]

В окрестности решения системы 5.1.3 условие (5.1.6) выполняется при $\omega=1$, что соответствует классическому методу Ньютона; на большем удалении от решения иногда оказывается полезным выбирать меньшую величину $\omega$, с тем чтобы условие (5.1.6) было выполнено. Указанная модификация классического метода Ньютона позволяет расширить область сходимости метода по отношению к выбору начального приближения $\mathbf{x}^{0}$.

Метод Ньютона обладает определенными достоинствами в отношении его сходимости, однако слабым местом этого метода является необходимость дифференцирования функций $f_{i}$. Если аналитическое дифференцирование оказывается затруднительным, то часто используется вариант метода Ньютона с матрицей Якоби, приближенно вычисляемой с помощью разностных отношений
\[
\frac{\partial f_{i}(\mathbf{x})}{\partial x_{j}} \approx \frac{f_{i}\left(\mathbf{x}+h \mathbf{e}_{j}\right)-f_{i}(\mathbf{x})}{h} .
\]

где $h$ – достаточно малый шаг, а $\mathbf{e}_{i}$ – единичный вектор с единицей на $j$-м месте. Таким образом, для проведения одной итерации нам в общем случае потребуется вычислять вектор $\mathbf{f}$ $(n+1)$ раз.

Сходимость метода Ньютона (равно как и других методов) существенно зависит от выбора начального приближения решения $\mathbf{x}^{0}$. При этом уравнения (5.1.3) часто имеют не одно, а несколько решений. Тем самым возникает вопрос, как отыскать все решения данной системы. Такого рода задача разрешима, вообще говоря, лишь в некоторых специальных случаях, например при $n=1$ и для функции $f_{1}$ в виде полинома. В случае функций $f_{i}$ общего вида чаще всего выбирается несколько различных начальных приближений $\mathbf{x}^{0}$, после чего для каждого из них последовательно применяется схема Ньютона. При этом для

выбранного начального приближения мы практически получаем один из следующих результатов:

1. итерационный процесс сходится к физически приемлемому решению (новому или уже найденному ранее);
2. итерационный процесс сходится к физически неприемлемому решению;
3. итерационный процесс расходится, т. е. значения переменных с какого-то момента в ходе итераций превышают установленные пределы;
4. итерационный процесс не сходится с предписанной точностью при числе итераций, не превосходящем заданного.

Рис. 5.1. Диаграмма успешности применения метода Ньютона в зависимости от случайно выбранного начального приближения. $K$ – порядковый номер: начального приближения. На вертикальной оси приведены номера полученных решений (см. табл. 5.1 при $N=4$ ). Символ д означает расходимость метода: от соответствующего начального приближения. Для решений, отмеченных в табл. 5.1 звездочкой, два симметричных решения на рисунке представлены отдельно.

Проиллюстрируем указанный подход на примере, приведенном в гл. 4 в качестве задачи 8 для случая $N=2,3,4$. В табл. 5.1 представлены решения для двух, трех и четырех связанных между собой реакторов (т. е. для $n=4,6,8$ в системе (5.1.3)). Число решений при этом оказывается соответственно равным 3, 7 и 15 [5.2]. В случае цепочки из большего числа реакторов количество различных стационарных решений будет еще больше.

Начальные приближения $\mathbf{x}^{0}$ для метода Ньютона выбирались случайным образом из интервалов
\[
x_{i}^{0} \in(0,10), \quad i=1,2, \ldots, n .
\]

Примеры расчета итераций приведены в табл. 5.2. На рис. 5.1 изображена диаграмма успешности этого подхода для случая четырех связанных между собой реакторов $(N=4)$.

Таблица 5.1. Стационарные решения в случае модели линейного каскада связанных друг с другом реакторов (задача 8: $A=2, B=6, D_{1}=0,4$, $D_{2}=4$ ). В первой строке каждой ячейки указана величина концентрации $X$, во второй – величина концентрации $Y$.

Таблица 5.2. Поведение итераций в процессе вычисления стационарных решений методом Ньютона (задача 8: $N=4, A=2, B=6, D_{1}=0,4$, $D_{2}=4$ ).

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

Запишем систему (5.1.3) в векторной форме
\[
\mathbf{f}(\mathbf{x})=\mathbf{0}
\]

при этом значение параметра $\boldsymbol{\alpha}$ мы считаем фиксированным и поэтому в уравнениях его не указываем. Идея метода введения искусственного параметра заключается в построении функции $\mathbf{H}(\tau, \mathbf{x})$, зависящей от некоторого вещественного параметра $\tau$, причем обычно $\tau \in[0,1]$. Функция $\mathbf{H}$ выбирается таким образом, чтобы решение уравнения
\[
\mathbf{H}(0, \mathbf{x})=\mathbf{0}
\]

можно было построить по возможности просто, например, не прибегая к процессу итераций (это решение мы обозначим через $\mathbf{x}_{0}$ ). Далее, уравнение
\[
\mathbf{H}(1, \mathbf{x})=0
\]

должно иметь те же корни, что и уравнение (5.1.9), например
\[
\mathbf{H}(1, \mathbf{x})=\mathbf{f}(\mathbf{x}) .
\]

Если существует непрерывная зависимость $\mathbf{x}(\tau)$, описывающая решение уравнения
\[
\mathbf{H}(\tau, \mathbf{x}(\tau))=\mathbf{0}, \quad \tau \in[0,1],
\]

с условием $\mathbf{x}(0)=\mathbf{x}_{0}$, то $\mathbf{x}(1)$ является решением исходного уравнения (5.1.9). В принципе мы можем «протянуть» решение $\mathbf{x}$ от $\tau=0$ (когда $\mathbf{x}=\mathbf{x}_{0}$ ) до $\tau=1$ двумя способами. Первый из них заключается в последовательном использовании какого-нибудь итерационного метода (например, метода Ньютона) в узлах сетки $\tau_{0}=0, \tau_{1}, \tau_{2}, \ldots, \tau_{m}=1$. В качестве начального приближения процесса итераций для $\tau=\tau_{j}$ можно выбрать решение $\mathbf{x}\left(\tau_{j-1}\right)$. Второй способ основан на дифференцировании тождества (5.1.13) по $\tau$, т. е. зависимость $\mathbf{x}(\tau)$ ищется как решение системы дифференциальных уравнений
\[
\frac{d \mathbf{x}}{d \tau}=-\left[\mathbf{H}_{\mathbf{x}}^{\prime}(\tau, \mathbf{x})\right]^{-1} \mathbf{H}_{\tau}^{\prime}(\tau, \mathbf{x})
\]

с начальным условием $\mathbf{x}(0)=\mathbf{x}_{0}$.
Укажем две возможности выбрать функцию $\mathbf{H}$ :

и
\[
\mathbf{H}(\tau, \mathbf{x})=\mathbf{f}(\mathbf{x})+(\tau-1) \mathbf{f}\left(\mathbf{x}_{0}\right)
\]
\[
\mathbf{H}(\tau, \mathbf{x})=(1-\tau)\left(\mathbf{x}-\mathbf{x}_{0}\right)+\tau \mathbf{f}(\mathbf{x}) .
\]

Представление (5.1.15) тесно связано с методом Ньютона. (Чтобы убедиться в этом, нужно, взяв Н в виде (5.1.15), решить. дифференциальное уравнение (5.1.14) методом Эйлера с шагом $h=1$.) Более подробный анализ различных вариантов метода введения искусственного параметра читатель может найти в книгах $[5.1,5.3]$. Описание некоторых других подходов и их связь с итерационными методами рассмотрены в работе [5.4].