Критерии устойчивости систем обыкновенных дифференциальных уравнений

Аннотация

Представлены критерии устойчивости в смысле Ляпунова систем обыкновенных дифференциальных уравнений, основанные на преобразованиях разностных схем. Целью преобразований является получение зависимости величины возмущения решения в произвольный момент времени от возмущения начальных данных.

Полный текст

1. Введение. Анализ устойчивости в смысле Ляпунову систем обыкновенных дифференциальных уравнений необходимо выполнять в различных областях современной науки и техники, в частности, при управлении летательными аппаратами, технологическими процессами, создании высотных строительных конструкций, в синергетике (см. [3]), при оценке динамики загрязнения экосистемы водного объекта (см. [4]). Необходимость выполнять анализ устойчивости в режиме реального времени требует разработки критериев, допускающих компьютерную реализацию.

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

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

$\frac{dY}{dt}\mathrm{<mo>=</mo>}A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}Y\mathrm{,}Y\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{Y}_0\mathrm{.}$ (1)

Предполагается, что для (1) выполнены все условия существования и единственности решения в области

$R_1\mathrm{<mo>=</mo><mo>\{</mo>}{t}_0\le t\mathrm{<mo><</mo>}\infty \mathrm{<mo>;</mo>}\tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo>}\parallel {\tilde{Y}}_0-Y_0\parallel \le {\delta }_1\mathrm{,}{\delta }_1\mathrm{<mo>></mo>0<mo>\}</mo>.}$

Элементы $\mathrm{<mo\; stretchy="false">(</mo>}n\times n\mathrm{<mo\; stretchy="false">)</mo>}$ -матрицы $A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ $—$ функции $a_{kj}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$, непрерывно дифференцируемые в $R_1$. Ниже используются каноническая норма матрицы и согласованная с ней норма вектора:

$\parallel A\parallel \mathrm{<mo>=</mo>}\underset{1\le k\le n}{{\displaystyle \max }}{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^n}\mathrm{<mo>|</mo>}{a}_{kj}\mathrm{<mo>|</mo>,}\parallel Y\parallel \mathrm{<mo>=</mo>}\underset{1\le k\le n}{{\displaystyle \max }}\mathrm{<mo>|</mo>}{y}_k\mathrm{<mo>|</mo>.}$

Для произвольно фиксированного значения $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$ всюду ниже предполагается, что величины $t$, $i$, $h$ связаны соотношениями $t\mathrm{<mo>=</mo>}{t}_{i+1}$, $h\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{t}_{i+1}-t_0\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}$, $i\mathrm{<mo>=</mo>0,1,}\dots$.

Для простоты изложения и наглядности преобразования разностных схем в качестве приближенного метода решения системы (1) используется метод Эйлера. Точное значение величины возмущения решения системы (1) методом Эйлера в форме с остаточным членом на произвольном промежутке $\mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}t\mathrm{<mo\; stretchy="false">]</mo>}$ определяется из соотношения

${\tilde{Y}}_{i+1}-Y_{i+1}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_i-Y_i\mathrm{<mo\; stretchy="false">)</mo>}+Q_{E\mathrm{,}i}\mathrm{.}$ (2)

Таким образом, величина возмущения на текущем шаге выражается через величину возмущения на предыдущем шаге. Выражая по аналогии величину возмущения на $i$ -м шаге через величину возмущения на $\mathrm{<mo\; stretchy="false">(</mo>}i-\mathrm{1<mo\; stretchy="false">)</mo>}$ -м шаге и подставляя в соотношение (2), имеем

${\tilde{Y}}_{i+1}-Y_{i+1}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-1}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\times \mathrm{<mo\; stretchy="false">(</mo>}{\tilde{Y}}_{i-1}-Y_{i-1}\mathrm{<mo\; stretchy="false">)</mo>}+\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{Q}_{E\mathrm{,}i-1}+Q_{E\mathrm{,}i}\mathrm{.}$ (3)

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

${\tilde{Y}}_{i+1}-Y_{i+1}\mathrm{<mo>=</mo>}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>}+L_i\mathrm{,}$ (4)

где

