Решение $\mathbf{x}(t)$ системы обыкновенных дифференциальных уравнений
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}(t, \mathbf{x})
\]

с начальным условием $\mathbf{x}(0)=\mathbf{x}^{0}$ (за редкими исключениями, когда нам удается решить уравнение (5.7.1) аналитически) обычно приходится находить численными методами. За последние 40 лет разработан целый ряд таких методов, позволяющих строить решения задач Коши для обыкновенных дифференциальных уравнений, и в настоящее время соответствующие программы являются частью математического обеспечения практически любой ЭВМ.

Подчеркнем, что некритическое использование этих программ может в некоторых задачах приводить к неправильным результатам (или вовсе не приводить к результатам). Это может случиться, в частности, при интегрировании систем с сильно неустойчивыми траекториями или систем, содержащих малые параметры.

В этом параграфе будет приведен лишь краткий обзор указанных методов и рассмотрены некоторые практические аспекты их использования. Читателей, которых данная проблематика заинтересует более глубоко, мы отсылаем к обширной библиографии по этому вопросу (см., например, [5.7], [16*], [18*], $[19 *])$.

5.7.1. Одношаговые методы

Общая особенность численных методов решения задачи Коши (или методов численного интегрирования) состоит в том, что решение ищется в виде некоторой дискретно определенной

функции, заданной на сетке, состоящей из узлов $t_{0}=0, t_{1}$, $t_{2}, \ldots$ с шагом $h_{j}=t_{j+1}-t_{i}>0$. Методы интегрирования можно разделить на две группы: одношаговые и многошаговые.

В одношаговых методах для нахождения приближенного решения в точке $t_{i+1}$ используется аппроксимация решения лишь в одной предшествующей узловой точке $t_{j}$, т. е.
\[
\mathbf{x}^{j+1}=\mathbf{x}^{j}+h_{j} \Phi\left(t_{j}, \mathbf{x}^{j}, h_{j}\right) .
\]

Вид функции Ф зависит от конкретного задания правых частей дифференциальных уравнений (5.7.1) и от выбранного метода. Начиная с $\mathbf{x}^{0}=\mathbf{x}(0)$ при вычислениях можно использовать рекуррентную процедуру (5.7.2), продвигаясь шаг за шагом по оси $t$. Для так называемых неявных методов функция Ф зависит также от вектора $\mathbf{x}^{j+1}$, так что алгоритм перехода от $t_{j}$ к $t_{j+1}$ должен включать в себя какую-либо итерационную процедуру для решения системы нелинейных уравнений относительно компонент $\mathbf{x}^{j+1}$.

5.7.1.1. Методы с использованием разложения Тейлора

Эти методы используются лишь в случаях, когда правые части уравнений (5.7.1) можно легко продифференцировать аналитически и когда размерность системы не слишком велика. Соотношение (5.7.2) при этом принимает вид
\[
\mathbf{x}^{j+1}=\mathbf{x}^{\prime}+h_{j}\left(\mathbf{x}^{j}\right)^{\prime}+\frac{h_{l}^{2}}{2 !}\left(\mathbf{x}^{l}\right)^{\prime \prime}+\frac{h_{j}^{3}}{3 !}\left(\mathbf{x}^{j}\right)^{\prime \prime \prime}+\ldots+\frac{h_{j}^{p}}{p !}\left(\mathbf{x}^{j}\right)^{(p)},
\]

где использовано обозначение $\left(\mathbf{x}^{j}\right)^{\prime}=\mathbf{f}\left(t_{i}, \mathbf{x}^{i}\right)$, а последующие производные получаются путем дифференцирования уравнений (5.7.1) по переменной $t$ и последующей подстановки значений $t_{i}, \mathbf{x}^{i}$. Так, для $\left(\mathbf{x}^{i}\right)^{\prime \prime}$ имеем
\[
\left(\mathbf{x}^{j}\right)^{\prime \prime}=\frac{\partial \mathrm{f}\left(t_{j}, \mathbf{x}^{j}\right)}{\partial t}+\frac{\partial \mathrm{f}\left(t_{i}, \mathbf{x}^{i}\right)}{\partial \mathbf{x}}\left(\mathbf{x}^{j}\right)^{\prime} .
\]

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

