Рассмотрим систему

где
\[
\dot{x}=a x+b y+P(x, y), \quad \dot{y}=c x+d y+Q(x, y),
\]
\[
P(x, y)=\sum_{i+j \geqslant 2} Q_{i j} x^{i} y^{j}, \quad Q(x, y)=\sum_{i+j \geqslant 2} b_{i j} x^{i} y^{j}
\]

и коэффициенты можно считать зависящими от некоторого параметра $\lambda$ (параметров $\lambda_{i}$ ).
В начале координат будет устойчивый фокус, если
\[
q \equiv\left|\begin{array}{ll}
a & b \\
c & d
\end{array}\right|>0, \quad p \equiv-(a+d)>0, \quad p^{2}-4 q<0 .
\]

Следуя Ляпунову и Андронову $[17,2]$, строим в окрестности начала координат функцию
\[
\psi\left(\rho_{0}, p\right)=\left(e^{-\frac{2 \pi}{\omega} p}-1\right) \rho_{0}+\alpha_{2} \rho_{0}^{2}+\alpha_{3} \rho_{0}^{3}+\ldots .
\]

Теорема (Ляпунова). Первый не равный нулю коэффициент в разложении функции $\psi\left(\rho_{0}, 0\right)$ непременно нечетного номера.

Если $p=0$, то первый не равный нулю коэффициент $\alpha_{i}$ называется ляпуновской величиной: $\alpha_{3} \equiv L_{1}$ – первая ляпуновская величина; если $\alpha_{3}=0, \alpha_{5}
eq 0$, то $\alpha_{5} \equiv L_{2}$ – вторая ляпуновская величина, и т. д.

Рассмотрение функции (2) в зависимости от параметров позволяет сделать исчерпывающие заключения относительно характера траекторий в окрестности состояния равновесия $x=y=0$ при различных значениях параметров. Отличные от нуля корни функции соответствуют предельным циклам. Вычисление ляпуновских величин в общем случае системы $n$ уравнений с характеристическим уравнением, имеющим пару чисто мнимых корней и остальные корни с отрицательной действительной частью, с помощью процедуры, указанной Ляпуновым, приводится к вычислению тех же величин $\alpha_{i}$ для некоторой системы второго порядка вида (1), выводимой из предложенной системы $n$ уравнений.

Первая ляпуновская величина, выраженная через коэффициенты системы (1), приведена в $[4,8,11,13]$. Вторая ляпуновская величина вычислена в [21]. Для системы трех и четырех уравнений общего вида в [11] изложено развернутое приведение к каноническому виду и вычислена первая ляпуновская величина через коэффициенты преобразованной системы (в случае чётырех уравнений отдельно для случая, когда вторая пара корней комплексная и когда она действительная). Для системы $n$ уравнений в случае, когда разложения правых частей не содержат членов второго порядка, в [11] дано выражение первой ляпуновской величины в виде интеграла по кривым вспомогательной консервативной системы без приведения исходной системы к каноническому виду. Там же даны аналогичные выражения первой ляпуновской величины и для общего вида систем двух, трех и четырех уравнений.

Опишем применительно к системе (1) простейшие возможные случаи.
А. Пусть $p\left(\lambda_{0}\right)=0$ и $L_{1}\left(\lambda_{0}\right)<0$. При переходе через границу $p=0$ от положительных значений к отрицательным появляется единственный устойчивый предельный цикл. При обратном изменении параметра устойчивый предельный цикл стягивается в точку. Переход через границу $p=0$ соответствует возникновению области неустойчивости внутри устойчивого предельного цикла, которая, однако, остается сколь угодно малой при достаточно малом нарушении условия устойчивости и стягивается в точку при обратном изменении параметра. Практически система при малом нарушении условия устойчивости будет вести себя как устойчивая (граница области устойчивости – безопасная).
Б. Пусть $p\left(\lambda_{0}\right)=0$ и $L_{1}\left(\lambda_{0}\right)>0$. При переходе через границу $p=0$ от положительных значений к отрицательным к состоянию равновесия стягивается единственный неустойчивый предельный цикл. При обратном изменении параметра из состояния равновесия появляется неустойчивый предельный цикл. Переход через границу $p=0$ соответствует исчезновению области устойчивости внутри неустойчивого предельного цикла; изображающая точка при этом срывается с состояния равновесия и уходит за пределы рассматриваемой окрестности состояния равновесия (граница области устойчивости опасная).