$L_i\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^i}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r}}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{Q}_{E\mathrm{<mo\; stretchy="false">(</mo>}r-\mathrm{1<mo\; stretchy="false">)</mo>}}+Q_{E\mathrm{,}i}\mathrm{,}\underset{i\to \infty }{{\displaystyle \lim }}{L}_i\mathrm{<mo>=</mo>}\overrightarrow{0}$

(см. [1]). Переходя в (4) к пределу при $h\to 0$, что равносильно $i\to \infty$, получим для всех $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$ тождество

$\tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>,}$

смысл которого в том, что для произвольного $t$ величина возмущения равна бесконечному матричному произведению, умноженному на возмущение начальных данных. Отсюда следует критерий устойчивости (5) и асимптотической устойчивости (6) системы (1) в форме необходимых и достаточных условий:

$\parallel \underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\parallel \le {\tilde{c}}_1\mathrm{,}{\tilde{c}}_1\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}$ (5)

$\underset{t\to \infty }{{\displaystyle \lim }}\parallel \underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\parallel \mathrm{<mo>=</mo>0.}$ (6)

Критерии (5), (6) позволяют определить характер устойчивости, асимптотической устойчивости либо неустойчивости систем линейных обыкновенных дифференциальных уравнений без представления решения в аналитической форме, преобразования правой части системы, построения функций Ляпунова. Мультипликативная форма выражений под знаком предела предоставляет возможность запрограммировать вычисление этих выражений и тем самым компьютеризировать анализ устойчивости.

Математически обосновано (см. [8]), что необходимая в процессе программирования замена бесконечного матричного произведения на конечное произведение сохраняет достоверность анализа устойчивости по предложенным критериям. Проведено исследование зависимости достоверности компьютерного анализа устойчивости от погрешности разностного решения системы.

Аналогичные критерии устойчивости строятся на основе методов Эйлера$–$Коши, Рунге$–$Кутта и Адамса. Использование данных критериев обеспечивает более высокую достоверность анализа устойчивости в силу улучшения оценки погрешности от отбрасывания остаточных членов.

Если матрица $A$ в (1) не зависит от времени, то критерии (5), (6), соответственно, примут вид

$\parallel \underset{m\to \infty }{{\displaystyle \lim }}{B}^{2^m}\parallel \le {\tilde{c}}_2\mathrm{,}{\tilde{c}}_2\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}$

где $B\mathrm{<mo>=</mo>}E+hA$, $i+{\mathrm{1<mo>=</mo>2}}^m$, и

$\underset{t\to \infty }{{\displaystyle \lim }}\parallel \underset{m\to \infty }{{\displaystyle \lim }}{B}^{2^m}\parallel \mathrm{<mo>=</mo>0.}$

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

В случае устойчивости (асимптотической устойчивости) системы (1) нелинейная добавка в правую часть (1) влечет ниже описываемые видоизменения в преобразования разностной схемы. Система (1) преобразуется к виду

$\frac{dY}{dt}\mathrm{<mo>=</mo>}A\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}Y+F\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>,}Y\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{Y}_0\mathrm{.}$ (7)

Предполагается, что нелинейная функция $F\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}$ определена, непрерывна и непрерывно дифференцируема на отрезке $\mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}t\mathrm{<mo\; stretchy="false">]</mo>}$ при любом выборе $t\mathrm{<mo>=</mo>}const$, $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$.

Величина возмущения определяется из соотношения

${\tilde{Y}}_{i+1}-Y_{i+1}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_i-Y_i\mathrm{<mo\; stretchy="false">)</mo>}+h\mathrm{<mo\; stretchy="false">(</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-F\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}+{\Theta }_{E\mathrm{,}i}\mathrm{,}$ (8)

где ${\Theta }_{E\mathrm{,}i}$ $—$ погрешность метода Эйлера на шаге разностной схемы.

Рекуррентное преобразование правой части (8) влечет

${\tilde{Y}}_{i+1}-Y_{i+1}\mathrm{<mo>=</mo>}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>}+D_{E\mathrm{,}i}+S_{E\mathrm{,}i}\mathrm{,}$ (9)

где

