\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{color}
\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}

\topmargin=-1.8cm
\oddsidemargin=-15mm
\evensidemargin=-15mm
\textheight=24.5cm
\textwidth=19cm
\tolerance=1000

\parskip=5pt plus 4pt minus 2pt
\tolerance=9000
%
\begin{document}
\section{Сингулярное разложение матрицы}
Сингулярное разложение (singular value decomposition)~--- представление произвольной матрицы $m\times n$ в виде
%
\begin{equation}
M=U\Sigma V,
\label{SVD}
\end{equation}
%
где $U$ и $V$~--- унитарные матрицы $m\times m$ и $n\times n$, соответственно, $\Sigma$~--- диагональная матрица $m\times n$.
Числа, стоящие на диагонали $\Sigma$, называются сингулярными числами матрицы $M$. Единичные вектора $u$ и $v$ такие, что $Mv=\sigma u$ и $M^*u=\sigma v$,
$\sigma$~--- сингулярное число $M$, называются, соответственно, левыми и правыми сингулярными векторами матрицы $M$.
Если матрица $M$ является вещественной, то её сингулярные числа также вещественны, а $U$ и $V$~--- ортогональные матрицы. Так как матрица $\Sigma$ и, соответственно, сингулярное разложение определены с точностью до перестановки сингулярных чисел, то можно потребовать упорядочения сингулярных чисел на диагонали $\Sigma$ по невозрастанию. Такое разложение является единственным.

Если матрица $M$ является квадратной, то сингулярное разложение имеет очень простой геометрический смысл. Действие любой матрицы на вектор можно представить в виде трёх последовательных операций: поворот/отражение матрицей $V$, растяжение/сжатие по координатным осям матрицей $\Sigma$, второй поворот/отражение матрицей $U$. Таким образом, матрица $M$ переводит сферу единичного радиуса в эллипсоид с полуосями, равными сингулярным числам и направленными вдоль левых сингулярных векторов. Правые сингулярные вектора~--- это, соответственно, прообразы полуосей эллипсоида.

Рассмотрим сингулярное разложение вещёственной матрицы $2\times 2$:
%
\begin{equation}
M=U\Sigma V \Rightarrow
\begin{pmatrix}
a&b\\c&d
\end{pmatrix}
=
\begin{pmatrix}
\cos\phi_2&-\sin\phi_2\\
\sin\phi_2&\cos\phi_2
\end{pmatrix}
\begin{pmatrix}
\sigma_1&0\\
0&\sigma_2
\end{pmatrix}
\begin{pmatrix}
\cos\phi_1&-\sin\phi_1\\
\sin\phi_1&\cos\phi_1
\end{pmatrix}.
\label{SVD2x2}
\end{equation}
%
Последовательное действие этих матриц на единичную окружность и сингулярные векторы показано на рис.~\ref{SVD_geom}.
В данном разложении не используются матрицы отражения, что приводит к возможности сингулярным числам быть отрицательными. Однако, из общих соображений понятно, что матрица эволюции непрерывного потока не может содержать отражений. Для матрицы \eqref{sol_linsys} это можно показать непосредственно.
%
\begin{figure}[!htb]
\centerline{\includegraphics[width=0.8\textwidth]{SVD.eps}}
\caption{Геометрический смысл сингулярного разложения матрицы.}
\label{SVD_geom}
\end{figure}
%