Особенности в поведении динамической системы вблизи тех точек границы $p=0$, где безопасная часть границы переходит в опасную и где, следовательно, первый ляпуновский коэффициент $L_{1}\left(\lambda_{0}\right)$ обращается в нуль, определяются знаком второй ляпуновской величины $L_{2}\left(\lambda_{0}\right)$, для вычисления которой необходимо учесть в разложениях правых частей уравнения (1) члены до пятого порядка включительно. В возможностях, которые здесь возникают, можно ориентироваться, рассматривая функцию (2). Нетрудно показать [10], что если в ряду коэффициентов $p, \alpha_{3} \equiv L_{1}(\lambda), \alpha_{5} \equiv L_{2}(\lambda)$ имеются одна или две перемены знака, то в малой окрестности начала координат будут существовать один или два корня функции $\psi\left(\rho_{0}, p\right)$ и соответственно один или два предельных цикла на фазовой плоскости вокруг начала координат.

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

В общем случае системы $n$ уравнений с характеристическим уравнением, имеющим одну пару чисто мнимых корней и остальные корни с отрицательной действительной частью, все сказанное о поведении траекторий системы (1) будет
Рис. 1.

справедливо по отношению к некоторому двумерному многообразию в фазовом пространстве системы, заполненному траекториями и содержащему состояние равновесия и предельные циклы (если последние существуют). Для остальных траекторий это многообразие будет элементом притяжения. Изложенное иллюстрирует (применительно к системе трех уравнений) рис. 1. В случае $L_{1}\left(\lambda_{0}\right)<0$ при изменении знака $p$ неустойчивый предельный цикл, расположенный на двумерном многообразии, стягивается к состоянию равновесия. При этом исчезает пространственная область устойчивости, заполненная траекториями, идущими к состоянию равновесия (область стягивается к паре сепаратрис, идущих в состояние равновесия).

Рассмотрим поведение динамической системы в малой окрестности точки на границе области устойчивости, в которой первая ляпуновская величина обращается в нуль. Рис. $2 a, 6$ дают разбиения плоскости параметров (параметрами могут быть, например, сами величины $p$ и $L_{1}$ ) для случаев $L_{2}>0$ и $L_{2}<0$ в точке смыкания опасной и безопасной границ. Жирной линией отмечена безопасная часть границы $p=0$, тонкой – опасная. Пунктиром отмечена бифуркационная кривая двойных циклов. Ее точкам соответствуют двойные циклы (полуустойчивые), возникшие из сгущения траекторий. При смещении параметров с кривой двойных циклов двойной предельный цикл либо исчезает, либо разделяется на два (устойчивый и неустойчивый). Если $L_{2}>0$, то при обходе точки $p=0, L_{1}=0$, начиная с области $p>0$,
(a)
(6)

Рис. 2. Разбиение плоскости параметров в окрестности точки, где опасная часть границы переходит в безопасную.
$L_{1}<0$ в направлении против часовой стрелки, последовательно осуществляются бифуркации: рождение устойчивого предельного цикла из состояния равновесия, слияние неустойчивого предельного цикла с устойчивым, исчезновение двойного цикла, рождение неустойчивого цикла. Если $L_{2}<0$, то бифуркации осуществляются в таком порядке: рождение устойчивого, рождение неустойчивого, слияние устойчивого с неустойчивым и исчезновение двойного цикла. Штриховкой покрыта область существования двух циклов. Оценка расположения бифуркационной кривой двух циклов в зависимости от значения $L_{2}$ дана в [21]. Асимптотическое представление кривой двойных циклов (в обозначениях разд. 1) дается формулой
\[
\sqrt{q} L_{1}^{2}+8 \pi L_{2} p=0 \quad\left(L_{1} L_{2}<0\right) .
\]