$D_{E\mathrm{,}i}\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^i}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r}}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}h\mathrm{<mo\; stretchy="false">(</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}{t}_{r-1}\mathrm{,}{\tilde{Y}}_{r-1}\mathrm{<mo\; stretchy="false">)</mo>}-F\mathrm{<mo\; stretchy="false">(</mo>}{t}_{r-1}\mathrm{,}{Y}_{r-1}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}+h\mathrm{<mo\; stretchy="false">(</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-F\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}$

$S_{E\mathrm{,}i}\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^i}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r}}b\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}{\Theta }_{E\mathrm{<mo\; stretchy="false">(</mo>}r-\mathrm{1<mo\; stretchy="false">)</mo>}}+{\Theta }_{E\mathrm{,}i}\mathrm{.}$

По аналогии с линейным случаем доказывается, что $\underset{i\to \infty }{{\displaystyle \lim }}{S}_{E\mathrm{,}i}\mathrm{<mo>=</mo>}\overrightarrow{0}$ (см. [1]).

Выполняя предельный переход в (9) при $i\to \infty$ для всех $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$ получим соотношение

$\tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>}+\underset{i\to \infty }{{\displaystyle \lim }}{D}_{E\mathrm{,}i}\mathrm{.}$

Отсюда следуют критерии устойчивости и асимптотической устойчивости системы (7):

$\parallel \tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>}\parallel \le {\tilde{c}}_3\mathrm{,}{\tilde{c}}_3\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}$

$\underset{t\to \infty }{{\displaystyle \lim }}\parallel \tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>}E+hA\mathrm{<mo\; stretchy="false">(</mo>}{t}_{i-\mathcal{l}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{Y}}_0-Y_0\mathrm{<mo\; stretchy="false">)</mo>}\parallel \mathrm{<mo>=</mo>0.}$

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

Далее рассмотрим задачу Коши для нелинейной системы

$\frac{dY}{dt}\mathrm{<mo>=</mo>}F\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>,}Y\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{Y}_0\mathrm{,}$ (10)

 где функция $F\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}$ непрерывно дифференцируема по $t$ в области

$R_2\mathrm{<mo>=</mo><mo>\{</mo>}{t}_0\le t\mathrm{<mo><</mo>}\infty \mathrm{<mo>;</mo>}\tilde{Y}\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>,}Y\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>:</mo>}\parallel {\tilde{Y}}_0-Y_0\parallel \le {\delta }_2\mathrm{,}{\delta }_2\mathrm{<mo>></mo>0<mo>\}</mo>.}$

Предполагается, что для (10) в области $R_2$ выполнены все условия существования и единственности решения.

На произвольном промежутке $\mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}t\mathrm{<mo\; stretchy="false">]</mo>}$ точное значение величины возмущения определяется из соотношения

${\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\left(1+h\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo>}}{{\tilde{y}}_{k\mathrm{,}i}-y_{k\mathrm{,}i}}\right)\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,}i}-y_{k\mathrm{,}i}\mathrm{<mo\; stretchy="false">)</mo>}+w_{k\mathrm{,}i}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (11)

На основе рекуррентного преобразования (11) имеет место соотношение

${\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}+L_{\mathrm{0,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}$ (12)

где

$L_{\mathrm{0,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^i}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}{w}_{k\mathrm{<mo\; stretchy="false">(</mo>}r-\mathrm{1<mo\; stretchy="false">)</mo>}}+w_{ki}\mathrm{,}{D}_i^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,}i}-y_{k\mathrm{,}i}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{.}$

В рассматриваемых условиях

$\underset{i\to \infty }{{\displaystyle \lim }}{L}_{\mathrm{0,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>0}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}$

(см. [5]). Критерии устойчивости и асимптотической устойчивости имеют вид

$\left|\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\right|\le {\tilde{c}}_4\mathrm{,}{\tilde{c}}_4\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$ (13)

$\underset{t\to \infty }{{\displaystyle \lim }}\left|\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\right|\mathrm{<mo>=</mo>0}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (14)

Так как

${\tilde{y}}_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-y_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}$

для всех $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$, то критерии (13), (14) можно представить в равносильной форме