Если обозначить через $\tilde{\mathbf{x}}(t)$ точное решение дифференциальных уравнений (5.7.1), удовлетворяющее начальному условию $\mathbf{x}\left(t_{j}\right)=\mathbf{x}^{i}$, то значение $\tilde{\mathbf{x}}\left(t_{j+1}\right)$ будет отличаться от $\mathbf{x}^{j+1}$, найденного по формуле (5.7.3), на остаточный член тейлоровского

разложения, т. е. на член вида ( $h_{j} h$ )
\[
\frac{h^{p+1}}{(p+1) !} \tilde{\mathbf{x}}^{(p+1)}(\tilde{t}) .
\]

Эта разность называется погрешностью аппроксимации на одном шаге, или, иначе, локальной погрешностью аппроксимации. Обозначим теперь через $\mathbf{x}(t)$ решение уравнения (5.7.1), удовлетворяющее условию $\mathbf{x}\left(t_{0}\right)=\mathbf{x}^{0}$. Тогда глобальной погрешностью аппроксимации при фиксированном $t=\bar{t}$ мы называем выражение
\[
E=\left\|\mathbf{x}^{j}-\mathbf{x}(\bar{t})\right\|,
\]

где $j=\left(\bar{t}-t_{0}\right) / h$. Если при $h \rightarrow 0$ локальная погрешность аппроксимации есть $O\left(h^{p+1}\right)$, то глобальная погрешность аппроксимации имеет порядок малости на 1 меньше $\left(E=O\left(h^{p}\right)\right.$ ). В связи с этим схема (5.7.3) называется схемой $p$-го порядка. Ниже для всех численных схем мы будем указывать только порядок глобальной погрешности аппроксимации для $h=$ const.
5.7.1.2. Варианты метода Рунге — Кутты

Наиболее распространенными из всех являются различные варианты метода Рунге — Кутты. При этом функция Ф в соотношении (5.7.2) задается с помощью расчетной схемы
\[
\mathbf{x}^{i+1}=\mathbf{x}^{i}+\gamma_{1} \mathbf{k}_{1}+\gamma_{2} \mathbf{k}_{2}+\ldots+\gamma_{m+1} \mathbf{k}_{m+1},
\]

где векторы $\mathbf{k}_{i}$ вычисляются рекуррентно по формулам
\[
\begin{array}{l}
\mathbf{k}_{1}=h \mathbf{f}\left(t_{j}, \mathbf{x}^{j}\right), \\
\mathbf{k}_{2}=h \mathbf{f}\left(t_{j}+\alpha_{1} h, \quad \mathbf{x}^{j}+\beta_{11} \mathbf{k}_{1}\right), \\
\mathbf{k}_{3}=h \mathbf{f}\left(t_{j}+\alpha_{2} h, \quad \mathbf{x}^{j}+\beta_{21} \mathbf{k}_{1}+\beta_{22} \mathbf{k}_{2}\right), \ldots \ldots . \cdots \\
\cdots \cdots \\
\mathbf{k}_{m+1}=h \mathbf{f}\left(t_{j}+\alpha_{m} h, \quad \mathbf{x}^{i}+\beta_{m 1} \mathbf{k}_{1}+\ldots+\beta_{m m} \mathbf{k}_{m}\right) .
\end{array}
\]

Для наглядности сопоставим соотношениям (5.7.5), (5.7.6) следующую схему:

Таблица 5.18. Варианты метода Рунге — Кутты
a) усовершенствованный метод Эйлера, $O\left(h^{2}\right)$
b) модифицированный метод Эйлера, $O\left(h^{2}\right)$
c) метод Хюна, $O\left(h^{3}\right)$
d) метод Кутты, $O\left(h^{3}\right)$
е) стандартный вариант метода, $O\left(h^{4}\right)$
f) метод Мерсона, $O\left(h^{4}\right)$
g) метод Нюстрёма, $O\left(h^{5}\right)$

В табл. 5.18 представлены схемы различных вариантов метода Рунге — Кутты, соответствующие разным значениям коэффициентов $\alpha, \beta$ и $\gamma$. Для каждого из них указывается также порядок погрешности аппроксимации (см. п. 5.7.1.1). Наиболее простым вариантом метода Рунге — Кутты является метод Эйлера
\[
\mathbf{x}^{j+1}=\mathbf{x}^{j}+h \mathbf{f}\left(t_{j}, \mathbf{x}^{j}\right),
\]