Более сложные ситуации возникают, если первая не обращающаяся в нуль на границе $p=0$ ляпуновская величина будет $L_{k}$, где $k>2$. В этом случае при малых изменениях параметров $\lambda_{i}$ в окрестности состояния равновесия в начале координат может появиться более двух предельных циклов. Ситуации здесь будут похожими на уже рассмотренные при $L_{1}
eq 0$ или $L_{2}
eq 0$ на границе $p=0$ и будут зависеть от знака первой отличной от нуля ляпуновской величины и от того, сколько будет циклов (четное или нечетное число). Сложность, возникающая в этом случае, связана в первую очередь с тем, что технические трудности вычисления последовательных ляпуновских величин с ростом их номера быстро возрастают. Эту трудность можно обойти, используя ЭВМ для алгебраического вычисления ляпуновских величин.
С. Д. Щуко $[29,30,31]$ разработала алгоритм, позволяющий последовательно вычислять на ЭВМ ляпуновские величины (каждую следующую в предположении равенства нулю всех предыдущих). В качестве исходной для вычисления ляпуновских величин рассматривается каноническая система уравнений
\[
\dot{x}=-y+\sum_{m=2}^{n} P_{m}(x, y), \quad \dot{y}=x+\sum_{m=2}^{n} Q_{m}(x, y),
\]

где $P_{m}$ и $Q_{m}$ – однородные полиномы степени $m$, взятые из разложений правых частей системы (1) после приведения ее к каноническому виду, наибольшие степени которых, $n$, согласуются с порядковым номером вычисляемой ляпуновской величины.

Для вычислений удобны методы, позволяющие осуществлять внутренний контроль вычислений в силу симметрии получающихся выражений. Это реализуется при переходе к комплексным переменным [14].
Полагая $z=x+i y, \bar{z}=x-i y$, получаем
\[
\frac{d \bar{z}}{d z}=-\frac{\bar{z}+Z(z, \bar{z})}{z+\bar{Z}(z, \bar{z})},
\]

где
\[
\begin{array}{l}
Z(z, \bar{z})=\sum_{m=2}^{n} \sum_{k=0}^{m} a_{k m} z^{k} \bar{z}^{m-k}, \\
\bar{Z}(z, \bar{z})=\sum_{m=2}^{n} \sum_{k=0}^{m} \bar{a}_{k m} \bar{z}^{k} z^{m-k} .
\end{array}
\]

Наличие центра в начале координат равносильно существованию голоморфного интеграла
\[
\Phi(z, \bar{z})=\sum_{m=0}^{\infty} \sum_{k=0}^{m} \Phi_{k, m-k} z^{k} \bar{z}^{m-k}
\]

уравнения
\[
L(\Phi) \equiv \frac{\partial \Phi}{\partial z} \cdot \frac{d z}{d t}+\frac{\partial \Phi}{\partial \bar{z}} \cdot \frac{d \bar{z}}{d t} .
\]

Функции $\Phi(z, \bar{z})$ сопоставляется последовательность $(m+1)$-мерных комплексных векторов $\Phi_{m}$ с компонентами, являющимися коэффициентами однородных полиномов степени $m$ из (5):
\[
\Phi_{m}=\left\{\Phi_{0, m} ; \Phi_{1, m-1} ; \ldots ; \Phi_{m-1,1} ; \Phi_{m, 0}\right\} .
\]

Для того чтобы в начале координат был цеңтр системы (3), необходимо и достаточно [18], чтобы векторное уравнение (6) имело решение $\Phi=\left\{\Phi_{m}\right\}$, удовлетворяющее условиям
\[
\begin{array}{c}
\Phi_{0,0}=\Phi_{0,1}=\Phi_{1,0}=\Phi_{0,2}=\Phi_{2,0}=0 ; \quad \Phi_{k, k}=0, \quad k \geqslant 2 ; \\
\Phi_{1,1}=1 ; \quad \Phi_{k, m-k}=\bar{\Phi}_{m-k, k} .
\end{array}
\]

Векторы $\Phi_{m}$ находятся по рекуррентной формуле
\[
\Phi_{m}=-C_{m, m}^{-1}\left\{\sum_{l=2}^{m-1} C_{m, l} \Phi_{l}\right\}, \quad m>2 .
\]