$\left|\frac{{\tilde{y}}_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-y_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}{{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}}\right|\le {\tilde{c}}_4\mathrm{,}{\tilde{c}}_4\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$ (15)

$\underset{t\to \infty }{{\displaystyle \lim }}\left|\frac{{\tilde{y}}_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-y_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}{{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}}\right|\mathrm{<mo>=</mo>0}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (16)

При компьютерной реализации критериев (15), (16) требуется находить решение системы обыкновенных дифференциальных уравнений с высокой степенью точности, что необходимо для получения достоверной оценки анализа устойчивости. Поэтому целесообразно использовать методы Рунге$–$Кутта 4-го, Батчера 6-го и Дормана$–$Принса 8-го порядка и кусочно-интерполяционный метод с итерационным уточнением [2].

С целью получения критериев устойчивости в логарифмической и аддитивной форме выполним следующее преобразование под знаком предела в (13):

${\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\exp \left\{\ln {\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\right\}\text{или}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\exp \left\{{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\ln \mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\right\}\mathrm{.}$

Соответственно величина возмущения определяется из соотношения

${\tilde{y}}_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}-y_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\exp \left\{\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\ln \mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\right\}\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$

В результате критерии устойчивости и асимптотической устойчивости преобразуются к виду

$\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\ln \mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\le {\tilde{c}}_5\mathrm{,}{\tilde{c}}_5\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

$\underset{t\to \infty }{{\displaystyle \lim }}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\ln b\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}-\infty \forall k\in \overline{\mathrm{1,}n}\mathrm{.}$

С учетом того, что $h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}$ $—$ бесконечно малая и выполняется соотношение:

$\frac{\ln \mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}{h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}}\to 1\forall \mathcal{l}\le i\mathrm{,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

имеет место аддитивная форма критериев:

$\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\le {\tilde{c}}_5\mathrm{,}{\tilde{c}}_5\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

$\underset{t\to \infty }{{\displaystyle \lim }}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \sum _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}-\infty \mathrm{.}$

Далее приводится вывод критериев устойчивости системы (10) по характеру поведения правой части. Критерии должны быть равносильны (15), (16) и допускать программную реализацию.

На фиксированном промежутке $\mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}t\mathrm{<mo\; stretchy="false">]</mo>}$ значение правой части системы (10) для возмущенного и невозмущенного решения определяется из соотношений

$f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{{\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-{\tilde{y}}_{k\mathrm{,}i}-{\tilde{q}}_{k\mathrm{,}i}}{h}\mathrm{,}{f}_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{y_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,}i}-q_{k\mathrm{,}i}}{h}\mathrm{,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

где ${\tilde{q}}_{k\mathrm{,}i}$, $q_{k\mathrm{,}i}$ $—$ остаточные члены формулы Тейлора для $k$ -х компонентов решения.

Соответственно, начальное значение правой части системы (10) для возмущенного и невозмущенного решения имеет вид

$f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{{\tilde{y}}_{k\mathrm{,1}}-{\tilde{y}}_{k\mathrm{,0}}-{\tilde{q}}_{k\mathrm{,0}}}{h}\mathrm{,}{f}_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{y_{k\mathrm{,1}}-y_{k\mathrm{,0}}-q_{k\mathrm{,0}}}{h}\mathrm{,}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$

Отношение возмущения правой части (10) к начальному возмущению определяется из соотношения

$\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-{\tilde{y}}_{k\mathrm{,}i}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{y}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,}i}\mathrm{<mo\; stretchy="false">)</mo>}-w_{ki}}{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-{\tilde{y}}_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{y}_{k\mathrm{,1}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}-w_{k\mathrm{,0}}}\mathrm{,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

где $w_{k\mathrm{,}i}\mathrm{<mo>=</mo>}{\tilde{q}}_{k\mathrm{,}i}-q_{k\mathrm{,}i}$.

Сгруппируем в правой части разности возмущенного и невозмущенного решения соответствующие одной точке разностной схемы:

$\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,}i}-y_{k\mathrm{,}i}\mathrm{<mo\; stretchy="false">)</mo>}-w_{k\mathrm{,}i}}{b\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-y_{k\mathrm{,1}}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}-w_{k\mathrm{,0}}}\mathrm{,}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (17)

