Сингулярное разложение (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$ переводит сферу единичного радиуса в эллипсоид с полуосями, равными сингулярным числам и направленными вдоль левых сингулярных векторов. Правые сингулярные вектора~--- это, соответственно, прообразы полуосей эллипсоида.
Последовательное действие этих матриц на единичную окружность и сингулярные векторы показано на рис.~\ref{SVD_geom}.
В данном разложении не используются матрицы отражения, что приводит к возможности сингулярным числам быть отрицательными. Однако, из общих соображений понятно, что матрица эволюции непрерывного потока не может содержать отражений. Для матрицы \eqref{sol_linsys} это можно показать непосредственно.
где $\Vert\cdot\Vert$ обозначает евклидову норму. То есть, {\em длина произвольного вектора $\mathbf{x}$ изменяется матрицей $M$ минимум в $\sigma_2$ раз и максимум в $\sigma_1$ раз}.
\section{Связь собственных и сингулярных чисел}
Устойчивость или неустойчивость движения определяется не сингулярными числами, а собственными. Для двухмерной матрицы \eqref{SVD2x2} собственные числа можно выразить через сингулярные. Запишем характеристическое уравнение
Знак детерминанта определяет, являются ли собственные числа комплексными (комплексно-сопряжённая пара) или вещественными. Для любых сингулярных чисел, если сумма углов $\phi_1$ и $\phi_2$ равна $\pi/2+n\pi$, $n$~--- произвольное целое, детерминант будет отрицательным. И наоборот, если сумма углов кратна $\pi$, детерминант неотрицателен для любых сингулярных чисел.
В случае, если действие матрицы $M$ сохраняет площадь ($\operatorname{Det} M=1$), существуют инвариантные кривые второго порядка, которые под действием матрицы переходят сами в себя. В этой главе мы будем искать инвариантный эллипс.
Общее уравнение эллипса с центром в начале координат имеет вид:
где $a$ и $b$~--- полуоси эллипса, а $\theta$~--- угол между осью абсцисс и полуосью $a$. Так как матрица $M$~--- линейный оператор, достаточно найти один инвариантный эллипс, остальные можно получить путём масштабирования. Будем искать эллипс с площадью, равной площади единичного круга и большой полуосью $k\ge1$. В этом случае параметры в уравнении \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'$ с новыми параметрами
По условию инвариантности $k'=k$. Тогда, из выражения для $C'$ следует, что $\sin2\theta'=\sin2\theta$. Это уравнение имеет два решения, $\theta'=\theta$ и $\theta'=\pi/2-\theta$. Первое решение не имеет смысла, так как эллипс переходит сам в себя, что возможно только при $s=1$. Таким образом, остаётся только второе решение. Изменение угла ориентации не зависит от $k$, поэтому мы теперь можем найти наклон инвариантного эллипса $\phi_0$
Можно показать, что при возведении в квадрат условие \eqref{thetacrit} эквивалентно условию $D<0$, где $D$ определяется формулой \eqref{det_ch_eq}. Таким образом, параметры инвариантного эллипса определяются через сингулярное разложение как
Однако, мы нашли только два параметра эллипса, в то время как в сингулярном разложении их три. Существует ещё один параметр, а именно угол поворота по инвариантному эллипсу, $\psi$. Чтобы его найти, преобразуем систему отсчёта так, чтобы инвариантный эллипс превратился в круг. Для этого необходимо последовательно провести поворот на угол $-\theta$, а затем растяжение-сжатие осей на коэффициент $k$. Такое преобразование координат осуществляется матрицей
%
\begin{equation}
Q=S^{-1}(k)R(-\theta),\quad
S(k)=
\begin{pmatrix}
k & 0\\
0 & 1/k
\end{pmatrix},\quad
R(\theta)=
\begin{pmatrix}
\cos\theta&-\sin\theta\\
\sin\theta&\cos\theta
\end{pmatrix}
\end{equation}
%
В новой системе координат матрица $M$ должна быть просто матрицей поворота на угол $\psi$
Условием существования решения уравнений \eqref{psidef} является, естественно, $D\le0$. Угол поворота $\psi$ лежит в той же четверти, что и $\phi_1+\phi_2$, и в явном виде определяется как
где $\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)$
Решение системы \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\to0}\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} может быть переписано как
Величину $\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} имеет вид
Собственные числа матрицы эволюции: $e^{\lambda_1t}$ и $e^{\lambda_2t}$.
Поток \eqref{linsys} в общем случае не сохраняет площадь. Отношение площади в момент времени $t$ к площади в начальный момент времени определяется детерминантом матрицы эволюции $\operatorname{Det}{G(t)}$
Таким образом, поток \eqref{linsys} сохраняет площадь, тогда и только тогда, когда след матрицы скоростей равен нулю. Это равносильно условию равенства нулю суммы собственных чисел матрицы скоростей.
Точное выражение для максимального сингулярного числа $\sigma_1$ матрицы эволюции \eqref{sol_linsys} (и для показателя Ляпунова на конечное время, соответственно) весьма громоздко, поэтому приведём только его
где $F(t)$~--- ограниченная периодическая функция с периодом $w/2$. Воспользовавшись определением \eqref{lyap_def_sigma}, нетрудно определить показатель Ляпунова системы \eqref{linsys}.