имеющий первый порядок точности (одновременно он является методом, основанным на разложении Тейлора). Чтобы получить у метода Рунге — Кутты порядок погрешности $O\left(h^{p}\right)$ для $p=1,2,3,4$, необходимо выполнить соответственно $1,2,3,4$ подстановки в правую часть решаемого дифференциального уравнения. При $p=5$ требуется уже минимум шесть подстановок; вообще говоря, для $p>4$ всегда оказывается необходимым использовать больше чем $p$ подстановок. Методы, имеющие порядок погрешности больше 4, используются поэтому достаточно редко — их преимущества проявляются лишь при высоких требованиях к точности вычислений.

5.7.1.3. Апостериорная оценка погрешности аппроксимации

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

ностей. Для решения задачи Коши $\mathbf{x}(t)$ при фиксированном $t$ (см. формулу (5.7.4)) погрешность аппроксимацни зависит от:
a) характера исследуемой задачи (вида дифференциальных уравнений, начальных условий и т. д.),
б) применяемого метода и его порядка,
в) используемого шага $h$ — д данном случае мы считаем его постоянным.

Если же мы фиксировали задачу и метод расчета, то в таком случае погрешность аппроксимации решения (при фиксированном $t$ ) зависит только от величины шага $h$. Предположим, что
\[
\lim _{h \rightarrow+0} \frac{E(h)}{h^{p}}=C
eq 0 .
\]

Более точно, предположим, что
\[
E(h)=C h^{p}+O\left(h^{p+1}\right) .
\]

Здесь член $O\left(h^{p+1}\right)$ убывает быстрее, чем предыдущий, и при достаточно малых $h$ мы можем им пренебречь. Подсчитаем две различных аппроксимации решения $Q_{1}$ и $Q_{2}$ (при одном $t$ ) для двух различных шагов $h_{1}$ и $h_{2}$, полагая, например, $h_{1} / h_{2}=2$ $\left(Q_{s}\right.$ здесь играют роль переменной $\mathbf{x}^{i}$ в (5.7.4), однако для различных шагов индекс $j$ также будет различным). Обозначим через $Q$ точное решение задачи. Для вычислений с шагом $h_{1}$ мы имеем
\[
Q_{1} \approx Q+C h_{1}^{p},
\]

а в случае шага $h_{2}$
\[
Q_{2} \approx Q+C h_{2}^{p} .
\]

Если рассматривать уравнения (5.7.10) и (5.7.11) как точные, то они составляют систему двух линейных алгебраических уравнений относительно неизвестных $Q$ и $C$. Значение $Q$, найденное из этой системы, обозначим $Q_{12}$ :
\[
Q_{12}=\frac{\left(\frac{h_{1}}{h_{2}}\right)^{p} Q_{2}-Q_{1}}{\left(\frac{h_{1}}{h_{2}}\right)^{p}-1} .
\]

Величина $Q_{12}$ не равна в точности $Q$, поскольку мы пренебрегли слагаемым $O\left(h^{p+1}\right.$ ) в формуле (5.7.9). Соотношение (5.7.12) и связанный с ним подход называются экстраполяцией Pичардсона и используются для апостериорной оценки погрешности аппроксимации:
\[
E_{1}=\left\|Q_{1}-Q_{12}\right\|, \quad E_{2}=\left\|Q_{2}-Q_{12}\right\| .
\]

Если принять максимально допустимую погрешность аппоксимации равной $\varepsilon$, то мы можем выбрать решение $Q_{1}$, если $E_{1}<\varepsilon$, или решение $Q_{2}$, если $E_{2}<\varepsilon$. Если же не выполнено ни одно из этих условий, то нам придется повторить вычисления для меньшего значения шага $h_{3}$, выбираемого в соответствии с оценкой
\[
h_{3} \leqslant h_{1}\left(\frac{\varepsilon}{E_{1}}\right)^{1 / p} .
\]

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

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