Здесь $C_{m, l}=\left\{c_{r s}^{(m, l)}\right\}$ – матрица из $m+1$ строк и $l+1$ столбцов, причем
\[
c_{r s}^{(m, l)}=(l-s) a_{r-s, m-l+1}+s \bar{a}_{m-l-r+s, m-l+1},
\]

ғде $r$ – номер строки, $s$ – номер столбца, $C_{m, m}^{1-}=\left\{c_{r s}^{-1}\right\}$ қвадратная диагональная $(m+1)$-матрица с элементами вида
\[
c_{r s}^{-1}=\left\{\begin{array}{ccc}
-\frac{i \delta_{r s}}{m-2 s} & \text { при } & m
eq 2 s \\
0 & \text { при } & m=2 s
\end{array} ; \delta_{r s}-\right.\text { символ Кронекера. }
\]

Ляпуновские величины имеют тогда вид
\[
L_{k} \equiv \alpha_{2 k+1}=i \sum_{l=2}^{2 k-1} T_{2 k, l} \Phi_{l},
\]

где $T_{2 k, l}$ – средняя строка матрицы $C_{2 k, l}$. Поскольку все элементы матриц, входящих в (7), содержат в качестве множителя $i$, вектор $\Phi_{m}$ не будет иметь мнимого множителя.

При реализации алгоритма на ЭВМ нужно учитывать необходимость проведения тождественных алгебраических преобразований над полиномами, коэффициенты которых являются обыкновенными дробями, и недопустимость приближенного представления дробей как в записи исходной информации, так и на всех промежуточных этапах. Поэтому для реализации алгоритма строится арифметика обыкновенных дробей, сохраняющая целочисленность числителя и знаменателя и обеспечивающая возможность выполнения в целых
числах операций умножения, сложения и сокращения обыкновенных дробей без использования арифметических операций, реализованных в системе команд ЭВМ, а также перевод числителя из двоичной систем́ы в десятичную [31].

Как видно из (7), для построения вектора Ф $_{m}$ необходимо хранить в оперативной памяти машины предыдущие векторы, что предъявляет определенные требования к объему памяти машины.

Разработанный алгоритм был применен С. Д. Щуко для вычисления ляпуновских величин некоторых систем вида (3).
a) Для системы
\[
\dot{x}=-y+P_{2}(x, y), \quad \dot{y}=x+Q_{2}(x, y),
\]

приведенной к виду (4), где
\[
Z(\boldsymbol{x}, \bar{z})=\alpha \bar{z}^{2}+\beta z \bar{z}+\gamma z^{2},
\]

получены три последовательные ляпуновские величины
\[
\begin{array}{c}
\alpha_{3}=\bar{\alpha} \bar{\beta}-\alpha \beta \\
\alpha_{5}=\frac{2}{3}\left(\bar{\alpha}^{2} \beta \bar{\gamma}-\alpha^{2} \bar{\beta} \gamma\right)+\frac{2}{3}\left(\bar{\beta}^{3} \gamma-\beta^{3} \bar{\gamma}\right)+\bar{\alpha} \beta^{2} \bar{\gamma}-\alpha \bar{\beta}^{2} \gamma \\
\alpha_{7}=8\left(\alpha^{2} \beta \bar{\beta}^{2} \gamma-\bar{\alpha}^{2} \bar{\beta} \beta^{2} \bar{\gamma}\right)+\frac{5}{4}\left(\bar{\alpha} \beta^{2} \bar{\gamma}^{2} \gamma-\alpha \bar{\beta}^{2} \gamma^{2} \bar{\gamma}\right)+ \\
+\frac{5}{8}\left(\bar{\beta}^{3} \gamma^{2} \bar{\gamma}-\beta^{3} \bar{\gamma}^{2} \gamma\right)+\frac{53}{4}\left(\alpha \beta \bar{\beta}^{3} \gamma-\bar{\alpha} \bar{\beta} \beta^{3} \bar{\gamma}\right)+\frac{69}{8}\left(\bar{\beta} \beta^{4} \bar{\gamma}-\beta \bar{\beta}^{4} \gamma\right) .
\end{array}
\]
б) Для системы
\[
\dot{x}=-y+P_{3}(x, y), \quad \dot{y}=x+Q_{3}(x, y),
\]