Преобразуем выражения в скобках числителя согласно (12). В результате получим

${\tilde{y}}_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}-y_{k\mathrm{,<mo\; stretchy="false">(</mo>}i+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-1}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-y_{k\mathrm{,1}}\mathrm{<mo\; stretchy="false">)</mo>}+L_{\mathrm{1,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}$

где

$L_{\mathrm{1,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>2}}^i}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}{w}_{k\mathrm{,<mo\; stretchy="false">(</mo>}r-\mathrm{1<mo\; stretchy="false">)</mo>}}+w_{k\mathrm{,}i}\mathrm{.}$

В отличие от (12) в качестве начальных значений возмущенного и невозмущенного решения взяты ${\tilde{y}}_{k\mathrm{,1}}$ и $y_{k\mathrm{,1}}$. По аналогии имеем

${\tilde{y}}_{k\mathrm{,}i}-y_{k\mathrm{,}i}\mathrm{<mo>=</mo>}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>1}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}+L_{2i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}$

где

$L_{\mathrm{2,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\displaystyle \sum _{r\mathrm{<mo>=</mo>1}}^{i-1}}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-r-1}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}-1}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}{w}_{k\mathrm{,<mo\; stretchy="false">(</mo>}r-\mathrm{1<mo\; stretchy="false">)</mo>}}+w_{k\mathrm{,<mo\; stretchy="false">(</mo>}i-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{.}$

В результате соотношение (17) примет вид

$\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{\tilde{Y}}_i\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_i\mathrm{,}{Y}_i\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}$

$\mathrm{<mo>=</mo>}\frac{{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-1}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-y_{k\mathrm{,1}}\mathrm{<mo\; stretchy="false">)</mo>}-{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>1}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}+L_{\mathrm{1,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}-L_{\mathrm{2,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}-w_{k\mathrm{,}i}}{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k1}-y_{k1}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k0}-y_{k0}\mathrm{<mo\; stretchy="false">)</mo>}-w_{k0}}\mathrm{,}$

для всех $k\in \overline{\mathrm{1,}n}$. Выполняя предельный переход при $i\to \infty$ с учетом того, что

$\underset{i\to \infty }{{\displaystyle \lim }}{L}_{\mathrm{1,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>0,}\underset{i\to \infty }{{\displaystyle \lim }}{L}_{\mathrm{2,}i}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>0}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{,}$

получим

$\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}\tilde{Y}\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\frac{\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-1}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-y_{k\mathrm{,1}}\mathrm{<mo\; stretchy="false">)</mo>}-\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>1}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}}{\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,1}}-y_{k\mathrm{,1}}\mathrm{<mo\; stretchy="false">)</mo>}-\mathrm{<mo\; stretchy="false">(</mo>}{\tilde{y}}_{k\mathrm{,0}}-y_{k\mathrm{,0}}\mathrm{<mo\; stretchy="false">)</mo>}}$

для всех $t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>}$ и всех $k\in \overline{\mathrm{1,}n}$. Очевидно,

$\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^{i-1}}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>1}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>.}$

Следовательно,

$\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}\tilde{Y}\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}\underset{i\to \infty }{{\displaystyle \lim }}{\displaystyle \prod _{\mathcal{l}\mathrm{<mo>=</mo>0}}^i}\mathrm{<mo\; stretchy="false">(</mo>1}+h{D}_{i-\mathcal{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}k\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (18)

Предположим, что решение системы (10) устойчиво. Но тогда с необходимостью должно выполняться условие (13), а значит, с учетом (18), и следующее условие:

$\left|\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}\tilde{Y}\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\right|\le {\tilde{c}}_6\mathrm{,}{\tilde{c}}_6\mathrm{<mo>=</mo>}const\mathrm{,}\forall t\in \mathrm{<mo\; stretchy="false">[</mo>}{t}_0\mathrm{,}\infty \mathrm{<mo\; stretchy="false">)</mo>,}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (19)

Условие (19) согласно (13) является также и достаточным условием устойчивости.

Аналогично для асимптотической устойчивости системы (10) необходимо и достаточно чтобы выполнялось следующее соотношение:

$\underset{t\to \infty }{{\displaystyle \lim }}\left|\frac{f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}\tilde{Y}\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{,}Y\mathrm{<mo\; stretchy="false">)</mo>}}{f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{\tilde{Y}}_0\mathrm{<mo\; stretchy="false">)</mo>}-f_k\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}{Y}_0\mathrm{<mo\; stretchy="false">)</mo>}}\right|\mathrm{<mo>=</mo>0}\forall k\in \overline{\mathrm{1,}n}\mathrm{.}$ (20)