Очень часто случается, что решение $\mathbf{x}(t)$ в разных областях изменения независимой переменной $t$ обладает различными свойствами и, следовательно, для того чтобы выдержать предписанную точность, в каждой из этих областей нужно было бы использовать выбранный метод численного интегрирования с длиной шага, зависящей от области. Наиболее простой путь заключается в том, чтобы для всего промежутка интегрирования выбрать некую минимальную величину шага (найдя ее с помощью апостериорной оценки погрешности аппроксимации на всем рассматриваемом промежутке). Однако в большинстве областей интегрирование будет проводиться с излишней точностью и, следовательно, окажется неэкономичным. Одношаговые методы (см. (5.7.2)) позволяют нужным образом изменять величину шага интегрирования в зависимости от свойств решения. Простой алгоритм автоматического изменения шага можно получить на основе апостериорной оценки погрешности аппроксимации, вычисленной на каждом шаге интегрирования, однако такой путь неэффективен. Впервые эффективная схема вычислений с автоматическим изменением шага была разработана Мерсоном (см. табл. 5.18). В случае метода четвертого порядка нам требуется проводить пять вычислений правых частей дифференциальных уравнений, т. е. на одно больше, чем в случае стандартного метода четвертого порядка. Этот «избыток информации» используется для оценки погрешности на данном шаге:
\[
E^{\prime}=\left\|\left(2 \mathbf{k}_{1}-9 \mathbf{k}_{3}+8 \mathbf{k}_{4}-\mathbf{k}_{5}\right) / 30\right\| .
\]

Пусть $\varepsilon^{\prime}$-максимальная допустимая погрешность на одном шаге. Если $E^{\prime}<\varepsilon^{\prime}$, то полученный результат считается

приемлемым, если же $E^{\prime}>\boldsymbol{\varepsilon}^{\prime}$, то результат считается неприемлемым (и алгоритм не увеличивает $t$ ). В обоих случаях шаг для дальнейшего. интегрирования подсчитывается по формуле
\[
h^{\text {новый }}=\omega h^{\text {старый }}\left(\frac{\boldsymbol{\varepsilon}^{\prime}}{E^{\prime}}\right)^{0,2} .
\]

которая соответствует оценке (5.7.14) при $p=5$ (=1+ порядок метода). Значение параметра «безопасности» $\omega$ выбирается между 0,8 и 0,9 . Тем самым мы предупреждаем ситуации, при которых выбор шага $h=h^{\text {новый }}$ не обеспечивает предписанную заранее точность. Фактическая погрешность на всем промежутке интегрирования тесно связана с величиной $\varepsilon$ в тех случаях, когда дифференциальное уравнение не слишком «увеличивает» величину соответствующих погрешностей, т. е. когда погрешность на $M$ шагах приближенно равняется $M$-кратной погрешности, полученной на одном шаге.

Целую группу методов, аналогичных методу Мерсона, разработал в последнее время Фелберг. При этом правые части системы приходится вычислять большее число раз, чем минимально необходимо для схемы данного порядка.

5.7.1.5. Погрешности аппроксимации и погрешности округления
В рассматривавшихся до сих пор задачах мы пренебрегали погрешностями округления; при этом неявно предполагалось, что по порядку они меньше, чем погрешности аппроксимации. Это предположение несправедливо в случае повышенных требований к точности решения, т. е. при очень малых значениях шага $h$. Рассмотрим теперь схематически зависимости этих погрешностей от величины шага $h$. Если мы выбираем в качестве некоторой характерной величины погрешностей округления погрешности, возникающие при вычислениях на одном шаге метода интегрирования, а число шагов (при фиксированном $t$ ) равно $t / h$, то суммарная погрешность округления оказывается пропорциональной $h^{-1}$. 1) В то же время зависимость погрешности аппроксимации от шага $h$ определяется формулой (5.7.9). На рис. 5.21а графически изображен случай $p=1$ (метод Эйлера). Из рисунка видно, что существует нижняя граница итоговой погрешности $E^{*}$, меньше которой нельзя получить с помощью выбора шага $h$. Если мы хотим получить решение с более высокой точностью, чем $E^{*}$, то нам нужно либо повысить порядок метода

(см. рис. 5.21b), либо уменьшить погрешности округления (перейти к вычислениям с двойной точностью).