Перемножая матрицы, получаем систему из четырёх уравнений с четырьмя неизвестными
%
\begin{equation}
\begin{pmatrix}
a&b\\c&d
\end{pmatrix}
=
\begin{pmatrix}
\sigma_1\cos\phi_1\cos\phi_2-\sigma_2\sin\phi_1\sin\phi_2&
-\sigma_1\sin\phi_1\cos\phi_2-\sigma_2\cos\phi_1\sin\phi_2\\
\sigma_1\cos\phi_1\sin\phi_2+\sigma_2\sin\phi_1\cos\phi_2&
-\sigma_1\sin\phi_1\sin\phi_2+\sigma_2\cos\phi_1\cos\phi_2
\end{pmatrix}.
\label{SVD2x2eq}
\end{equation}
%
Введём обозначения
%
\begin{equation}
\begin{aligned}
\alpha&=a+d, &\beta&=a-d, &\gamma&=c+b, &\delta&=c-b,\\
\xi&=\sigma_1+\sigma_2, &\eta&=\sigma_1-\sigma_2, &\Phi&=\phi_1+\phi_2, &\Psi&=\phi_2-\phi_1.
\end{aligned}
\label{SVD2x2def}
\end{equation}
%
Вычитая и складывая уравнения \eqref{SVD2x2eq} и используя введённые обозначения \eqref{SVD2x2def}, получаем
%
\begin{equation}
\alpha=\xi\cos\Phi,\quad \beta=\eta\cos\Psi,\quad \gamma=\eta\sin\Psi,\quad \delta=\xi\sin\Phi.
\label{SVD2x2fineq}
\end{equation}
%
Решение системы \eqref{SVD2x2fineq} имеет вид
%
\begin{equation}
\xi=\sqrt{\alpha^2+\delta^2},\quad \eta=\sqrt{\beta^2+\gamma^2},\quad
\Phi=\operatorname{arctg2}{(\delta,\alpha)}, \quad \Psi=\operatorname{arctg2}{(\gamma,\beta)}.
\label{SVD2x2sol}
\end{equation}
%
Здесь $\operatorname{arctg2}{(y,x)}$~--- это угол между вектором $(x,y)$ и осью $x$. Он может быть определён как
%
\begin{equation}
\operatorname{arctg2}{(y,x)}=\left\{
\begin{aligned}
&\arctg{(y/x)}, &x\ge 0,\\
&\arctg{(y/x)}+\pi, &x<0.
\end{aligned}\right.
\label{arctg2}
\end{equation}
%
Окончательный вид решения
%
\begin{equation}
\begin{gathered}
\sigma_1=\frac{\sqrt{(a+d)^2+(c-b)^2}+\sqrt{(a-d)^2+(b+c)^2}}{2},\quad
\sigma_2=\frac{\sqrt{(a+d)^2+(c-b)^2}-\sqrt{(a-d)^2+(b+c)^2}}{2},\\
\phi_1=\frac{\operatorname{arctg2}{(c-b,\,a+d)}-\operatorname{arctg2}{(c+b,\,a-d)}}{2},\quad
\phi_2=\frac{\operatorname{arctg2}{(c-b,\,a+d)}+\operatorname{arctg2}{(c+b,\,a-d})}{2}.
\end{gathered}
\label{SVD2x2finsol}
\end{equation}
%
Из решения очевидно, что сингулярные числа упорядочены по невозрастанию, то есть $\sigma_1\ge\sigma_2$.
Произведение $\sigma_1\sigma_2$ определяет отношение конечной и начальной площадей и равно, естественно, $\operatorname{Det}M$.
Из определения сингулярного разложения следует, что
%
\begin{equation}
\sigma_1\geqslant\frac{\Vert M\mathbf{x}\Vert}{\Vert\mathbf{x}\Vert}\geqslant\sigma_2,
\label{SVD_mainprop}
\end{equation}
%
где $\Vert\cdot\Vert$ обозначает евклидову норму. То есть, {\em длина произвольного вектора $\mathbf{x}$ изменяется матрицей $M$ минимум в $\sigma_2$ раз и максимум в $\sigma_1$ раз}.

\section{Связь собственных и сингулярных чисел}
Устойчивость или неустойчивость движения определяется не сингулярными числами, а собственными. Для двухмерной матрицы \eqref{SVD2x2} собственные числа можно выразить через сингулярные. Запишем характеристическое уравнение
%
\begin{equation}
\operatorname{Det}{(M-\lambda)}=0 \Longrightarrow \lambda^2-(a+d)\lambda+ad-bc=0.
\label{ch_eq}
\end{equation}
%
Подставив значения $a$, $b$, $c$ и $d$ из \eqref{SVD2x2eq}
%
\begin{displaymath}
\begin{gathered}
a+d=\sigma_1(\cos\phi_1\cos\phi_2-\sin\phi_1\sin\phi_2)+\sigma_2(\cos\phi_1\cos\phi_2-\sin\phi_1\sin\phi_2)=(\sigma_1+\sigma_2)\cos{(\phi_1+\phi_2)},\\
ad-bc=\operatorname{Det}M=\sigma_1\sigma_2,
\end{gathered}
\end{displaymath}
%
получаем
%
\begin{equation}
\lambda^2-(\sigma_1+\sigma_2)\cos{(\phi_1+\phi_2)}\lambda+\sigma_1\sigma_2=0.
\label{ch_eq_sigma}
\end{equation}
%
Детерминант характеристического уравнения:
%
\begin{equation}
D=(\sigma_1+\sigma_2)^2\cos^2{(\phi_1+\phi_2)}-4\sigma_1\sigma_2.
\label{det_ch_eq}
\end{equation}
%
Знак детерминанта определяет, являются ли собственные числа комплексными (комплексно-сопряжённая пара) или вещественными. Для любых сингулярных чисел, если сумма углов $\phi_1$ и $\phi_2$ равна $\pi/2+n\pi$, $n$~--- произвольное целое, детерминант будет отрицательным. И наоборот, если сумма углов кратна $\pi$, детерминант неотрицателен для любых сингулярных чисел.
Сами собственные числа равны
%
\begin{equation}
\lambda_{1,2}=\frac{(\sigma_1+\sigma_2)\cos{(\phi_1+\phi_2)}\pm\sqrt{D}}{2}.
\label{lambda_ch_eq}
\end{equation}
%
Произведение собственных чисел равно произведению сингулярных и равно детерминанту матрицы.

\section{Инвариантный эллипс}
В случае, если действие матрицы $M$ сохраняет площадь ($\operatorname{Det} M=1$), существуют инвариантные кривые второго порядка, которые под действием матрицы переходят сами в себя. В этой главе мы будем искать инвариантный эллипс.

Общее уравнение эллипса с центром в начале координат имеет вид:
%
\begin{equation}
\begin{gathered}
Ax^2+By^2+Cxy=1,\\
A=\frac{\cos^2\theta}{a^2}+\frac{\sin^2\theta}{b^2},\qquad
B=\frac{\sin^2\theta}{a^2}+\frac{\cos^2\theta}{b^2},\\
C=\frac{b^2-a^2}{a^2b^2}\sin2\theta,
\end{gathered}
\label{ellbase}
\end{equation}
%
где $a$ и $b$~--- полуоси эллипса, а $\theta$~--- угол между осью абсцисс и полуосью $a$. Так как матрица $M$~--- линейный оператор, достаточно найти один инвариантный эллипс, остальные можно получить путём масштабирования. Будем искать эллипс с площадью, равной площади единичного круга и большой полуосью $k\ge 1$. В этом случае параметры в уравнении \eqref{ellbase} имеют вид
%
\begin{equation}
\begin{gathered}
A=\frac{\cos^2\theta}{k^2}+k^2\sin^2\theta,\qquad
B=\frac{\sin^2\theta}{k^2}+k^2\cos^2\theta,\\
C=-\frac{k^4-1}{k^2}\sin2\theta.
\end{gathered}
\label{ellpars0}
\end{equation}
%
Угол $\theta$ определяет направление большой полуоси. Сумма коэффициентов $A+B=(k^4+1)/k^2$ не зависит от ориентации эллипса и определяется только его формой.

Матрица $M=U\Sigma V$ последовательно поворачивает ($V$), растягивает-сжимает ($\Sigma$) и снова поворачивает ($U$) эллипс. При этом большая полуось может изменится только под действием растяжения-сжатия. Рассмотрим действие матрицы $\Sigma$. Так как, $\operatorname{Det}\Sigma=1$,
%
\begin{displaymath}
\Sigma=
\begin{pmatrix}
\sigma & 0\\
0& 1/\sigma
\end{pmatrix},\quad \sigma\ge 1.
\end{displaymath}
%
Действие $\Sigma$ преобразует координаты $x$ и $y$
%
\begin{gather}
x'=\sigma x,\qquad y'=y/\sigma,\\
x=x'/\sigma,\qquad y=\sigma y'.
\end{gather}
%
Подставив $x$ и $y$ в уравнение \eqref{ellbase} с учётом \eqref{ellpars0} получаем уравнение эллипса в новых координатах $x'$ и $y'$ с новыми параметрами
%
\begin{equation}
\begin{gathered}
A'=\frac{\cos^2\theta'}{k'^2}+k'^2\sin^2\theta'=\frac{\cos^2\theta}{\sigma^2k^2}+\frac{k^2}{\sigma^2}\sin^2\theta,\qquad
B'=\frac{\sin^2\theta'}{k'^2}+k'^2\cos^2\theta'=\frac{\sigma^2}{k^2}\sin^2\theta+\sigma^2k^2\cos^2\theta,\\
C'=-\frac{k'^4-1}{k'^2}\sin2\theta'=-\frac{k^4-1}{k^2}\sin2\theta.
\end{gathered}
\label{ellpars1}
\end{equation}
%
По условию инвариантности $k'=k$. Тогда, из выражения для $C'$ следует, что $\sin2\theta'=\sin2\theta$. Это уравнение имеет два решения, $\theta'=\theta$ и $\theta'=\pi/2-\theta$. Первое решение не имеет смысла, так как эллипс переходит сам в себя, что возможно только при $s=1$. Таким образом, остаётся только второе решение. Изменение угла ориентации не зависит от $k$, поэтому мы теперь можем найти наклон инвариантного эллипса $\phi_0$
%
\begin{equation}
\phi_0\xrightarrow{V}\theta=\phi_0+\phi_1\xrightarrow{\Sigma}\theta'=\frac{\pi}{2}-(\phi_0+\phi_1)\xrightarrow{U}\frac{\pi}{2}-(\phi_0+\phi_1)+\phi_2=\left\{
\begin{aligned}
&\phi_0\\
&\phi_0+\pi
\end{aligned}
\right.,
\label{phi0transform}
\end{equation}
%
откуда получаем
%
\begin{equation}
\phi_0=\frac{\phi_2-\phi_1}{2}\pm\frac{\pi}{4}.
\label{phi0}
\end{equation}
%
К выбору знака мы вернёмся ниже.

Вернёмся к поиску $k$. Так как $A+B$ является инвариантом относительно поворота, $k$ найдём из уравнения $A+B=A'+B'$
%
\begin{equation}
\frac{k^4+1}{k^2}=\frac{\sigma^4k^4+1}{\sigma^2k^2}\cos^2{\theta}+\frac{k^4+\sigma^4}{\sigma^2k^2}\sin^2\theta.
\label{inveq}
\end{equation}
%
После тривиальных преобразований получаем
%
\begin{equation}
\begin{gathered}
\sigma^2k^4+\sigma^2=(\sigma^4k^4+1)\cos^2\theta+(k^4+\sigma^4)(1-\cos^2\theta),\\
k^4(\sigma^2-1-(\sigma^4-1)\cos^2\theta)=\sigma^2(\sigma^2-1)-(\sigma^4-1)\cos^2\theta,\\
k^4(1-(\sigma^2+1)\cos^2\theta)=\sigma^2-(\sigma^2+1)\cos^2\theta.
\end{gathered}
\label{inveq1}
\end{equation}
%
%
\begin{equation}
k^4=\frac{\sigma^2-(\sigma^2+1)\cos^2\theta}{1-(\sigma^2+1)\cos^2\theta}=1+\frac{\sigma^2-1}{1-(\sigma^2+1)\cos^2\theta}=1+\frac{2(\sigma^2-1)}{1-\sigma^2-(\sigma^2+1)\cos2\theta}.
\label{inveq2}
\end{equation}
%
Так как $k\ge 1$, $1-\sigma^2-(\sigma^2+1)\cos2\theta> 0$, или
%
\begin{equation}
-\cos2\theta>\frac{\sigma^2-1}{\sigma^2+1},
\label{thetacrit}
\end{equation}
%
из чего следует, что $\cos2\theta<0$. Это условие позволяет определить знак в формуле \eqref{phi0}
%
\begin{equation}
\cos2\theta=\cos2(\phi_0+\phi_1)=\cos2\left(\frac{\phi_2-\phi_1}{2}\pm\frac{\pi}{4}+\phi_1\right)=\cos\left(\frac{\pi}{2}\pm(\phi_1+\phi_2)\right)=\mp\sin(\phi_1+\phi_2),
\end{equation}
%
Можно показать, что при возведении в квадрат условие \eqref{thetacrit} эквивалентно условию $D<0$, где $D$ определяется формулой \eqref{det_ch_eq}. Таким образом, параметры инвариантного эллипса определяются через сингулярное разложение как
%
\begin{equation}
\phi_0=\frac{\phi_2-\phi_1}{2}+\xi\frac{\pi}{4},\quad k=\sqrt[4]{1+\frac{2(\sigma^2-1)}{1-\sigma^2+\xi(\sigma^2+1)\sin(\phi_1+\phi_2)}},\quad \xi=\operatorname{sgn}\sin(\phi_1+\phi_2).
\end{equation}
%

\section{Показатель Ляпунова}
Пусть дана система нелинейных уравнений
%
\begin{equation}
\mathbf{\dot x}=\mathbf{f}(\mathbf{x},t),\quad
\mathbf{x}=(x_1,\dotsc,x_n),\quad
\mathbf{f}(\mathbf{x},t)=(f_1(x_1,\dotsc,x_n,t),\dots,f_n(x_1,\dotsc,x_n,t)).
\label{nonlinsys}
\end{equation}
%
Показатель Ляпунова в точке $\mathbf{x_0}$ определяется как
%
\begin{equation}
\Lambda(\mathbf{x_0})=\lim_{t\to\infty}\lim_{\Vert\delta \mathbf{x}(0)\Vert\to 0}\frac{\ln(\Vert\delta \mathbf{x}(t)\Vert/\Vert\delta \mathbf{x}(0)\Vert)}{t},
\label{lyap_def}
\end{equation}
%
где $\delta \mathbf{x}(t)=\mathbf{x_1}(t)-\mathbf{x_0}(t)$, $\mathbf{x_0}(t)$ и $\mathbf{x_1}(t)$~--- решения системы \eqref{nonlinsys}, $\mathbf{x_0}(0)=\mathbf{x_0}$. Данный предел существует и одинаков почти для любого выбора $\delta \mathbf{x}(0)$. Геометрический смысл показателя Ляпунова состоит в том, что две близлежащих траектории расходятся в пространстве в среднем экспоненциально по времени с показателем экспоненты, равным показателю Ляпунова.

В силу малости $\Vert\delta \mathbf{x}(0)\Vert$ и, соответственно, $\Vert\delta \mathbf{x}(t)\Vert$ можно линеаризовать
систему \eqref{nonlinsys} в окрестности траектории $\mathbf{x_0}(t)$, получив систему нестационарных линейных уравнений
%
\begin{equation}
\begin{pmatrix}
\delta\dot x_1\\
\hdotsfor{1}\\
\delta\dot x_n
\end{pmatrix}=J(t)
\begin{pmatrix}
\delta x_1\\
\hdotsfor{1}\\
\delta x_n
\end{pmatrix},
\label{linearization}
\end{equation}
%
где $J(t)$~--- якобиан системы \eqref{nonlinsys} вдоль траектории $\mathbf{x_0}(t)$
%
\begin{equation}
J(t)=
\begin{pmatrix}
\dfrac{\partial f_1(\mathbf{x_0}(t),t)}{\partial x_1}&\dots&\dfrac{\partial f_1(\mathbf{x_0}(t),t)}{\partial x_n}\\
\hdotsfor{3}\\
\dfrac{\partial f_n(\mathbf{x_0}(t),t)}{\partial x_1}&\dots&\dfrac{\partial f_n(\mathbf{x_0}(t),t)}{\partial x_n}
\end{pmatrix}.
\label{jacobian}
\end{equation}
%
Решение системы \eqref{linearization} можно представить в виде
%
\begin{equation}
\begin{pmatrix}
\delta x_1(t)\\
\hdotsfor{1}\\
\delta x_n(t)
\end{pmatrix}=G(t)
\begin{pmatrix}
\delta x_1(0)\\
\hdotsfor{1}\\
\delta x_n(0)
\end{pmatrix},
\label{evol_mat}
\end{equation}
%
где $G(t)$~--- матрица эволюции. Матрица эволюции подчиняется дифференциальному уравнению, которое получается подстановкой \eqref{evol_mat} в \eqref{linearization}
%
\begin{equation}
\dot G=JG,
\label{evol_mat_diffur}
\end{equation}
%
с начальным условием $G(0)=I$, где $I$~--- единичная матрица.

Максимальное значение $\lim\limits_{\Vert\delta \mathbf{x}(0)\Vert\to 0}\dfrac{\Vert\delta \mathbf{x}(t)\Vert}{\Vert\delta \mathbf{x}(0)\Vert}$ для системы \eqref{linearization} будет равно $\sigma_1(G(t))$~--- максимальному сингулярному числу матрицы $G(t)$. Если $\lim\limits_{t\to\infty}\dfrac{\sigma_2(G(t))}{\sigma_1(G(t))}=0$, где $\sigma_2(G(t))$~--- следующее по величине сингулярное число матрицы $G(t)$, то
определение \eqref{lyap_def} может быть переписано как
%
\begin{equation}
\Lambda=\lim_{t\to\infty}\frac{\ln\sigma_1(G(t))}{t}=\lim_{t\to\infty}\Lambda(t).
\label{lyap_def_sigma}
\end{equation}
%
Величину $\Lambda(t)$ будем называть {\em показателем Ляпунова за время $t$}. Это есть {\em отношение логарифма максимально возможного растяжения вектора ко времени}. Также определим {\em мгновенный показатель Ляпунова} $\Lambda_0$ как показатель Ляпунова системы стационарных линейных уравнений.
%
\begin{equation}
\begin{pmatrix}
\delta\dot x_1\\
\hdotsfor{1}\\
\delta\dot x_n
\end{pmatrix}=J(0)
\begin{pmatrix}
\delta x_1\\
\hdotsfor{1}\\
\delta x_n
\end{pmatrix}.
\label{linearization_stat}
\end{equation}
%
Мгновенный показатель Ляпунова определяет скорость экспоненциального разбегания траекторий в данной точке пространства в данный момент времени.

\section{Решение системы стационарных линейных уравнений}
Пусть задана система линейных уравнений с матрицей скоростей $J$, не зависящей от времени
%
\begin{equation}
\begin{pmatrix}
\dot x\\
\dot y
\end{pmatrix}=J
\begin{pmatrix}
x\\y
\end{pmatrix},\quad
J=
\begin{pmatrix}
a&b\\
c&d
\end{pmatrix}.
\label{linsys}
\end{equation}
%
Тогда, матрица эволюции $G$ системы \eqref{linsys} имеет вид
%
\begin{equation}
G(t)=\left\{
\begin{aligned}
&A_1e^{\lambda_1t}+A_2e^{\lambda_2t}, &D>0,\\
&e^{\lambda t}\left(I+\frac{K}{2}t\right), &D=0,\\
&e^{\lambda t}\left(I\cos{(\omega t)}+\frac{K}{\sqrt{-D}}\sin{(\omega t)}\right), &D<0,
\end{aligned}\right.
\label{sol_linsys}
\end{equation}
%
где
%
\begin{equation}
\begin{gathered}
I=\begin{pmatrix}1&0\\0&1\end{pmatrix},\quad
K=\begin{pmatrix}a-d&2b\\2c&d-a\end{pmatrix},\quad
A_1=\frac{\sqrt{D}I+K}{2\sqrt{D}},\quad
A_2=\frac{\sqrt{D}I-K}{2\sqrt{D}},\\
D=4bc+(a-d)^2,\quad \lambda=\frac{a+d}{2},\quad \omega=\frac{\sqrt{-D}}{2},\quad
\quad \lambda_1=\frac{a+d+\sqrt{D}}{2},\quad \lambda_2=\frac{a+d-\sqrt{D}}{2}.
\end{gathered}
\label{Amatrix}
\end{equation}
%
Собственные числа матрицы эволюции: $e^{\lambda_1t}$ и $e^{\lambda_2t}$.
Поток \eqref{linsys} в общем случае не сохраняет площадь. Отношение площади в момент времени $t$ к площади в начальный момент времени определяется детерминантом матрицы эволюции $\operatorname{Det}{G(t)}$
%
\begin{equation}
\frac{S(t)}{S(0)}=\operatorname{Det}{G(t)}=e^{(a+d)t}=e^{\operatorname{Tr}Jt}.
\label{square}
\end{equation}
%
Таким образом, поток \eqref{linsys} сохраняет площадь, тогда и только тогда, когда след матрицы скоростей равен нулю. Это равносильно условию равенства нулю суммы собственных чисел матрицы скоростей.

Точное выражение для максимального сингулярного числа $\sigma_1$ матрицы эволюции \eqref{sol_linsys} (и для показателя Ляпунова на конечное время, соответственно) весьма громоздко, поэтому приведём только его
асимптотику при $t\to\infty$:
%
\begin{equation}
\sigma_1(t)=
\left\{
\begin{aligned}
&\sqrt{\frac{(b+c)^2+(a-d)^2}{4bc+(a-d)^2}}e^{\lambda_1t},&D>0,\\
&|b-c|te^{\lambda t},&D=0,\\
&F(t)e^{\lambda t},&D<0,
\end{aligned}
\right.
\label{sigma_of_evol}
\end{equation}
%
где $F(t)$~--- ограниченная периодическая функция с периодом $w/2$. Воспользовавшись определением \eqref{lyap_def_sigma}, нетрудно определить показатель Ляпунова системы \eqref{linsys}.
\end{document}