Получена еще одна разновидность критериев устойчивости и асимптотической устойчивости (19), (20) систем обыкновенных дифференциальных уравнений, ориентированная на компьютерную реализацию.

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

3. Программный и численный эксперимент. Эксперимент проводился с помощью ПК на базе процессора Intel(R) Core(TM) i5-4460 в среде программирования Delphi. Написаны программы, реализующие конструкцию критериев (5), (15), (19). При анализе устойчивости линейной системы в программе циклически вычисляется частичное матричное произведение из левой части критерия (5) и через заданное количество шагов определяется и выводится на печать норма. В ходе анализа устойчивости нелинейной систем для каждого уравнения вычисляется значение выражения из левой части критериев (15), (19) и находится векторная норма.

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

Приближенные значения возмущенного и невозмущенного решения, входящие в конструкцию критериев (15), (19) находятся с помощью метода Рунге$–$Кутта 4-го порядка.

Пример 1. Исследуется на устойчивость система

${y^{\prime }}_1\mathrm{<mo>=</mo>}-y_1+\frac{y_2}{1+t^2}\mathrm{,}{y^{\prime }}_2\mathrm{<mo>=</mo>}-\frac{y_1}{1+t^2}-y_2\mathrm{.}$ (21)

Компьютерный анализ устойчивости выполняется на промежутке ${\mathrm{<mo\; stretchy="false">[</mo>0,10}}^4\mathrm{<mo\; stretchy="false">]</mo>}$ при значении шага разностной схемы ${10}^{-3}$. Результаты представлены в таблице 1. Во второй строке приводятся значения нормы, соответствующие критерию (5), в третьей строке $—$ критерию (15), в четвертой строке $—$ критерию (19).

Таблица 1. Результаты анализа устойчивости системы (21)

$t$	${10}^3$	$2\cdot {10}^3$	$5\cdot {10}^3$	$6\cdot {10}^3$	$9\cdot {10}^3$	${10}^4$
norma ${}_1$	$\mathrm{2,03}\cdot {10}^{-434}$	$\mathrm{1,03}\cdot {10}^{-868}$	$\mathrm{1,35}\cdot {10}^{-2171}$	$\mathrm{6,84}\cdot {10}^{-2606}$	$\mathrm{8,95}\cdot {10}^{-3909}$	$\mathrm{4,54}\cdot {10}^{-4343}$
norma ${}_2$	$\mathrm{6,77}\cdot {10}^{-435}$	$\mathrm{3,44}\cdot {10}^{-869}$	$\mathrm{4,49}\cdot {10}^{-2172}$	$\mathrm{2,28}\cdot {10}^{-2606}$	$\mathrm{2,98}\cdot {10}^{-3909}$	$\mathrm{1,51}\cdot {10}^{-4343}$
norma ${}_3$	$\mathrm{5,08}\cdot {10}^{-435}$	$\mathrm{2,58}\cdot {10}^{-869}$	$\mathrm{3,37}\cdot {10}^{-2172}$	$\mathrm{1,71}\cdot {10}^{-2606}$	$\mathrm{2,23}\cdot {10}^{-3909}$	$\mathrm{1,13}\cdot {10}^{-4343}$

Представленные значения нормы по всем трем критериям монотонно стремятся к нулю, что является признаком асимптотической устойчивости. Время работы программы при анализе по критерию (5) составляет около 4 с., по критериям (15), (19) $—$ около 1,5 с.

Пример 2. Исследуется на устойчивость решение системы (см. [7])

${y^{\prime }}_1\mathrm{<mo>=</mo>}-y_2\sqrt{y_1^2+y_2^2}\mathrm{,}{y^{\prime }}_2\mathrm{<mo>=</mo>}{y}_1\sqrt{y_1^2+y_2^2}\mathrm{.}$ (22)