Замечание. На рис. 5.21 схематически показано поведение погрешностей при $h \rightarrow 0$. В практических расчетах обычно ис-

Рис. 5.21. Зависимость погрешности аппроксимации (1), погрешности округления (2) и полной погрешности (3) от шага $h$. $a$ ) $p=1, b$ ) $p=2$.

пользуются достаточно большие $h$-такие, чтобы можно было пренебречь погрешностями округления.

5.7.2. Многошаговые методы

Многошаговые методы в отличие от одношаговых используют значения решения и правых частей не в одной, а в нескольких предшествующих узловых точках. Линейный $k$-шаговый метод в случае постоянного шага $h$ описывается соотношением
\[
\begin{aligned}
& \alpha_{k} \mathbf{x}^{j+1}+\alpha_{k-1} \mathbf{x}^{i}+\ldots+\alpha_{0} \mathbf{x}^{j+1-k}= \\
= & h\left[\beta_{k} \mathrm{f}^{i+1}+\beta_{k-1} \mathbf{f}^{j}+\ldots+\beta_{0} \mathrm{f}^{i+1-k}\right],
\end{aligned}
\]

где $\alpha_{k}
eq 0, \alpha_{0}^{2}+\beta_{0}^{2}>0$ и использовано обозначение $\mathbf{f}^{s}=\mathbf{f}\left(t_{s}\right.$, $\mathbf{x}^{s}$ ). Существует целый ряд различных многошаговых методов. В частном случае $\alpha_{k}=1, \quad \alpha_{k-1}=-1, \quad \alpha_{i}=0$ при $i=0,1, \ldots, k-2$ соотношение (5.7.17) называется методом Адамса. При $\beta_{k}=0$ эти методы называются явными, поскольку (при известных $\mathbf{x}^{j}, \ldots, \mathbf{x}^{j-k+1}$ ) мы можем найти $\mathbf{x}^{j+1}$ из формулы (5.7.17) неитерационным способом. Такой вариант метода Адамса называют методом Адамса — Бэшфорта; его коэффициенты приведены в табл. 5.19. В случае $\beta_{k}
eq 0$ методы, описываемые соотношением (5.7.17), оказываются неявными, поскольку $\mathbf{x}^{i+1}$ входит и в правую часть формулы (5.7.17). Неявный вариант метода Адамса называют методом Адамса — Моултона, а его

Таблица 5.19. Коэффициенты метода Адамса — Бэшфорта в случае погрешности аппроксимации $p$-го порядка

Таблица 5.20. Коэффициенты метода Адамса — Моултона в случае погрешности аппроксимации $p$-го порядка

коэффициенты представлены в табл. 5.20. При этом для нахождения $\mathbf{x}^{j+1}$ (по известным $\mathbf{x}^{j}, \ldots, \mathbf{x}^{j-k+1}$ ) мы должны уже использовать какую-нибудь итерационную процедуру, для которой мы располагаем достаточно хорошим начальным приближением $\mathbf{x}^{i}$, это начальное приближение может определяться с помощью какого-либо явного метода. Так возникает метод типа предиктор — корректор.

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

ния. Такого рода алгоритм, основанный на методах Адамса Бэшфорта — Моултона и подробно разработанный Гиром, очень часто используется в приложениях.

5.7.3. Методы интегрирования «жестких» систем

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

Рассмотрим простой пример. Пусть нам нужно решить дифференциальное уравнение
\[
x^{\prime}=-1000 x, \quad x(0)=1
\]

методом Эйлера. Точное решение этой задачи имеет вид $x(t)=$ $=e^{-1000 t}$, и, следовательно, при $t>0,01$ оно практически равно нулю. Методом Эйлера (5.7.8) находим
\[
x^{j+1}=(1-1000 h) x^{j}, \quad x^{0}=1,
\]

и, следовательно,
\[
x^{i}=(1-1000 h)^{i} .
\]

Легко видеть, что при $h>0,002$ получим $\left|x^{j}\right| \rightarrow \infty$ при $j \rightarrow \infty$, и значит, для аппроксимации практически постоянной (равной нулю) функции нам нужно использовать чрезвычайно малый шаг. В случае «жестких» систем подобными свойствами обладают все явные методы интегрирования, так что их применение весьма неэффективно. При этом оказывается, что с практической точки зрения для интегрирования жестких систем удобны лишь так называемые $A$-устойчивые методы. Относительно $A$-устойчивости данного метода можно просто судить по его поведению при решении скалярного дифференциального уравнения
\[
\frac{d x}{d t}=\lambda x,
\]

