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

где
$$\dot{x}=ax+by+P(x,y),\dot{y}=cx+dy+Q(x,y),$$
$$P(x,y)=\sum _{i+j⩾2}{Q}_{ij}{x}^i{y}^j,Q(x,y)=\sum _{i+j⩾2}{b}_{ij}{x}^i{y}^j$$

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

Следуя Ляпунову и Андронову $[17,2]$, строим в окрестности начала координат функцию
$$\psi ({\rho }_0,p)=(e^{-\frac{2\pi }{\omega }p}-1){\rho }_0+{\alpha }_2{\rho }_0^2+{\alpha }_3{\rho }_0^3+\dots .$$

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

Если $p=0$, то первый не равный нулю коэффициент ${\alpha }_i$ называется ляпуновской величиной: ${\alpha }_3\equiv L_1$ — первая ляпуновская величина; если ${\alpha }_3=0,{\alpha }_{5}eq0$, то ${\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 ({\rho }_0,p)$ и соответственно один или два предельных цикла на фазовой плоскости вокруг начала координат.

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

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

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

Рассмотрим поведение динамической системы в малой окрестности точки на границе области устойчивости, в которой первая ляпуновская величина обращается в нуль. Рис. $2a,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(L_1{L}_2<0).$$

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

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

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

где
$$\begin{array}{l}Z(z,\overline{z})=\sum _{m=2}^n\sum _{k=0}^m{a}_{km}{z}^k{\overline{z}}^{m-k},\\ \overline{Z}(z,\overline{z})=\sum _{m=2}^n\sum _{k=0}^m{\overline{a}}_{km}{\overline{z}}^k{z}^{m-k}.\end{array}$$

Наличие центра в начале координат равносильно существованию голоморфного интеграла
$$\mathrm{\Phi }(z,\overline{z})=\sum _{m=0}^{\mathrm{\infty }}\sum _{k=0}^m{\mathrm{\Phi }}_{k,m-k}{z}^k{\overline{z}}^{m-k}$$

уравнения
$$L(\mathrm{\Phi })\equiv \frac{\partial \mathrm{\Phi }}{\partial z}\cdot \frac{dz}{dt}+\frac{\partial \mathrm{\Phi }}{\partial \overline{z}}\cdot \frac{d\overline{z}}{dt}.$$

Функции $\mathrm{\Phi }(z,\overline{z})$ сопоставляется последовательность $(m+1)$-мерных комплексных векторов ${\mathrm{\Phi }}_m$ с компонентами, являющимися коэффициентами однородных полиномов степени $m$ из (5):
$${\mathrm{\Phi }}_m=\{{\mathrm{\Phi }}_{0,m};{\mathrm{\Phi }}_{1,m-1};\dots ;{\mathrm{\Phi }}_{m-1,1};{\mathrm{\Phi }}_{m,0}\}.$$

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

Векторы ${\mathrm{\Phi }}_m$ находятся по рекуррентной формуле
$${\mathrm{\Phi }}_m=-C_{m,m}^{-1}\left\{\sum _{l=2}^{m-1}{C}_{m,l}{\mathrm{\Phi }}_l\right\},m>2.$$

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

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

Ляпуновские величины имеют тогда вид
$$L_k\equiv {\alpha }_{2k+1}=i\sum _{l=2}^{2k-1}{T}_{2k,l}{\mathrm{\Phi }}_l,$$

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

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

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

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

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

получены три последовательные ляпуновские величины
$$\begin{array}{c}{\alpha }_3=\overline{\alpha }\overline{\beta }-\alpha \beta \\ {\alpha }_5=\frac{2}{3}({\overline{\alpha }}^2\beta \overline{\gamma }-{\alpha }^2\overline{\beta }\gamma )+\frac{2}{3}({\overline{\beta }}^3\gamma -{\beta }^3\overline{\gamma })+\overline{\alpha }{\beta }^2\overline{\gamma }-\alpha {\overline{\beta }}^2\gamma \\ {\alpha }_7=8({\alpha }^2\beta {\overline{\beta }}^2\gamma -{\overline{\alpha }}^2\overline{\beta }{\beta }^2\overline{\gamma })+\frac{5}{4}(\overline{\alpha }{\beta }^2{\overline{\gamma }}^2\gamma -\alpha {\overline{\beta }}^2{\gamma }^2\overline{\gamma })+\\ +\frac{5}{8}({\overline{\beta }}^3{\gamma }^2\overline{\gamma }-{\beta }^3{\overline{\gamma }}^2\gamma )+\frac{53}{4}(\alpha \beta {\overline{\beta }}^3\gamma -\overline{\alpha }\overline{\beta }{\beta }^3\overline{\gamma })+\frac{69}{8}(\overline{\beta }{\beta }^4\overline{\gamma }-\beta {\overline{\beta }}^4\gamma ).\end{array}$$
б) Для системы
$$\dot{x}=-y+P_3(x,y),\dot{y}=x+Q_3(x,y),$$

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