Первоначально исследуется нулевое решение системы (22) при численных значениях шага и промежутка из предыдущего примера. Начальные значения компонент возмущенного решения ${\tilde{y}}_{10}{\mathrm{<mo>=</mo>10}}^{-5}$, ${\tilde{y}}_{20}\mathrm{<mo>=</mo>2}\cdot {10}^{-5}$. Результаты анализа устойчивости представлены в таблице 2. Во второй строке результаты анализа по критерию (15), в третьей строке $—$ по критерию (19).

Таблица 2. Результаты анализа устойчивости нулевого решения системы (22)

$t$	${10}^3$	$2\cdot {10}^3$	$5\cdot {10}^3$	$6\cdot {10}^3$	$9\cdot {10}^3$	10 ${}^4$
norma ${}_1$	1,01	1,02	1,05	1,06	1,08	1,09
norma ${}_2$	1,01	1,02	1,05	1,06	1,08	1,09

Значения нормы ограничены константой, что свидетельствует об устойчивости.

Далее выполняется анализ устойчивости ненулевого решения системы (22): $y_{10}\mathrm{<mo>=</mo>}{y}_{20}\mathrm{<mo>=</mo>0,01}$, величина возмущения начальных данных $\delta {\mathrm{<mo>=</mo>10}}^{-5}$.

Наблюдается монотонный рост значений нормы по обоим критериям, что свидетельствует о неустойчивости. Время работы программы при анализе нулевого решения около $\mathrm{1,5}$ с, ненулевого $—$ около $3$ с. Полученные результаты анализа устойчивости системы (22) находятся в полном соответствии с представленными в [7].

Пример 3. Исследуется модель периодической реакции Белоусова$–$Жаботинского (см. [6])

${y^{\prime }}_1\mathrm{<mo>=</mo>77,27<mo\; stretchy="false">(</mo>}{y}_2+y_1\mathrm{<mo\; stretchy="false">(</mo>1}-\mathrm{8,375}\cdot {10}^{-6}{y}_1-y_2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}{y^{\prime }}_2{\mathrm{<mo>=</mo>77,27}}^{-1}\mathrm{<mo\; stretchy="false">(</mo>}{y}_3-y_2\mathrm{<mo\; stretchy="false">(</mo>1}+y_1\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>,}{y^{\prime }}_3\mathrm{<mo>=</mo>0,161<mo\; stretchy="false">(</mo>}{y}_1-y_3\mathrm{<mo\; stretchy="false">)</mo>.}$ (23)

Первоначально анализ устойчивости выполняется при начальных условиях $y_1\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1}$, $y_2\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>2}$, $y_3\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>3}$ на промежутке ${\mathrm{<mo\; stretchy="false">[</mo>0,10}}^3\mathrm{<mo\; stretchy="false">]</mo>}$ с шагом ${10}^{-5}$, величина возмущения для каждого решения соответственно равна ${\delta }_1{\mathrm{<mo>=</mo>10}}^{-5}$, ${\delta }_2\mathrm{<mo>=</mo>2}\cdot {10}^{-5}$, ${\delta }_3\mathrm{<mo>=</mo>3}\cdot {10}^{-5}$.

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

Таблица 3. Результаты анализа устойчивости ненулевого решения системы (22)

$t$	${10}^3$	$2\cdot {10}^3$	$5\cdot {10}^3$	$6\cdot {10}^3$	$9\cdot {10}^3$	${10}^4$
norma ${}_1$	14,99	28,63	73,92	89,84	139,39	156,46
norma ${}_2$	8,00	14,83	36,49	44,46	69,25	77,79

Далее выполняется анализ устойчивости системы (23) при начальных условиях $y_1\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>4}$, $y_2\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1,1}$, $y_3\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>4}$. В отличие от предыдущего случая резкий скачок значений нормы происходит в самом начале промежутка исследования системы (23) (таблица 3).

Таблица 4. Результаты анализа устойчивости решения системы (23) при начальных условиях $y_1\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1}$, $y_2\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>2}$, $y_3\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>3}$