приведенной к виду (4), где
\[
Z(z, \bar{z})=\alpha \tilde{z}^{3}+\beta \bar{z}^{2} z+\gamma z^{2} \bar{z}+\delta z^{3},
\]

получено шесть последовательных ляпуновских величин
\[
\begin{aligned}
\alpha_{3}= & k_{3}(\beta-\bar{\beta}), \quad \alpha_{5}=k_{5}(\alpha \gamma-\bar{\alpha} \bar{\gamma}), \\
\alpha_{7}= & k_{7}\left[\left(\alpha^{2} \delta-\bar{\alpha}^{2} \bar{\delta}\right)+\frac{8}{3}(\alpha \bar{\gamma} \delta-\bar{\alpha} \gamma \bar{\delta})+\left(\gamma^{2} \delta-\bar{\gamma}^{2} \bar{\delta}\right)\right], \\
\alpha_{9}= & k_{9}\left[(\alpha \bar{\gamma} \delta-\bar{\alpha} \gamma \bar{\delta})+\frac{1}{3}\left(\gamma^{2} \bar{\delta}-\bar{\gamma}^{2} \delta\right)\right] \bar{\beta}, \\
\alpha_{11}= & k_{11}\left[(\alpha \bar{\gamma} \delta-\bar{\alpha} \gamma \bar{\delta})\left(\frac{207}{4} \tilde{\alpha} \bar{\gamma}+\delta \bar{\delta}\right)+\right. \\
& \left.+\left(\gamma^{2} \bar{\delta}-\bar{\gamma}^{2} \delta\right)\left(\frac{1}{3} \delta \bar{\delta}-134 \bar{\alpha} \bar{\gamma}+\frac{605}{12} \gamma \bar{\gamma}\right)\right], \\
\alpha_{13} \equiv & 0, \quad k_{j}=\text { const }>0,
\end{aligned}
\]

В случаях а) и б) обращение в нуль найденных ляпуновских величин дает необходимые и достаточные условия центра для систем (8) и (9), полученные впервые соответственно Каптейном [15] и К. Е. Малкиным [18].

У систем (8) и (9) при изменении коэффициентов (в том числе и линейных членов) из состояния равновесия $(0,0)$ не может появиться более трех (для системы (8)) или более пяти (для системы (9)) предельных циклов [10, 22].
в) Для системы
\[
\dot{x}=-y+P_{5}(x, y), \quad \dot{y}=x+Q_{5}(x, y),
\]

приведенной к виду (4), где
\[
Z(z, \bar{z})=\alpha \bar{z}^{5}+\beta \bar{z}^{4} z+\gamma \bar{z}^{3} z^{2}+\delta \bar{z}^{2} z^{3}+\varepsilon \bar{z} z^{4}+\zeta z^{5},
\]

получены первые три ляпуновские величины
\[
\begin{array}{l}
\alpha_{3}=k_{3}(\gamma-\bar{\gamma}), \quad \alpha_{5}=k_{5}[(\alpha \varepsilon-\bar{\alpha} \bar{\varepsilon})+(\beta \delta+\bar{\beta} \bar{\delta})], \\
\alpha_{7}=k_{7}\left[(\alpha \beta \zeta-\bar{\alpha} \bar{\beta} \bar{\zeta})+\frac{9}{2}\left(\alpha \delta^{2}-\bar{\alpha} \bar{\delta}^{2}\right)+\frac{3}{2}(\bar{\alpha} \beta \bar{\delta}-\alpha \bar{\beta} \delta)+\right. \\
+\frac{9}{2}(\alpha \bar{\delta} \zeta-\bar{\alpha} \delta \bar{\zeta})+\frac{9}{2}(\beta \bar{\delta} \varepsilon-\bar{\beta} \delta \bar{\varepsilon})+3(\beta \bar{\varepsilon} \zeta-\bar{\beta} \varepsilon \bar{\zeta})+ \\
\left.+\frac{3}{2}\left(\delta^{2} \bar{\varepsilon}-\bar{\delta}^{2} \varepsilon\right)+\frac{5}{2}(\delta \varepsilon \bar{\zeta}-\bar{\delta} \bar{\varepsilon} \xi)\right], \quad k_{j}=\mathrm{const}>0 . \\
\end{array}
\]
г) В общем случае системы
\[
\dot{x}=-y+P(x, y), \quad \dot{y}=x+Q(x, y),
\]