где $\lambda$-в общем случае комплексное число с отрицательной вещественной частью. Численная схема, с помощью которой строится последовательность значений $x^{j}$ при постоянном шаге интегрирования $h$, называется $A$-устойчивой, если в рекуррентной формуле
\[
x^{l+1}=P(h \lambda) x^{i}
\]

множитель $P$ удовлетворяет соотношению
\[
|P(h \lambda)|<1
\]

при произвольном $h>0$. Это определение по существу представляет собой требование, чтобы $x^{j} \rightarrow 0$ при произвольном $h>0$. Указанное требование можно усилить; например, для $L$-устойчивой схемы необходимо потребовать, чтобы $|P(h \lambda)| \rightarrow$ $\rightarrow 0$ при $h \rightarrow \infty$.

За последние 15 лет был разработан целый ряд численных методов для решения жестких систем дифференциальных уравнений. Эти методы можно ориентировочно разделить на три группы. К первой относится модифицированный алгоритм Гира, основанный на многошаговой схеме типа предиктор — корректор. Этот метод используется очень широко и включается в стандартное программное обеспечение многих ЭВМ. Ко второй группе относятся так называемые полунеявные варианты метода Рунге — Кутты. В случае автономной системы
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}(\mathbf{x})
\]

одна такая группа методов описывается соотношениями (ср. формулы (5.7.5), (5.7.6))
\[
\begin{array}{l}
\mathbf{k}_{1}=h_{j}\left[\mathbf{I}-h_{j} a_{1} \mathbf{J}\left(\mathbf{x}^{j}\right)\right]^{-1} \mathbf{f}\left(\mathbf{x}^{j}\right), \\
\mathbf{k}_{2}=h_{j}\left[\mathbf{I}-h_{j} a_{2} \mathbf{J}\left(\mathbf{x}^{j}+c_{1} \mathbf{k}_{1}\right)\right]^{-1} \mathbf{f}\left(\mathbf{x}^{j}+b_{1} \mathbf{k}_{1}\right), \\
\mathbf{x}^{i+1}=\mathbf{x}^{j}+\gamma_{1} \mathbf{k}_{1}+\gamma_{2} \mathbf{k}_{2},
\end{array}
\]

где через $\mathbf{J}=\frac{d \mathbf{f}}{d \mathbf{x}}$ обозначена матрица Якоби.
Коэффициенты $a_{1}, a_{2}, b_{1}, c_{1}, \gamma_{1}$ и $\gamma_{2}$ приведены в табл. 5.21.

Таблица 5.21. Коэффициенты полунеявных вариантов метода Рунге- Кутты (5.7.24)

Все эти методы $A$-устойчивы, в чем можно убедиться непосредственной подстановкой в уравнение (5.7.20). В этих методах не нужно прибегать к итерационным процедурам, но необходимо вычислять матрицу Якоби J. Кроме того, на каждом шаге приходится дважды решать систему линейных алгебраических уравнений относительно неизвестных $\mathbf{k}_{1}$ и $\mathbf{k}_{2}$.

K третьей группе принадлежат методы, в которых используются вторые производные решения. Одним из таких методов является схема Линнигера — Уиллоуби (см. [5.31]).

5.7.4. Системы дифференциальных и алгебраических уравнений

Довольно часто встречаются системы, включающие как дифференциальные, так и алгебраические уравнения (термин «алгебраический» используется здесь для обозначения уравнения, не являющегося дифференциальным); возьмем, скажем, систему
\[
\begin{array}{l}
\frac{d \mathbf{x}}{d t}=\mathbf{f}(t, \mathbf{x}, \mathbf{u}), \\
\mathbf{0}=\mathbf{g}(t, \mathbf{x}, \mathbf{u}),
\end{array}
\]

где $\mathbf{x} \in \mathrm{R}^{n}, \mathbf{f} \in \mathrm{R}^{n}, \mathbf{u} \in \mathrm{R}^{m}, \mathbf{g} \in \mathrm{R}^{m}$. Системы такого типа возникают, к примеру, в моделях, описывающих химическую кинетику.