$t$	$18$	$19$	$20$	$21$	$22$	$23$
norma ${}_1$	$\mathrm{61,79}$	$\mathrm{4175,04}$	$\mathrm{3,92}\cdot {10}^3$	$\mathrm{8,59}\cdot {10}^5$	$\mathrm{9,40}\cdot {10}^5$	$\mathrm{3,83}\cdot {10}^6$
norma ${}_2$	$\mathrm{26,13}$	$\mathrm{85,55}$	$\mathrm{2,22}\cdot {10}^3$	$\mathrm{7,53}\cdot {10}^5$	$\mathrm{6,88}\cdot {10}^5$	$\mathrm{2,03}\cdot {10}^6$
$t$	$24$	$25$	$26$	$27$	$28$	$29$
norma ${}_1$	$\mathrm{6,29}\cdot {10}^4$	$\mathrm{5,36}\cdot {10}^4$	$\mathrm{4,56}\cdot {10}^4$	$\mathrm{3,88}\cdot {10}^4$	$\mathrm{3,30}\cdot {10}^4$	$\mathrm{2,80}\cdot {10}^4$
norma ${}_2$	$\mathrm{1,18}\cdot {10}^5$	$\mathrm{1,01}\cdot {10}^5$	$\mathrm{8,55}\cdot {10}^4$	$\mathrm{7,27}\cdot {10}^4$	$\mathrm{6,17}\cdot {10}^4$	$\mathrm{5,24}\cdot {10}^4$

Таблица 5. Результаты анализа устойчивости решения системы (23) при начальных условиях $y_1\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>4}$, $y_2\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>1,1}$, $y_3\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo><mo>=</mo>4}$

$t$	$1$	$2$	$3$	$4$	$5$	$6$
norma ${}_1$	$\mathrm{3,21}\cdot {10}^4$	$\mathrm{3,68}\cdot {10}^5$	$\mathrm{4,31}\cdot {10}^5$	$\mathrm{3,03}\cdot {10}^4$	$\mathrm{2,57}\cdot {10}^4$	$\mathrm{2,19}\cdot {10}^4$
norma ${}_2$	$\mathrm{1,62}\cdot {10}^4$	$\mathrm{3,10}\cdot {10}^5$	$\mathrm{2,97}\cdot {10}^5$	$\mathrm{4,69}\cdot {10}^4$	$\mathrm{3,99}\cdot {10}^4$	$\mathrm{3,40}\cdot {10}^4$
$t$	$7$	$8$	$9$	$10$	$11$	$12$
norma ${}_1$	$\mathrm{1,86}\cdot {10}^4$	$\mathrm{1,58}\cdot {10}^4$	$\mathrm{1,35}\cdot {10}^4$	$\mathrm{1,14}\cdot {10}^4$	$\mathrm{9,78}\cdot {10}^3$	$\mathrm{8,32}\cdot {10}^3$
norma ${}_2$	$\mathrm{2,89}\cdot {10}^4$	$\mathrm{2,46}\cdot {10}^4$	$\mathrm{2,09}\cdot {10}^4$	$\mathrm{1,78}\cdot {10}^4$	$\mathrm{1,51}\cdot {10}^4$	$\mathrm{1,28}\cdot {10}^4$

В целом наблюдаются ограниченные колебания значений нормы, что в соответствии с критериями (15), (19) свидетельствует об устойчивости решения системы (23). Время работы программы в обоих случаях около $7.5$ мин. При анализе устойчивости жестких нелинейных систем нередко возникают трудности при вычислении возмущенного и невозмущенного решения, входящих в конструкцию критериев. Для их вычисления целесообразно использовать специализированные программы для решения жестких систем (см. [6]) или кусочно-интерполяционный метод, основанный на приближении решения и правой части системы полиномами Лагранжа с числовыми коэффициентами (см. [8]). Приближения решения и правой части системы находятся с высокой точностью, кроме того существенно сокращается время на вычисление требуемых приближений. На практике это дает возможность проводить анализ устойчивости на промежутках существенно большей длины и устанавливать асимптотические свойства исследуемых систем обыкновенных дифференциальных уравнений.

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