где $P(x, y)$ и $Q(x, y)$ – функции, аналитические в окрестности начала координат, разложения которых начинаются с членов порядка не меньше двух, для вычисления $\alpha_{3}$ в разложениях правых частей нужно удерживать члены до третьего порядка включительно, а для вычисления $\alpha_{5}$ – до пятого включительно. После приведения (10) к виду (4) получаем
\[
\begin{aligned}
Z(z, \bar{z})= & \alpha \bar{z}^{2}+\beta \bar{z} z+\gamma z^{2}+\delta \bar{z}^{3}+\varepsilon \bar{z}^{2} z+\varphi \bar{z} z^{2}+\varphi z^{3}+ \\
& +x \bar{z}^{4}+\lambda \bar{z}^{3} z+\mu \bar{z}^{2} z^{2}+
u \bar{z} z^{3}+\pi z^{4}+\rho \bar{z}^{5}+\sigma \bar{z}^{4} z+ \\
& +\tau \bar{z}^{3} z^{2}+\theta \bar{z}^{2} z^{3}+\zeta \bar{z} z^{4}+\chi z^{5} .
\end{aligned}
\]

Первые две ляпуновские величины будут
\[
\begin{array}{l}
\alpha_{3}=-\alpha \beta+\bar{\alpha} \bar{\beta}+\varepsilon-\bar{\varepsilon}, \\
\alpha_{5}=-\frac{2}{3}\left(\alpha^{2} \bar{\beta} \gamma-\bar{\alpha}^{2} \beta \bar{\gamma}\right)+\left(\bar{a} \beta^{2} \bar{\gamma}-\alpha \bar{\beta}^{2} \gamma\right)+\frac{2}{3}\left(\bar{\beta}^{3} \gamma-\beta^{3} \gamma\right)+ \\
+(\alpha \bar{\beta} \bar{\delta}-\bar{\alpha} \beta \delta)+\frac{2}{3}(\alpha \gamma \delta-\bar{a} \bar{\gamma} \bar{\delta})+2\left(\beta^{2} \delta-\bar{\beta}^{2} \bar{\delta}\right)+ \\
+\frac{5}{3}(\bar{\beta} \gamma \delta-\beta \bar{\gamma} \bar{\delta})+(\bar{\varepsilon}-\varepsilon) \bar{\alpha} \bar{\beta}+\frac{4}{3}(\bar{\varepsilon}-\varepsilon) \gamma \bar{\gamma}+ \\
+(\alpha \beta \varphi-\bar{\alpha} \bar{\beta} \bar{\varphi})+\frac{2}{3}(\alpha \gamma \bar{\varphi}-\bar{\alpha} \bar{\gamma} \varphi)+(\beta \bar{\gamma} \varphi-\bar{\beta} \gamma \bar{\varphi})+ \\
+\frac{2}{3}(\alpha \bar{\gamma} \psi-\bar{\alpha} \gamma \bar{\psi})+\frac{1}{3}(\beta \gamma \bar{\psi}-\bar{\beta} \bar{\gamma} \psi)+\left(\tilde{\delta}_{\bar{\varphi}}-\delta \varphi\right)+ \\
+(\varepsilon-\bar{\varepsilon}) \varepsilon-\frac{4}{3}(\gamma x-\bar{\gamma} \bar{x})+(\bar{\alpha} \lambda+\alpha \bar{\lambda})+ \\
+2(\bar{\beta} \bar{\lambda}-\beta \lambda)+(\beta \bar{\mu}-\bar{\beta} \mu)+\frac{1}{3}(\gamma \bar{v}-\bar{\gamma} v)+(\tau-\bar{\tau}) . \\
\end{array}
\]