Для решения системы (5.7.25), (5.7.26) могут быть использованы два подхода. Первый из них состоит в интегрировании уравнений (5.7.25) с помощью какого-либо общеупотребительного метода, причем при каждом вычислении правых частей $\mathbf{f}$ решается система $m$ нелинейных уравнений (5.7.26) относительно неизвестной и. Начальная оценка и для используемой итерационной процедуры получается из предыдущего этапа вычислений правой части. Для начального условия
\[
t=0: \mathbf{x}(0)=\mathbf{x}^{0}, \quad \mathbf{u}(0)=\mathbf{u}^{0},
\]

естественно, должно иметь место соотношение
\[
\mathbf{g}\left(0, \mathbf{x}^{0}, \mathbf{u}^{0}\right)=\mathbf{0} .
\]

Второй подход заключается в преобразовании системы алгебраических уравнений в систему дифференциальных уравнений. С этой целью продифференцируем уравнения (5.7.26) по переменной $t$, и в результате мы получим систему
\[
\frac{\partial \mathbf{g}}{\partial \mathbf{u}} \frac{d \mathbf{u}}{d t}=-\frac{\partial \mathbf{g}}{\partial \mathbf{x}} \mathbf{f}-\frac{\partial \mathbf{g}}{\partial t} .
\]

Для системы $n+m$ дифференциальных уравнений (5.7.25), (5.7.29) мы имеем начальные условия (5.7.27).

При этом мы предполагали, что матрица $\partial \mathrm{g} / \partial \mathbf{u}$ является регулярной вдоль соответствующей траектории. Если же она не регулярна, то возникают трудности в обоих подходах: непрерывное решение $\mathbf{u}(t)$ системы (5.7.26) может не существовать (или быть неединственным). В таких случаях необходимо оценить полученные результаты с физической точки зрения, a также проанализировать обоснованность модели (5.7.25), (5.7.26). Рассмотренный выше случай — самый простой; более сложные случаи обсуждаются, например, в [5.34].

5.7.5. Интегрирование дифференциальных уравнений с запаздыванием

Дифференциальные уравнения с запаздыванием возникают при рассмотрении целого ряда математических моделей процессов переноса, теории автоматического регулирования, а также различного рода биологических объектов. Дифференциальные уравнения с $k$ различными запаздываниями $\tau_{1}>0, \ldots$ $\ldots, \tau_{k}>0$ можно представить в виде
\[
\frac{d \mathbf{x}}{d t}=\mathbf{f}\left(t, \mathbf{x}(t), \mathbf{x}\left(t-\tau_{1}\right), \ldots, \mathbf{x}\left(t-\tau_{k}\right)\right) .
\]

Обозначим теперь $\tau=\max \tau_{i}$. Тогда для интегрирования уравнения (5.7.30) при $t>0$ нам потребуются начальные условия на промежутке $t \in[-\tau, 0]:$ мы должны знать
\[
\mathbf{x}(t)=\boldsymbol{\varphi}(t), \quad t \in[-\tau, 0] .
\]

Итак, в отличие от систем без запаздывания, для которых достаточно знать вектор начальных условий в один момент времени $t$, в данном случае для того, чтобы однозначно определить решение, мы должны его знать на целом промежутке изменения $t$. При этом изменение в методе интегрирования состоит лишь в том, что для вычисления правой части нужно использовать значения решения в точках $t-\tau_{i}$. Эти значения получаются с помощью интерполяции из таблицы значений $(t, \mathbf{x}(t)$ ), систематически накапливаемых в ходе интегрирования. Такая таблица (память) заполняется сначала на основе знания начальных условий (5.7.31), а затем в процессе интегрирования постепенно обновляется (слишком старые значения изымаются из таблицы). Если мы используем метод интегрирования с автоматическим выбором шага, то указанная таблица оказывается неравномерной. А именно, в областях изменения $t$, где $\mathbf{x}$ меняется быстрее, таблица оказывается более плотной, это

обеспечивает сохранение точности в процессе интерполяции. Примером системы с временны́м запаздыванием служит задача 6, в которой функция роста задается формулой (P6-4).