получено шесть последовательных ляпуновских величин
$$\begin{array}{rl}{\alpha }_3=& k_3(\beta -\overline{\beta }),{\alpha }_5=k_5(\alpha \gamma -\overline{\alpha }\overline{\gamma }),\\ {\alpha }_7=& k_7[({\alpha }^2\delta -{\overline{\alpha }}^2\overline{\delta })+\frac{8}{3}(\alpha \overline{\gamma }\delta -\overline{\alpha }\gamma \overline{\delta })+({\gamma }^2\delta -{\overline{\gamma }}^2\overline{\delta })],\\ {\alpha }_9=& k_9[(\alpha \overline{\gamma }\delta -\overline{\alpha }\gamma \overline{\delta })+\frac{1}{3}({\gamma }^2\overline{\delta }-{\overline{\gamma }}^2\delta )]\overline{\beta },\\ {\alpha }_{11}=& k_{11}[(\alpha \overline{\gamma }\delta -\overline{\alpha }\gamma \overline{\delta })(\frac{207}{4}\tilde{\alpha }\overline{\gamma }+\delta \overline{\delta })+\\ & +({\gamma }^2\overline{\delta }-{\overline{\gamma }}^2\delta )(\frac{1}{3}\delta \overline{\delta }-134\overline{\alpha }\overline{\gamma }+\frac{605}{12}\gamma \overline{\gamma })],\\ {\alpha }_{13}\equiv & 0,k_j=\text{ const }>0,\end{array}$$

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

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

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

получены первые три ляпуновские величины
$$\begin{array}{l}{\alpha }_3=k_3(\gamma -\overline{\gamma }),{\alpha }_5=k_5[(\alpha \epsilon -\overline{\alpha }\overline{\epsilon })+(\beta \delta +\overline{\beta }\overline{\delta })],\\ {\alpha }_7=k_7[(\alpha \beta \zeta -\overline{\alpha }\overline{\beta }\overline{\zeta })+\frac{9}{2}(\alpha {\delta }^2-\overline{\alpha }{\overline{\delta }}^2)+\frac{3}{2}(\overline{\alpha }\beta \overline{\delta }-\alpha \overline{\beta }\delta )+\\ +\frac{9}{2}(\alpha \overline{\delta }\zeta -\overline{\alpha }\delta \overline{\zeta })+\frac{9}{2}(\beta \overline{\delta }\epsilon -\overline{\beta }\delta \overline{\epsilon })+3(\beta \overline{\epsilon }\zeta -\overline{\beta }\epsilon \overline{\zeta })+\\ +\frac{3}{2}({\delta }^2\overline{\epsilon }-{\overline{\delta }}^2\epsilon )+\frac{5}{2}(\delta \epsilon \overline{\zeta }-\overline{\delta }\overline{\epsilon }\xi )],k_j=\mathrm{const}>0.\end{array}$$
г) В общем случае системы
$$\dot{x}=-y+P(x,y),\dot{y}=x+Q(x,y),$$

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

Первые две ляпуновские величины будут
$$\begin{array}{l}{\alpha }_3=-\alpha \beta +\overline{\alpha }\overline{\beta }+\epsilon -\overline{\epsilon },\\ {\alpha }_5=-\frac{2}{3}({\alpha }^2\overline{\beta }\gamma -{\overline{\alpha }}^2\beta \overline{\gamma })+(\overline{a}{\beta }^2\overline{\gamma }-\alpha {\overline{\beta }}^2\gamma )+\frac{2}{3}({\overline{\beta }}^3\gamma -{\beta }^3\gamma )+\\ +(\alpha \overline{\beta }\overline{\delta }-\overline{\alpha }\beta \delta )+\frac{2}{3}(\alpha \gamma \delta -\overline{a}\overline{\gamma }\overline{\delta })+2({\beta }^2\delta -{\overline{\beta }}^2\overline{\delta })+\\ +\frac{5}{3}(\overline{\beta }\gamma \delta -\beta \overline{\gamma }\overline{\delta })+(\overline{\epsilon }-\epsilon )\overline{\alpha }\overline{\beta }+\frac{4}{3}(\overline{\epsilon }-\epsilon )\gamma \overline{\gamma }+\\ +(\alpha \beta \phi -\overline{\alpha }\overline{\beta }\overline{\phi })+\frac{2}{3}(\alpha \gamma \overline{\phi }-\overline{\alpha }\overline{\gamma }\phi )+(\beta \overline{\gamma }\phi -\overline{\beta }\gamma \overline{\phi })+\\ +\frac{2}{3}(\alpha \overline{\gamma }\psi -\overline{\alpha }\gamma \overline{\psi })+\frac{1}{3}(\beta \gamma \overline{\psi }-\overline{\beta }\overline{\gamma }\psi )+({\tilde{\delta }}_{\overline{\phi }}-\delta \phi )+\\ +(\epsilon -\overline{\epsilon })\epsilon -\frac{4}{3}(\gamma x-\overline{\gamma }\overline{x})+(\overline{\alpha }\lambda +\alpha \overline{\lambda })+\\ +2(\overline{\beta }\overline{\lambda }-\beta \lambda )+(\beta \overline{\mu }-\overline{\beta }\mu )+\frac{1}{3}(\gamma \overline{v}-\overline{\gamma }v)+(\tau -\overline{\tau }).\end{array}$$