news.shamcode.ru | Математика равновесия: как уравнение Ляпунова держит весь мир в узде

Представьте себе FPV-дрон, который должен мгновенно стабилизироваться после резкого маневра уклонения, чтобы нанести точный удар. Или крылатую ракету, огибающую рельеф местности на сверхмалой высоте, где любая ошибка в расчетах приведет к удару о землю. Или гражданскую ступень Falcon 9, совершающую посадку на морскую платформу.

Что объединяет эти процессы? За каждым из них стоит незримая, но железная логика теории устойчивости.

Если вы когда-нибудь занимались теорией управления (Control Theory), вы знаете этот страх: система может пойти вразнос. Чтобы этого не случилось, инженеры и математики используют один мощный инструмент — Уравнение Ляпунова. Это тот самый «серый кардинал» линейной алгебры, который отвечает на главный вопрос жизни, Вселенной и всего такого (в инженерном смысле): «упадет или не упадет?»

Гравитация против матриц: реальный пример

Чтобы понять, зачем нам это уравнение, давайте спустимся с небес на землю (или наоборот, попытаемся взлететь).

Допустим, мы пишем автопилот для удержания ракеты в вертикальном положении.
В механике Ньютона мир описывается уравнениями второго порядка (, где ускорение — это вторая производная координаты $\ddot{z}$ ). Но в теории управления мы не любим вторые производные. Мы любим системы первого порядка.

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

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

$\begin{cases}\frac{dx}{dt} = -y - x^3 \\\frac{dy}{dt} = x - y^3\end{cases}$

Как понять, устойчива ли эта ракета? Вернется ли она в нулевое положение () после порыва ветра?

Попытка №1: линеаризация (путь ленивого инженера)

Любой нормальный инженер сначала попытается систему линеаризовать. Мы отбрасываем «малые» кубические члены ( и ) и смотрим на матрицу Якоби линейной части:

$A = \begin{pmatrix}0 & -1 \\1 & 0\end{pmatrix}$

Теперь найдем собственные числа ( $\lambda$ ) этой матрицы. Это стандартный метод: если у всех $\lambda$ вещественная часть отрицательная ( $\text{Re}(\lambda) < 0$ ), система устойчива.

Решаем характеристическое уравнение $\det(A - \lambda I) = 0$ :

$(-\lambda)(-\lambda) - (-1)(1) = \lambda^2 + 1 = 0$

И получаем:

$\lambda_{1,2} = \pm i$

Вещественная часть равна нулю ( $\text{Re} = 0$ ).

Это катастрофа. В теории устойчивости это называется «критический случай». Линейный анализ говорит нам, что система — это «центр» (вечные колебания по кругу), но он врет. Он слеп к малым нелинейным добавкам. А ведь именно эти и решат судьбу ракеты: либо погасят колебания, либо раскачают их до взрыва.

Нам нужен инструмент мощнее.

Попытка №2: метод Ляпунова (путь опытного исследователя)

Ляпунов предложил гениальную идею: не надо решать сложные дифференциальные уравнения. Надо смотреть на энергию.

Давайте введем функцию , которая будет аналогом полной энергии отклонения нашей ракеты. Самый простой вариант — сумма квадратов ошибок (похоже на потенциальную энергию пружины):

Очевидно, что везде, кроме точки равновесия.
Теперь самое интересное: как меняется эта «энергия» со временем? Найдем полную производную $\dot{V}$ вдоль траекторий нашей системы:

$\frac{dV}{dt} = \frac{\partial V}{\partial x} \cdot \frac{dx}{dt} + \frac{\partial V}{\partial y} \cdot \frac{dy}{dt}$

Подставляем наши уравнения:

$\dot{V} = (2x) \cdot (-y - x^3) + (2y) \cdot (x - y^3)$

Раскрываем скобки:

$\dot{V} = -2xy - 2x^4 + 2yx - 2y^4$

Линейные слагаемые и сокращаются. Это и есть та самая «вечная» осцилляция, которую предсказала матрица . Но посмотрите, что осталось:

$\dot{V} = -2x^4 - 2y^4$

Так как и всегда положительны, то их сумма с минусом всегда отрицательна.

Вердикт: производная энергии отрицательна ( $\dot{V} < 0$ ). Это значит, что, куда бы ни полетела наша ракета, её «энергия» отклонения будет неизбежно убывать, пока не станет нулем.
Система асимптотически устойчива. Ракета не упадет.

В чем подвох?

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

И вот здесь на сцену выходит Уравнение Ляпунова. Оно позволяет не угадывать. Оно позволяет вычислить эту функцию (точнее, ее матрицу) алгебраически, если мы знаем динамику системы.

А если система и так линейная? Зачем тогда страдать?

Вы можете спросить: «Хорошо, для хитрых нелинейных случаев это нужно. Но если я управляю обычным квадрокоптером, который прекрасно описывается линейными уравнениями, зачем мне уравнение Ляпунова? Я же могу просто посчитать собственные числа матрицы (полюса системы). Если они отрицательные — дрон не упадет. Всё?»

Нет, не всё.

Собственные числа отвечают на бинарный вопрос: Да/Нет. (Упадет или не упадет).
Уравнение Ляпунова отвечает на качественный вопрос: Как именно полетит?

Представьте, что вы настраиваете PID-регулятор для того самого дрона.

Вариант А: Дрон возвращается в горизонт медленно, как в киселе. (Устойчив? Да).
Вариант Б: Дрон возвращается мгновенно, но его трясет так, что камера отваливается. (Устойчив? Формально да).

Собственные числа в обоих случаях могут быть отрицательными («устойчив»). Но качество управления — разное. Решение уравнение Ляпунова (матрица ) дает нам интегральную метрику качества. Значение функции в начальный момент времени — это предсказание того, сколько суммарно «боли» (отклонений и затрат энергии) испытает система, пока не успокоится.

Как этим рулить на практике?

Давайте «на пальцах» разберем, как уравнение Ляпунова вшито в мозг любого современного автопилота. У нас есть дрон. Его состояние описывается вектором :

$x = \begin{bmatrix} \text{угол крена} \\ \text{скорость крена} \end{bmatrix}$

Сама по себе физика дрона неустойчива (он хочет перевернуться). Это описывается уравнением $\dot{x} = Ax$ .
Чтобы он летел ровно, мы добавляем управление (напряжение на моторах). Мы используем принцип отрицательной обратной связи:

«Если накренился влево — газуй левыми моторами».

Матрица — это наши настройки (коэффициенты усиления).

Теперь уравнение полета выглядит так:

$\dot{x} = Ax - BKx = (A - BK)x$

И вот тут инженер сталкивается с проблемой.
Мы можем выбрать «злые» коэффициенты , чтобы дрон был резким как пуля.

Но тогда матрица системы может стать такой, что малейший шум датчиков вызовет резонанс. Уравнение Ляпунова выступает здесь как судья:

Мы говорим: «я хочу, чтобы энергия ошибки убывала вот с такой скоростью» (задаем матрицу ).
Подставляем нашу замкнутую систему в уравнение Ляпунова:
Находим .

Если элементы матрицы получаются гигантскими, это значит, что «чаша» энергии очень пологая. Дрону потребуется много времени или дикие усилия моторов, чтобы вернуться в равновесие. Это плохой контроллер. Если компактная — контроллер отличный.

Именно на этой логике строится LQR (Linear Quadratic Regulator) — алгоритм, который управляет всем: от головок жестких дисков до ориентации МКС. Он автоматически ищет такой баланс, чтобы минимизировать решение уравнения Ляпунова.

В чем проблема? Почему нельзя просто «посмотреть»?

Казалось бы, зачем нам сложные матрицы? Хочешь узнать, устойчива ли ракета — запусти симуляцию на 10 минут модельного времени. Если не упала — значит, устойчива.

Но симуляция — это «черный ящик». Она говорит «что» произошло, но не говорит «почему». А главное, она не дает гарантий. Если ракета не упала за 10 минут, упадет ли она за 11? А если ветер подует чуть сильнее?

Уравнение Ляпунова дает алгебраический критерий. Оно превращает вопрос о будущем (динамика, время) в вопрос о структуре (алгебра, матрицы). Решив его один раз, мы получаем математическую гарантию устойчивости для любых начальных условий.

О чем эта статья?

Я взял за основу хорошую, но, к сожалению, краткую англоязычную статью про Lyapunov equation, полностью перевел её и, самое главное, насытил дополнительным контентом:

💻 Python-код: реализация методов и визуализация.
🧠 Интуиция: объяснения «на пальцах», зачем нам это нужно.
🎨 Визуализация: от «шарика в чаше» до векторных полей.
⚙️ Алгоритмическая магия: разбор того, как уйти от чудовищной сложности $\mathcal{O}(n^6)$ к элегантной $\mathcal{O}(n^3)$ с помощью разложения Шура.

Зачем это здесь?

Это пре-релиз статьи для русской Википедии.

Я выкладываю материал на суд сообщества Хабра. Я хочу, чтобы мы вместе сделали лучший материал по этой теме в рунете.

Моя просьба к вам:

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

Итак, наливайте кофе. Мы отправляемся в мир матриц, дифференциальных уравнений, устойчивости и красивых алгоритмов.

Уравнение Ляпунова

Уравнение Ляпунова — это линейное матричное уравнение, играющее фундаментальную роль в теории автоматического управления, дифференциальных уравнениях и анализе устойчивости динамических систем.

Оно названо в честь математика Ляпунова, который в 1892 году в своей докторской диссертации «Общая задача об устойчивости движения» заложил основы современной теории устойчивости.

Памятник Ляпунову в Харькове, где была защищена его легендарная диссертация. Ляпунов был председателем Харьковского математического общества

Уравнение связывает динамику системы (матрицу ) с её энергетическими характеристиками (матрицы и ). Существует две основные формы уравнения: для непрерывного и для дискретного времени.

Непрерывное уравнение Ляпунова

Для линейной системы, описываемой дифференциальным уравнением $\dot{x} = Ax$ , непрерывное уравнение Ляпунова записывается как:

Дискретное уравнение Ляпунова

Для дискретной системы $x_{k+1} = A x_k$ (часто называемое уравнением Стейна) оно имеет вид:

Расшифровка обозначений

* (System Matrix): Квадратная матрица $n \times n,$ описывающая динамику системы (закон, по которому система изменяет свое состояние).

* (Unknown Matrix): Искомая эрмитова (или симметричная) матрица $n \times n$ . В контексте теории устойчивости она определяет функцию Ляпунова которую можно интерпретировать как «обобщенную энергию» системы.

* (Given Matrix): Известная эрмитова матрица $n \times n$ . Она задает скорость диссипации (рассеивания) энергии. Обычно выбирается положительно определенной ().

* (Hermitian Transpose): Эрмитово-сопряженная матрица к (транспонирование с комплексным сопряжением). Если все элементы матрицы — вещественные числа (что чаще всего встречается в инженерных задачах), то заменяется на обычное транспонирование .

Происхождение уравнения

Форма уравнения Ляпунова — это не случайный набор символов, а прямой результат дифференцирования «энергии» системы.

Вывод для непрерывного времени

Представим, что мы хотим проверить устойчивость системы $\dot{x} = Ax.$ Согласно методу Ляпунова, система устойчива, если её «энергия» со временем убывает.

Пусть энергия задается квадратичной формой:

где — положительно определенная матрица.

Найдём скорость изменения этой энергии по времени ( $\frac{d V}{dt}$ ):

$\frac{d}{dt} (x^H P x) = \dot{x}^H P x + x^H P \dot{x}$

Подставим уравнение динамики $\dot{x} = Ax$ и (сопряженное $\dot{x}^H = x^H A^H$ ):

$\frac{d V}{dt} = (x^H A^H) P x + x^H P (Ax) = x^H (A^H P + P A) x$

Чтобы система была устойчивой, энергия должна убывать, то есть производная должна быть отрицательной. Мы требуем, чтобы скорость убывания была равна некоторой заданной величине, определяемой матрицей :

$\frac{d V}{dt} = -x^H Q x$

Сравнивая два выражения, получаем требование к матрицам:

или, в более привычной форме:

(Примечание: Если эрмитова, то , и порядок умножения или зависит от конвенции вектора-строки или столбца, но суть — сумма двух слагаемых — остается неизменной)

Вывод для дискретного времени

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

$V(x_{k+1}) - V(x_k) < 0$

Подставим $x_{k+1} = A x_k$ в формулу энергии :

$V(x_{k+1}) = x_{k+1}^H P x_{k+1} = (A x_k)^H P (A x_k) = x_k^H A^H P A x_k$

Теперь запишем разность энергий:

$\Delta V = x_k^H A^H P A x_k - x_k^H P x_k = x_k^H (A^H P A - P) x_k$

Приравниваем эту разность к отрицательной величине диссипации и получаем уравнение:

$A^H P A - P = -Q \quad \Rightarrow \quad A^H P A - P + Q = 0$

Физический и геометрический смысл

Физическая аналогия: чаша и трение

Представьте, что состояние системы - это положение шарика, катающегося в чаше.

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

Геометрический смысл: эллипсоиды

Геометрически уравнение (где — константа) задает в пространстве состояний -мерный эллипсоид.

В устойчивой системе все траектории движения $\dot{x}$ должны пересекать границы этих эллипсоидов снаружи внутрь.

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

Фазовый портрет асимптотически устойчивой системы.

Красные контуры представляют линии уровня функции Ляпунова , форма которых определяется решением уравнения .

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

Применение к анализу устойчивости

Этот раздел описывает «сердце» применения уравнения Ляпунова. Главная ценность приведенных ниже теорем заключается в том, что они превращают сложную задачу анализа дифференциальных уравнений (где нужно предсказывать будущее системы) в задачу линейной алгебры (где нужно просто решить уравнение).

Вместо того чтобы моделировать систему и смотреть, «упадет» она или нет, мы проводим алгебраический «стресс-тест».

Основная идея: «обратный подход»

Обычно метод функций Ляпунова работает так: мы угадываем функцию энергии и проверяем, убывает ли она. Угадать такую функцию для сложной системы — искусство.

Уравнение Ляпунова позволяет действовать наоборот:

1. Мы заранее требуем, чтобы энергия убывала с определенной скоростью (задаем матрицу диссипации , чаще всего просто единичную матрицу ).

2. Мы спрашиваем математику: «какая форма "чаши" обеспечит такое убывание для нашей системы ?»

3. Мы решаем уравнение относительно .

4. Вердикт: если полученная матрица оказывается положительно определенной (, то есть чаша имеет дно), то система устойчива. Если же не является таковой (например, у неё есть отрицательные собственные числа), то система неустойчива.

Случай непрерывного времени

Для системы, описываемой уравнением $\dot{x} = Ax$ .

Теорема

Система $\dot{x} = Ax$ является глобально асимптотически устойчивой (все траектории сходятся к нулю) тогда и только тогда, когда для любой заданной симметричной положительно определенной матрицы (например, ) решение уравнения

является симметричным и положительно определенным.

Простой пример (скалярный случай)

Пусть наша система — это просто число: $\dot{x} = a x$ .

* Устойчивый случай ():

Уравнение Ляпунова превращается в .

Возьмем (положительное число).

$2(-2)p = -1 \implies -4p = -1 \implies p = 0.25$

Результат: .

Это «правильная» энергия . Система устойчива.

* Неустойчивый случай ():

$2(2)p = -1 \implies 4p = -1 \implies p = -0.25$

Результат: .

Энергия перевернута (холм вместо ямы). Положительно определенного решения не существует. Система неустойчива.

Геометрический смысл

В многомерном случае уравнение проверяет углы между векторами.

* Вектор — это скорость потока.

* Вектор — это нормаль к эллипсоидам энергии.

* Уравнение гарантирует, что во всех точках пространства угол между скоростью и нормалью тупой (косинус <0). Это заставляет поток течь «вниз» по уровням энергии к центру.

Случай дискретного времени

В отличие от непрерывного времени, где система «течет» как вода (плавно изменяя состояние), в дискретном времени система совершает мгновенные «прыжки» или шаги:

$t=0, 1, 2, \dots$

Уравнение динамики:

$x_{k+1} = A x_k$

Здесь матрица — это инструкция для прыжка: «возьми вектор , умножь его на и перенеси в новую точку $x_{k+1}$ ».

Теорема

Линейная система $x_{k+1} = A x_k$ является глобально асимптотически устойчивой (со временем затухает в ноль) тогда и только тогда, когда для любой положительно определенной матрицы решение уравнения:

является положительно определенным ().

Условие устойчивости самой матрицы здесь иное, чем в непрерывном случае:

все собственные числа должны лежать внутри единичного круга на комплексной плоскости ( $|\lambda| < 1$ ).

Почему уравнение выглядит иначе?

В непрерывном мире мы требовали, чтобы скорость изменения энергии была отрицательной ( $\dot{V} < 0$ ). В дискретном мире понятий «скорость» или «производная» нет. Вместо этого мы сравниваем энергию «завтра» и энергию «сегодня».

1. Энергия сегодня: $V_{old} = x_k^T P x_k$ .

2. Энергия завтра (после прыжка):

$V_{new} = x_{k+1}^T P x_{k+1} = (A x_k)^T P (A x_k) = x_k^T (A^T P A) x_k$

3. Изменение энергии за один шаг:

$\Delta V = V_{new} - V_{old} = x_k^T (A^T P A) x_k - x_k^T P x_k$

$\Delta V = x_k^T (\underbrace{A^T P A - P}_{\text{Должно быть } -Q}) x_k$

Мы требуем, чтобы $\Delta V$ было отрицательным (система теряет энергию на каждом шаге). Отсюда и получается уравнение:

Матрица "энергии завтра" минус Матрица "энергии сегодня" равна минус Диссипация.

Простой пример (скалярный случай)

Пусть система — это банковский счет, где сумма умножается на коэффициент каждый месяц: $x_{k+1} = a x_k$ .

Уравнение Ляпунова для скаляров: $p a^2 - p = -q \Rightarrow p(a^2 - 1) = -q$ .

* Устойчивый случай (инфляция съедает деньги)

Пусть (деньги уменьшаются вдвое). Зададим скорость потерь .

$p(0.5^2 - 1) = -1 \implies p(0.25 - 1) = -1 \implies -0.75 p = -1 \implies p \approx 1.33$

Результат: Мы нашли «правильную чашу». Система устойчива.

* Неустойчивый случай (депозит растет): Пусть (деньги удваиваются).

$p(2^2 - 1) = -1 \implies 3p = -1 \implies p = -0.33$

Результат: . Решение отрицательное (перевернутая горка). Устойчивости нет.

Геометрический смысл: прыжки по эллипсам

Представьте себе стадион с беговыми дорожками в форме эллипсов. Самый маленький эллипс — в центре.

* Каждый эллипс — это уровень энергии .

В непрерывной системе бегун плавно смещается с внешней дорожки на внутреннюю по спирали. В дискретной системе бегун телепортируется.

Решить уравнение Ляпунова — значит подобрать такую форму дорожек (матрицу $P$), чтобы при любом прыжке из любой точки бегун гарантированно приземлялся на дорожку, которая находится ближе к центру, чем та, с которой он прыгнул.

Визуализация прыжков

На графике мы покажем траекторию системы. Точки — это положения системы на шагах . Пунктирная линия показывает последовательность. Фон — линии уровня функции Ляпунова. Обратите внимание: в отличие от непрерывного случая, траектория не обязана касаться линий уровня под определенным углом. Она просто должна "перепрыгнуть" через линию уровня внутрь.

Визуализация: сравнение устойчивой и неустойчивой систем

Для иллюстрации давайте создадим график, где решим уравнение Ляпунова для двух случаев.

1. Случай А (устойчивый): Спираль, скручивающаяся внутрь. Матрица задает эллипсы, которые «обнимают» траектории.

2. Случай Б (неустойчивый): Спираль, раскручивающаяся наружу. Если мы попытаемся решить уравнение с , мы получим матрицу , которая не будет положительно определенной (это будет седловая поверхность, будет неопределенной знака), что сигнализирует о неустойчивости.

Методы решения: между красотой и эффективностью

Здесь начинается самое интересное.

Глядя на уравнение , любой, кто проходил курс линейной алгебры, скажет: «Эй, так это же линейное уравнение! Давайте просто выразим ».

И формально он будет прав. Но дьявол, как всегда, кроется в деталях — а точнее, в вычислительной сложности.

Если мы попробуем решить это уравнение «в лоб», как обычную систему , нас ждет катастрофа. Для матрицы размером $n \times n$ сложность такого решения составит $\mathcal{O}(n^6)$ .
Чтобы вы понимали масштаб трагедии: для матрицы $100 \times 100$ (что в инженерии считается небольшим размером) количество операций будет порядка $10^{12}$ .

Ваш ноутбук не скажет вам спасибо.

К счастью, умные люди (Бартелс, Стюарт, Китагава) придумали, как использовать структуру уравнения, чтобы снизить сложность до приемлемых $\mathcal{O}(n^3)$ .

Давайте разберем два подхода:

«наивный» (аналитический) и «профессиональный» (численный).

1. Аналитический подход (или «Ловушка Кронекера»)

Этот метод обожают математики за его красоту и ненавидят программисты за прожорливость.

Идея проста: давайте превратим матрицу в длинный вектор-столбец $\operatorname{vec}(P)$ , просто поставив её столбцы друг под другом. Тогда матричное уравнение превращается в классическую систему линейных уравнений:

$\mathbf{M} \cdot \operatorname{vec}(P) = -\operatorname{vec}(Q)$

Для формирования матрицы $\mathbf{M}$ используется произведение Кронекера ( $\otimes$ ).

Врезка: Что такое произведение Кронекера?

Если обычное умножение матриц — это «строка на столбец», то произведение Кронекера — это «матрица в матрице», своего рода математический фрактал.

Представьте, что у нас есть маленькая матрица размером $2 \times 2$ и любая матрица . Произведение $A \otimes B$ получается так: мы берем матрицу и заменяем каждый её элемент $a_{ij}$ на целую матрицу , умноженную на этот элемент.

$A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} \quad \Rightarrow \quad A \otimes B = \begin{pmatrix} 1 \cdot B & 2 \cdot B \\ 3 \cdot B & 4 \cdot B \end{pmatrix}$

Визуально это выглядит так:

Если тоже размером $2 \times 2$ , то результат раздуется до $4 \times 4$ :

$\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} \otimes \begin{pmatrix} x & y \\ z & w \end{pmatrix} = \left(\begin{array}{cc|cc} x & y & 2x & 2y \\ z & w & 2z & 2w \\ \hline 3x & 3y & 4x & 4y \\ 3z & 3w & 4z & 4w\end{array}\right)$

Именно этот оператор позволяет «распаковать» матричное уравнение в одну длинную систему линейных уравнений. Но цена за это удобство — квадратичный взрыв размеров. Матрица $100 \times 100$ превращается в монстра $10\,000 \times 10\,000$ .

В непрерывном времени () уравнение трансформируется в:

$(I_n \otimes A^T + A^T \otimes I_n) \operatorname{vec}(P) = -\operatorname{vec}(Q)$

В дискретном времени () — в:

$(A^T \otimes A^T - I_{n^2}) \operatorname{vec}(P) = -\operatorname{vec}(Q)$

В чем подвох?

В размере матрицы $\mathbf{M}$ .

Если исходная система имела размерность , то новая система имеет размерность . Матрица коэффициентов при этом становится монстром размером $n^2 \times n^2$ .

При эта матрица содержит элементов.
Решение системы методом Гаусса требует куба от размерности: .

Поэтому этот метод хорош только для доказательства теорем или для совсем крошечных систем (например, ).

2. Взгляд физика: интегралы и ряды

Если система устойчива, решение уравнения Ляпунова можно записать в явном виде, который имеет глубокий физический смысл.

Это сумма (или интеграл) всей «истории» энергии системы.

Непрерывное время

Если матрица гурвицева (все собственные числа имеют отрицательную вещественную часть, $\text{Re}(\lambda) < 0$ ), решение можно представить как интеграл:

$P = \int_0^{\infty} e^{A^T \tau} Q e^{A \tau} d\tau$

Что это значит?
Это накопленная энергия отклика системы на начальное возмущение.
Проверить это легко. Вспомните скалярный случай (где ). Его решение:

$p = \frac{-q}{2a} = \int_0^{\infty} q e^{2a\tau} d\tau$

Матричная формула — это прямое обобщение этого интеграла.

Дискретное время

Если матрица шуровская (все собственные числа лежат внутри единичного круга, $|\lambda| < 1$ ), интеграл заменяется на бесконечную сумму:

$P = \sum_{k=0}^{\infty} (A^T)^k Q A^k$

Это матричный аналог геометрической прогрессии. Для скаляра (или ) решение выглядит как:

$p = \frac{q}{1 - a^2} = \sum_{k=0}^{\infty} q a^{2k}$

Но как считать это на компьютере быстро?
Интегралы и ряды считать долго. Векторизация убивает память.
Здесь на сцену выходит «тяжелая артиллерия» численных методов — алгоритм Бартелса — Стюарта, основанный на разложении Шура. О нем и пойдет речь дальше.

Сравнение методов решения уравнения Ляпунова.

Сверху. График демонстрирует «проклятие размерности». Наивный метод решения через векторизацию (красная линия, $O\left(n^6\right)$ ) становится вычислительно невозможным уже для небольших матриц, в то время как специализированные алгоритмы (зеленая линия, $O\left(n^3\right)$ ) остаются эффективными.

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

3. Численное решение (алгоритм Бартелса — Стюарта)

В инженерной практике — будь то стабилизация ракеты, управление роботом, дроном или балансировка энергосети — размерности матриц легко достигают $1000 \times 1000$ .
Как мы выяснили выше, «наивный» метод здесь потребует $10^{18}$ операций. Если ваш компьютер выполняет миллиард операций в секунду, ждать решения придется около 30 лет.

Нам нужно что-то побыстрее.

Бартелс и Стюарт предложили алгоритм, который стал «золотым стандартом». Именно он работает под капотом, когда вы вызываете lyap в MATLAB, scipy.linalg.solve_continuous_lyapunov в Python или функции библиотеки LAPACK.

Главная идея: «удобная система координат»

Сложность решения системы линейных уравнений напрямую зависит от формы матрицы.

С плотной (заполненной числами) матрицей работать тяжело и дорого.
С квазитреугольной матрицей (где почти половина элементов — нули) уравнения решаются мгновенно, методом простой подстановки.

Суть алгоритма: давайте повернем нашу систему координат так, чтобы страшная матрица стала милой и как можно ближе к треугольной.

Пошаговый разбор (для уравнения )

Шаг 1. Разложение Шура (Schur Decomposition)

Это фундамент метода. Мы ищем такую ортогональную матрицу , которая приведет к верхнему квазитреугольному виду .

: Ортогональная матрица (). Она отвечает за переход в новый базис.
: Верхняя квазитреугольная матрица. На диагонали такой матрицы стоят блоки $1\times1$ (вещественные собственные числа) и $2\times2$ (вещественные блоки, соответствующие комплексно-сопряженным собственным числам). Ниже — нули.

Как это считается?
Матрицы и не падают с неба. Для их нахождения используется итерационный QR-алгоритм (обычно с двойным сдвигом Фрэнсиса). Он последовательно «вычищает» поддиагональные элементы, пока матрица не примет нужную форму.

Если вы хотите разобраться, как именно работает эта магия под капотом (от вращений Гивенса до сходимости), я подробно разобрал механику QR-алгоритма в своем курсе: Численные методы: QR-алгоритм и разложение Шура. Рекомендую заглянуть, если хотите понимать, что происходит внутри scipy.linalg.schur.

Алгоритм приведения к верхнетреугольной форме.

Сначала матрица приводится к верхне-гессенберговской форме

Матрица имеет верхне-гессенбергову форму, если $a_{i j}=0$ , при $i \geq j+2$ .

$H=\left[\begin{array}{lllll}* & * & * & * & * \\* & * & * & * & * \\0 & * & * & * & * \\0 & 0 & * & * & * \\0 & 0 & 0 & * & *\end{array}\right]$

Мы строим отражение так, чтобы оно не трогало первую строку и первый столбец.

Шаг 1. Первый столбец

Вместо того чтобы занулять элементы, начиная со 2-го ( $a_{21}$ ), мы будем занулять, начиная с 3-го ( $a_{31}, a_{41} \dots$ ).
Элемент $a_{21}$ останется (это будет поддиагональ).

Вектор для построения Хаусхолдера берется из части первого столбца:

$\mathbf{x} = \begin{pmatrix} a_{21} \\ a_{31} \\ \vdots \\ a_{n1} \end{pmatrix} \quad \xrightarrow{P_1} \quad \begin{pmatrix} \|\mathbf{x}\| \\ 0 \\ \vdots \\ 0 \end{pmatrix}$

Матрица преобразования будет иметь блочный вид:

$H_1 = \begin{pmatrix} 1 & 0 & \dots \\ 0 & & \\ \vdots & & P_1 & \\ 0 & & \end{pmatrix}$

Где — матрица Хаусхолдера размера $(n-1) \times (n-1)$ .

Что происходит при умножении?

Слева (): действует только на строки $2 \dots n$ . Она создает нули в позициях и ниже.
$H_1 A = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ \mathbf{*} & * & * \\ \mathbf{0} & * & * \end{pmatrix}$
(Жирным выделен измененный первый столбец).
Справа (): Теперь мы умножаем результат на .

Из-за структуры $\begin{pmatrix} 1 & 0 \\ 0 & P_1 \end{pmatrix}$ , это действие смешивает столбцы $2 \dots n$ , но не трогает
первый столбец.
Поэтому наши нули $(\mathbf{0})$ остаются на месте!

Шаг 2. Второй столбец

Теперь мы работаем со вторым столбцом. Мы хотим занулить элементы ниже $a_{32}$ (то есть $a_{42}, a_{52} \dots$ ).
Мы строим матрицу , которая оставляет нетронутыми уже первые две строки и столбца.

$H_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & P_2 \end{pmatrix}$

Умножение слева создаст нули под $a_{32}$ . Умножение справа перемешает столбцы $3 \dots n$ , не затрагивая столбец 1 и 2.

Итог алгоритма

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

$H = \begin{pmatrix} * & * & * & * \\ * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \end{pmatrix}$

Именно с этой формы стартует QR-алгоритм Фрэнсиса (гонка за горбом), так как один шаг QR-итерации для такой матрицы стоит , а не .

Для наглядной демонстрации логики алгоритма привожу код на Питоне

import numpy as np

def householder_reflection(v):
    """Возвращает матрицу отражения, обнуляющую вектор v ниже первого элемента."""
    n = len(v)
    x = np.array(v, dtype=float)
    norm = np.linalg.norm(x)
    if norm == 0:
        return np.eye(n)
    
    # Выбор знака для точности
    sign = np.sign(x[0]) if x[0] != 0 else 1
    x[0] += sign * norm
    x /= np.linalg.norm(x)
    
    return np.eye(n) - 2 * np.outer(x, x)

def to_hessenberg(A):
    """Приводит матрицу A к форме Хессенберга подобием."""
    A_h = A.copy()
    n = A_h.shape[0]
    
    for k in range(n - 2):
        # 1. Берем вектор под диагональю k-го столбца
        # Нам нужно занулить элементы начиная с A[k+2, k]
        x = A_h[k+1:, k]
        
        # 2. Строим отражение
        P_small = householder_reflection(x)
        
        # 3. Расширяем до полного размера (блочная матрица)
        # | I_k+1   0   |
        # |  0    P_small |
        H_k = np.eye(n)
        H_k[k+1:, k+1:] = P_small
        
        # 4. Применяем подобие: H A H^T
        # Умножение слева (работа со строками)
        A_h = H_k @ A_h
        # Умножение справа (работа со столбцами)
        A_h = A_h @ H_k.T
        
    return A_h

# Тест
if __name__ == "__main__":
    np.set_printoptions(precision=4, suppress=True)
    A = np.random.randn(5, 5)
    H = to_hessenberg(A)
    
    print("Исходная матрица A:\n", A)
    print("\nФорма Хессенберга H:\n", H)
    
    # Проверка собственных чисел (должны совпадать)
    print("\nСобств. числа A:", np.sort(np.linalg.eigvals(A).real))
    print("Собств. числа H:", np.sort(np.linalg.eigvals(H).real))

Примечание. О двойном сдвиге Френсиса.

Двойной сдвиг Фрэнсиса (Francis Double Shift) — это гениальный алгоритмический трюк в численно линейной алгебре, который делает QR-алгоритм применимым на практике для вещественных матриц. Если говорить просто: это способ найти комплексные собственные числа вещественной матрицы, не выходя при этом в комплексную арифметику.

В чем проблема обычного QR-алгоритма?

Чтобы найти собственные числа, QR-алгоритм использует сдвиги (shifts) для ускорения сходимости. Обычно в качестве сдвига берется элемент из нижнего правого угла матрицы ( $a_{nn}$ ). Но у вещественных матриц часто бывают комплексные собственные числа (они всегда идут парами: и ).

Если мы хотим найти такое число, нам нужно сделать сдвиг на комплексное число $\sigma$ .
Как только мы вычитаем комплексное $\sigma$ из вещественной матрицы ( $A - \sigma I$ ), вся матрица становится комплексной.
Последствия: Объем памяти вырастает в 2 раза, количество операций умножения — в 4 раза. Это дорого и неэффективно.

Идея Фрэнсиса: «Ходим парами»

Джон Фрэнсис (John Francis) придумал, как обойти это ограничение.
Вместо одного шага со сдвигом $\sigma$ , давайте сделаем сразу два шага подряд: один со сдвигом $\sigma$ , а второй — с его комплексно-сопряженным $\bar{\sigma}$ .

Математически это эквивалентно работе с полиномом от матрицы:

$M = (A - \sigma I)(A - \bar{\sigma} I)$

Раскроем скобки:

$M = A^2 - (\sigma + \bar{\sigma})A + (\sigma \bar{\sigma}) I$

Магия:

Сумма сопряженных чисел $(\sigma + \bar{\sigma})$ — это вещественное число ( $2 \cdot \text{Re}(\sigma)$ ).
Произведение $(\sigma \bar{\sigma})$ — это тоже вещественное число ( $|\sigma|^2$ ).
Значит, матрица остается полностью вещественной, даже если сдвиги были мнимыми!

Реализация: «Погоня за горбом» (Bulge Chasing)

Фрэнсис придумал, как не перемножать матрицы целиком (это было бы долго, ). Он использовал теорему о неявном Q (Implicit Q Theorem).

Мы вычисляем только первый столбец воображаемой матрицы .
Это создает «возмущение» (горб, bulge) в верхней части матрицы Хессенберга.
Затем мы с помощью вращений Хаусхолдера или Гивенса «прогоняем» этот горб вниз по диагонали, пока он не исчезнет в правом нижнем углу.

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

Итог

Благодаря двойному сдвигу Фрэнсиса:

Алгоритм работает очень быстро (кубическая сходимость).
Мы остаемся в вещественных числах ( памяти вместо ).
На выходе получается та самая матрица Шура с блоками $1\times1$ (вещественные с.ч.) и $2\times2$ (комплексные пары), которая нужна нам для алгоритма Бартелса — Стюарта.

Именно этот алгоритм работает внутри функции scipy.linalg.schur и делает возможным быстрое решение уравнения Ляпунова.

Как это работает с геометрической точки зрения?
Алгоритм использует преобразования Хаусхолдера (матрицы отражения).
Первое отражение создает «возмущение» (ненулевой горб) в левом верхнем углу матрицы, кодирующее информацию о двойном сдвиге.
Серия последующих отражений «гонит» этот горб вниз по диагонали, восстанавливая треугольную структуру.
Этот процесс (Bulge Chasing) позволяет выполнить сложную математическую операцию над всей матрицей, затрагивая на каждом шаге только 3 строки и 3 столбца, что обеспечивает высочайшую скорость $\mathcal{O}(n^3)$ .
Матрица Хаусхолдера - это матрица вида $H \equiv H(v)=I-2 v v^*$ , где $v-n \times 1$ вектор и . Это также матрица отражения относительно плоскости с нормалью : $H x=x-2\left(v^* x\right) v$

Шаг 2. Трансформация уравнения

Берем наше уравнение:

Подставляем в него (и помним, что ):

Теперь применим «сэндвич-трюк»: умножим все уравнение слева на , а справа на . Благодаря свойству ортогональности (), внешние матрицы «аннигилируют», а внутренние группируются:

$\underbrace{T^T}_{\text{Треуг.}} \underbrace{(U^T P U)}_{\tilde{P}} + \underbrace{(U^T P U)}_{\tilde{P}} \underbrace{T}_{\text{Треуг.}} = -\underbrace{(U^T Q U)}_{\tilde{Q}}$

Шаг 3. Замена переменных

Мы получили новое уравнение относительно новой неизвестной матрицы $\tilde{P}$ :

$T^T \tilde{P} + \tilde{P} T = -\tilde{Q}$

Здесь:

$\tilde{P} = U^T P U$ — искомая матрица в повернутом базисе.
$\tilde{Q} = U^T Q U$ — пересчитанная правая часть (легко вычисляется).

Шаг 4. Решение квазитреугольной системы («Матричное домино»)

Это самый красивый момент. Казалось бы, уравнение $T^T \tilde{P} + \tilde{P} T = -\tilde{Q}$ выглядит так же сложно, как исходное. Но нет!

Поскольку — квазитреугольная матрица, мы можем находить элементы $\tilde{P}$ по одному (или блоками по 2 элемента), начиная с углов.

Мы вычисляем один элемент, подставляем его, вычисляем соседний — и так, по цепочке, распутываем весь клубок.

Вместо решения одной гигантской системы $N^2 \times N^2$ , мы решаем последовательность из крошечных уравнений.

Это снижает сложность с $\mathcal{O}(n^6)$ до $\mathcal{O}(n^3)$ .

В вычислительной математике уравнение такого типа называется уравнением Сильвестра.

Уравнение Сильвестра — это уравнение вида . Оно является более общим случаем уравнения Ляпунова (где не обязательно равно ).

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

Теперь разберем подробно, как выглядит это решение.

Итак, мы получили уравнение $T^T \tilde{P} + \tilde{P} T = -\tilde{Q}$ .

Матрица — квазитреугольная. Это значит, что мы можем найти матрицу $\tilde{P}$ последовательно, блок за блоком.

Мы идем с нижнего правого угла () к верхнему левому ().
Пусть мы хотим найти блок $\tilde{P}_{ij}$ (размером $1\times1$ , $1\times2$ , $2\times1$ или $2\times2$ ).

Запишем уравнение только для этого блока:

$T_{ii}^T \tilde{P}_{ij} + \tilde{P}_{ij} T_{jj} = \underbrace{-\tilde{Q}_{ij} - \sum_{k=i+1}^{N} T_{ki} \tilde{P}_{kj} - \sum_{k=j+1}^{N} \tilde{P}_{ik} T_{kj}}_{R_{ij} \text{ (Уже известно)}}$

Взгляните на правую часть ( $R_{ij}$ ). Все слагаемые в суммах содержат индексы или . Поскольку мы идем с конца ( $N \to 1$ ), эти значения нам уже известны! Мы просто переносим их в правую часть как константы.

Нам остается решить маленькое линейное уравнение относительно $\tilde{P}_{ij}$ :

$T_{ii}^T \tilde{P}_{ij} + \tilde{P}_{ij} T_{jj} = R_{ij}$

Здесь возможны всего три сценария:

Сценарий А: Скаляр + Скаляр ( $1 \times 1$ )
Если $T_{ii}$ и $T_{jj}$ — обычные числа (вещественные собственные числа):

$t_{ii} \cdot p_{ij} + p_{ij} \cdot t_{jj} = r_{ij}$ $p_{ij} (t_{ii} + t_{jj}) = r_{ij} \quad \Rightarrow \quad p_{ij} = \frac{r_{ij}}{t_{ii} + t_{jj}}$

Сценарий Б: Скаляр + Блок $2\times2$ (или наоборот)

Это уравнение Сильвестра для прямоугольной матрички $1\times2$ . Это система из 2-х линейных уравнений с 2-мя неизвестными. Решается элементарно методом Крамера или подстановкой.

Сценарий В: Блок $2\times2$ + Блок $2\times2$

Самый сложный случай (когда сталкиваются два комплексно-сопряженных собственных числа). Мы ищем матрицу $\tilde{P}_{ij}$ размером $2\times2$ (4 неизвестных).

Уравнение вида для матриц $2\times2$ можно развернуть через векторизацию (Кронекера) в систему $4\times4$ :

$(I_2 \otimes T_{ii}^T + T_{jj}^T \otimes I_2) \operatorname{vec}(\tilde{P}_{ij}) = \operatorname{vec}(R_{ij})$

Решить систему линейных уравнений $4\times4$ для компьютера — это наносекунды.

Шаг 5. Возвращение домой

Найдя $\tilde{P}$ , мы возвращаемся в исходные координаты обратным преобразованием:

$P = U \tilde{P} U^T$

Вуаля! Задача решена, процессор холоден, ракета летит стабильно.

4. Дискретный мир: Алгоритм Китагавы

Если в непрерывном времени мы решаем дифференциальные уравнения, то в цифровых системах управления (DSP, компьютерное зрение, экономика) мы живем в мире дискретных шагов.
Уравнение здесь меняется на:

Идеологически алгоритм решения остается прежним (разложение $\to$ трансформация $\to$ подстановка), но алгебра меняется. Этот метод часто называют алгоритмом Китагавы.

Шаг 1. Разложение Шура (Фундамент)

Здесь всё без изменений. Мы находим ортогональную матрицу и верхнюю квазитреугольную :

Шаг 2. Трансформация уравнения

Подставим разложение в наше дискретное уравнение .
Помним, что .

Снова применяем «сэндвич-трюк»: умножаем всё уравнение слева на и справа на .

В первом слагаемом () внешние матрицы и аннигилируют с домноженными, оставляя «начинку».
Во втором слагаемом () и правой части () образуются новые переменные в повернутом базисе.

$T^T \underbrace{(U^T P U)}_{\tilde{P}} T - \underbrace{(U^T P U)}_{\tilde{P}} = -\underbrace{(U^T Q U)}_{\tilde{Q}}$

В итоге получаем красивое уравнение для новой матрицы $\tilde{P}$ :

$T^T \tilde{P} T - \tilde{P} = -\tilde{Q}$

Шаг 3. Решение квазитреугольной системы

Теперь перед нами уравнение, где — верхняя квазитреугольная матрица. Как и в методе Бартелса — Стюарта, мы используем обратную подстановку, начиная с нижнего правого угла. Давайте посмотрим, во что вырождается это уравнение для самого последнего диагонального элемента . Из-за квазитреугольной структуры матрицы , элемент в углу зависит только от самого себя:

$t_{nn} \cdot \tilde{p}_{nn} \cdot t_{nn} - \tilde{p}_{nn} = -\tilde{q}_{nn}$

Выносим $\tilde{p}_{nn}$ за скобки:

$\tilde{p}_{nn} (t_{nn}^2 - 1) = -\tilde{q}_{nn}$

Отсюда мгновенно получаем решение:

$\tilde{p}_{nn} = \frac{-\tilde{q}_{nn}}{t_{nn}^2 - 1}$

В чем фундаментальное отличие?

В непрерывном случае: В знаменателе было $2t_{nn}$ . Деление на ноль происходило, если собственное число $\lambda = 0$ .
В дискретном случае: В знаменателе стоит $t_{nn}^2 - 1$ . Деление на ноль происходит, если $t_{nn}^2 = 1$ , то есть $|\lambda|=1$ . Это блестяще согласуется с теорией устойчивости: дискретная система теряет устойчивость, когда собственные числа выходят на единичную окружность. Алгоритм буквально «ломается» ровно в тот момент, когда система становится неустойчивой.

Шаг 4. Эффект домино

Вычислив угловой элемент $\tilde{p}_{nn}$ , мы поднимаемся к соседним (например, $\tilde{p}_{n-1, n}$ ). Уравнение для них снова становится линейным с одним неизвестным, так как все "нижние" значения уже найдены. Мы распутываем матрицу с конца к началу.

Шаг 5. Восстановление

Возвращаемся в исходные координаты:

$P = U \tilde{P} U^T$

Шпаргалка: Сравнение методов

Чтобы уложить все в голове, давайте сведем два метода в одну таблицу.

Характеристика	Непрерывный случай (Bartels — Stewart)	Дискретный случай (Kitagawa)
Уравнение
Физический смысл	Баланс скорости изменения энергии	Баланс потери энергии за один шаг
Трансформированное уравнение	$T^T \tilde{P} + \tilde{P} T = -\tilde{Q}$	$T^T \tilde{P} T - \tilde{P} = -\tilde{Q}$
Ядро алгоритма (для угла)	$\tilde{p}_{nn} = \frac{-\tilde{q}_{nn}}{2 t_{nn}}$	$\tilde{p}_{nn} = \frac{-\tilde{q}_{nn}}{t_{nn}^2 - 1}$
Опасная зона (Сингулярность)	$\lambda_i + \lambda_j = 0$ (Симметрия относительно мнимой оси)	$\lambda_i \lambda_j = 1$ (Симметрия относительно единичной окружности)
Сложность	$\mathcal{O}(n^3)$	$\mathcal{O}(n^3)$

Сравнение условий разрешимости (скалярный случай).

Графики показывают зависимость решения $y_{nn}$ от собственного числа матрицы $\lambda_{nn}$ .

Слева (Непрерывный): решение уходит в бесконечность (сингулярность) при $\lambda = 0$ . Это соответствует нейтральной устойчивости (система не меняет состояние).

Справа (Дискретный): решение имеет две точки разрыва при $\lambda = 1$ и $\lambda = -1$ . Это границы единичного круга на комплексной плоскости. Если собственное число попадает на эту границу, дискретное уравнение Ляпунова не имеет единственного решения.

Связь миров: от дискретного к непрерывному

На первый взгляд может показаться, что непрерывное уравнение () и дискретное () — это две разные математические вселенные. Но на самом деле они — близкие родственники.

Дискретная система — это, по сути, «покадровая съемка» непрерывного процесса. Если мы начнем уменьшать время между кадрами (шаг $\delta \to 0$ ), то дискретное уравнение должно плавно превратиться в непрерывное.

Давайте проверим это математически, используя старый добрый метод Эйлера.

1. От динамики к статике

Возьмем линейную непрерывную систему:

$\dot{\mathbf{x}}(t) = \mathbf{A} \mathbf{x}(t)$

Чтобы загнать её в компьютер, мы дискретизируем время с шагом $\delta$ . Заменим производную конечной разностью:

$\frac{\mathbf{x}_{k+1} - \mathbf{x}_k}{\delta} \approx \mathbf{A} \mathbf{x}_k$

Выразим состояние на следующем шаге:

$\mathbf{x}_{k+1} \approx \mathbf{x}_k + \delta \mathbf{A} \mathbf{x}_k = (\mathbf{I} + \delta \mathbf{A}) \mathbf{x}_k$

Мы получили классическую дискретную систему $\mathbf{x}_{k+1} = \mathbf{B} \mathbf{x}_k$ , где матрица перехода за один шаг:

$\mathbf{B} = \mathbf{I} + \delta \mathbf{A}$

2. Трансформация уравнения Ляпунова

Теперь возьмем дискретное уравнение Ляпунова для нашей новой матрицы $\mathbf{B}$ .

Важный физический нюанс:

В непрерывном уравнении матрица $\mathbf{Q}$ задает мощность диссипации (Джоули в секунду).

В дискретном уравнении правая часть — это энергия, потерянная за один шаг.

Поскольку шаг длится $\delta$ секунд, то энергия за шаг примерно равна $\mathbf{Q} \cdot \delta$ .

Поэтому правая часть дискретного уравнения должна быть масштабирована на $\delta$ :

$\mathbf{B}^T \mathbf{P} \mathbf{B} - \mathbf{P} = -\delta \mathbf{Q}$

Подставим сюда наше выражение для $\mathbf{B} = (\mathbf{I} + \delta \mathbf{A})$ :

$(\mathbf{I} + \delta \mathbf{A})^T \mathbf{P} (\mathbf{I} + \delta \mathbf{A}) - \mathbf{P} = -\delta \mathbf{Q}$

Помним, что $(\mathbf{I} + \delta \mathbf{A})^T = \mathbf{I} + \delta \mathbf{A}^T$ . Раскроем скобки, удерживая в голове, что $\mathbf{I}\mathbf{P} = \mathbf{P}$ :

$\left( \mathbf{P} + \delta \mathbf{P}\mathbf{A} + \delta \mathbf{A}^T \mathbf{P} + \delta^2 \mathbf{A}^T \mathbf{P} \mathbf{A} \right) - \mathbf{P} = -\delta \mathbf{Q}$

Матрицы $\mathbf{P}$ и $-\mathbf{P}$ уничтожаются. Остается:

$\delta (\mathbf{A}^T \mathbf{P} + \mathbf{P} \mathbf{A}) + \delta^2 (\mathbf{A}^T \mathbf{P} \mathbf{A}) = -\delta \mathbf{Q}$

3. Магия малых величин

Разделим всё уравнение на $\delta$ (считая, что $\delta > 0$ ):

$\mathbf{A}^T \mathbf{P} + \mathbf{P} \mathbf{A} + \delta (\mathbf{A}^T \mathbf{P} \mathbf{A}) = -\mathbf{Q}$

А теперь — финал. Устремим шаг времени к нулю ( $\delta \to 0$ ). Мы переходим от «слайд-шоу» к плавному кино.
Слагаемое $\delta (\mathbf{A}^T \mathbf{P} \mathbf{A})$ , содержащее множитель $\delta$ , исчезает.

В пределе получаем:

$\mathbf{A}^T \mathbf{P} + \mathbf{P} \mathbf{A} = -\mathbf{Q}$

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

Иллюстрация: сходимость на практике

Лучше один раз увидеть, чем сто раз продифференцировать.
На графике ниже показана ошибка сходимости. Мы берем одну и ту же систему, решаем для нее непрерывное уравнение (получаем «идеальную» ), а затем решаем дискретное уравнение с уменьшающимся шагом $\delta$ . График в логарифмическом масштабе — это прямая линия, уходящая вниз. Это означает, что разница между матрицами $||\mathbf{P}_{discrete} - \mathbf{P}_{continuous}||$ стремится к нулю линейно вместе с $\delta$ .

Секретное оружие: сжатие реальности

Уравнение Ляпунова умеет делать еще один фокус, за который его обожают инженеры больших систем. Оно позволяет сжимать гигабайтные модели в килобайтные формулы.

Представьте, что вы моделируете нагрев процессора. У вас есть сетка из миллиона узлов. Ваша матрица имеет размер $10^6 \times 10^6$ . Решать такое в реальном времени невозможно.

Мы решаем два уравнения Ляпунова:

Для Управляемости (Controllability Gramian ): показывает, какие состояния системы легко «раскачать» входом.
Для Наблюдаемости (Observability Gramian ): показывает, какие состояния сильнее всего влияют на выход.

Дальше работает магия метода, называемого сбалансированное усечение (Balanced Truncation). Мы оставляем только те состояния, которые одновременно и легко управляются, и хорошо наблюдаются (у них большие собственные числа в произведении ). Всё остальное — математический шум, который можно выкинуть.

Так из монструозной модели получается элегантная аппроксимация, которая работает в 1000 раз быстрее и почти так же точно.

Код для самостоятельных проектов и иллюстраций.

Я собрал все описанные методы в один класс.

P.S. Если этот материал оказался полезным, и вы хотите видеть больше подобных разборов для Википедии — дайте знать. Наука должна быть открытой и понятной.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import solve_continuous_lyapunov, solve_discrete_lyapunov, schur, norm, expm
import time

class LyapunovLab:
    """
    Лаборатория для визуализации и решения уравнений Ляпунова.
    Включает методы для анализа устойчивости, геометрии и алгоритмов.
    """

    def __init__(self):
        # Настройки стилей графиков
        plt.rcParams['figure.figsize'] = (10, 6)
        plt.rcParams['axes.grid'] = True
        plt.rcParams['grid.alpha'] = 0.3
        print("LyapunovLab инициализирована. Готовы к экспериментам.")

    def run_all(self):
        """Запускает все демо-примеры последовательно."""
        self.demo_stability_bowl()
        self.demo_geometric_vectors()
        self.demo_schur_structure()
        self.demo_benchmark_complexity()
        self.demo_discrete_convergence()
        plt.show()

    def demo_stability_bowl(self):
        """
        Демонстрация 1: Физическая аналогия (Чаша).
        Показывает функцию Ляпунова как поверхность, по которой скатывается состояние.
        """
        print("\n--- Демо 1: Физическая аналогия (Чаша) ---")
        # Устойчивая матрица (спиральный сток)
        A = np.array([[-0.5,  2.0],
                      [-2.0, -1.0]])
        Q = np.eye(2)
        
        # Решаем уравнение: A.T * P + P * A = -Q
        P = solve_continuous_lyapunov(A.T, -Q)
        
        # Сетка
        x = np.linspace(-2, 2, 30)
        y = np.linspace(-2, 2, 30)
        X, Y = np.meshgrid(x, y)
        
        # V(x) = x^T P x
        Z = P[0,0]*X**2 + (P[0,1]+P[1,0])*X*Y + P[1,1]*Y**2
        
        fig = plt.figure(figsize=(10, 7))
        ax = fig.add_subplot(111, projection='3d')
        
        ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8, edgecolor='none')
        ax.plot_wireframe(X, Y, Z, color='black', alpha=0.1, rstride=3, cstride=3)
        
        # Траектория "скатывания"
        from scipy.integrate import odeint
        t = np.linspace(0, 10, 200)
        path = odeint(lambda state, t: A @ state, [1.8, 1.8], t)
        path_z = [vec @ P @ vec for vec in path] # Высота траектории
        
        ax.plot(path[:,0], path[:,1], path_z, color='red', lw=2, label='Траектория системы')
        ax.scatter([0], [0], [0], color='red', s=50, label='Равновесие')
        
        ax.set_title("Функция Ляпунова: Энергетическая чаша\n$V(x) = x^T P x$")
        ax.set_xlabel('$x_1$'); ax.set_ylabel('$x_2$'); ax.set_zlabel('Energy')
        ax.legend()
        plt.tight_layout()

    def demo_geometric_vectors(self):
        """
        Демонстрация 2: Геометрия (Тупой угол).
        Показывает, что вектор скорости всегда смотрит внутрь эллипса энергии.
        """
        print("\n--- Демо 2: Геометрия векторов ---")
        A = np.array([[-1.0,  1.5],
                      [-1.5, -1.0]])
        Q = np.eye(2)
        P = solve_continuous_lyapunov(A.T, -Q)
        
        # Выбираем точку на эллипсе
        x0 = np.array([1.5, 1.0])
        # Нормируем точку под уровень энергии C=3
        current_energy = x0 @ P @ x0
        x0 = x0 * np.sqrt(3.0 / current_energy)
        
        # Векторы
        v_vel = A @ x0   # Скорость (Ax)
        v_norm = P @ x0  # Нормаль (Px)
        
        fig, ax = plt.subplots(figsize=(6, 6))
        
        # Рисуем эллипс
        limit = 2.5
        vals = np.linspace(-limit, limit, 100)
        X, Y = np.meshgrid(vals, vals)
        Z = P[0,0]*X**2 + (P[0,1]+P[1,0])*X*Y + P[1,1]*Y**2
        ax.contour(X, Y, Z, levels=[3.0], colors='black', linewidths=2)
        
        # Рисуем векторы
        ax.quiver(*x0, *v_norm, color='crimson', scale=10, label='Нормаль (Px)')
        ax.quiver(*x0, *v_vel, color='navy', scale=10, label='Скорость (Ax)')
        
        ax.set_title("Геометрия: Скорость против Нормали\nУгол > 90° обеспечивает устойчивость")
        ax.set_aspect('equal')
        ax.legend()
        
    def demo_benchmark_complexity(self):
        """
        Демонстрация 3: Сложность O(n^6) vs O(n^3).
        Сравнение наивного метода и метода Бартелса-Стюарта.
        """
        print("\n--- Демо 3: Бенчмарк (Битва алгоритмов) ---")
        
        # Наивный решатель (через Кронекера)
        def naive_solve(A_in, Q_in):
            n = A_in.shape[0]
            # vec(AX + XA^T) = (I x A + A x I) vec(X)
            # Обратите внимание: для XA^T это (A x I), для AX это (I x A)
            M = np.kron(np.eye(n), A_in) + np.kron(A_in, np.eye(n))
            return np.linalg.solve(M, -Q_in.reshape(-1)).reshape(n, n)

        sizes = range(2, 21, 2) # Небольшие размеры, т.к. наивный метод очень медленный
        times_naive = []
        times_smart = []
        
        print(f"{'N':<5} | {'Naive (s)':<10} | {'Smart (s)':<10}")
        for n in sizes:
            A = -np.eye(n) + 0.1*np.random.randn(n, n)
            Q = np.eye(n)
            
            t0 = time.time()
            naive_solve(A, Q)
            times_naive.append(time.time() - t0)
            
            t0 = time.time()
            solve_continuous_lyapunov(A.T, -Q)
            times_smart.append(time.time() - t0)
            
            print(f"{n:<5} | {times_naive[-1]:<10.4f} | {times_smart[-1]:<10.4f}")

        fig, ax = plt.subplots()
        ax.plot(sizes, times_naive, 'r-o', label='Naive $O(n^6)$')
        ax.plot(sizes, times_smart, 'g-s', label='Bartels-Stewart $O(n^3)$')
        ax.set_title("Взрыв сложности: Почему важен алгоритм")
        ax.set_xlabel("Размер матрицы N")
        ax.set_ylabel("Время (сек)")
        ax.legend()
        ax.set_yscale('log')

    def demo_schur_structure(self):
        """
        Демонстрация 4: Матрица Шура.
        Визуализация того, как A превращается в T.
        """
        print("\n--- Демо 4: Разложение Шура ---")
        # Создаем матрицу с комплексно-сопряженными числами
        T_true = np.array([[ -1,  2, 0, 0],
                           [ -2, -1, 0, 0], # Блок 2x2
                           [  0,  0,-3, 1],
                           [  0,  0, 0,-4]])
        # Запутываем её
        Q_orth, _ = np.linalg.qr(np.random.randn(4, 4))
        A = Q_orth @ T_true @ Q_orth.T
        
        T, U = schur(A, output='real')
        
        fig, axes = plt.subplots(1, 3, figsize=(12, 4))
        
        for ax, mat, title in zip(axes, [A, U, T], ['Исходная A', 'Ортогональная U', 'Форма Шура T']):
            im = ax.imshow(mat, cmap='Blues')
            ax.set_title(title)
            fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
            
            # Подсветим нули на T
            if title == 'Форма Шура T':
                # Рисуем треугольник нулей
                n = mat.shape[0]
                poly = patches.Polygon([[-0.5, 1.5], [-0.5, n-0.5], [n-2.5, n-0.5]], 
                                       closed=True, color='white', alpha=0.5, hatch='///')
                ax.add_patch(poly)
                ax.text(1, 2.5, "Zeros", color='gray', fontweight='bold', rotation=45)

        plt.suptitle("Визуализация разложения Шура ($A = U T U^T$)", y=1.02)
        plt.tight_layout()

    def demo_discrete_convergence(self):
        """
        Демонстрация 5: Связь времен.
        Показывает, что при delta -> 0 решение дискретного уравнения сходится к непрерывному.
        """
        print("\n--- Демо 5: Сходимость (Дискретное -> Непрерывное) ---")
        A = np.array([[-0.5, 2.0], [-2.0, -1.0]])
        Q = np.eye(2)
        
        # Истинное непрерывное решение
        P_cont = solve_continuous_lyapunov(A.T, -Q)
        
        deltas = np.logspace(-1, -4, 15)
        errors = []
        
        for d in deltas:
            # Дискретный аналог: x_{k+1} = (I + dA)x_k
            A_disc = np.eye(2) + d * A
            # Дискретная диссипация энергии ~ d * Q
            Q_disc = d * Q
            
            # Решаем дискретное (нужно передать A.T, т.к. scipy решает AXA^H - X = -Q)
            try:
                P_disc = solve_discrete_lyapunov(A_disc.T, -Q_disc)
                errors.append(norm(P_disc - P_cont, 'fro'))
            except:
                errors.append(np.nan)

        fig, ax = plt.subplots()
        ax.loglog(deltas, errors, 'm-o', lw=2, label='Ошибка $|P_{disc} - P_{cont}|$')
        ax.set_title("Сходимость дискретного решения к непрерывному")
        ax.set_xlabel("Шаг времени $\delta$ (log)")
        ax.set_ylabel("Ошибка (log)")
        ax.grid(True, which="both", linestyle='--')
        ax.legend()

# --- Точка входа ---
if __name__ == "__main__":
    lab = LyapunovLab()
    lab.run_all()

Этот код проверен и работает. Теперь вы можете экспериментировать самостоятельно!

Заключение

Мы начали с того, что уравнение Ляпунова — это ответ на вопрос «упадет или не упадет». Но, погрузившись глубже, мы увидели, что это фундамент гораздо большего здания.

Оно гарантирует безопасность там, где интуиция и линейная алгебра бессильны.
Оно оценивает качество управления, позволяя создавать идеальные автопилоты
Оно сжимает сложность, позволяя моделировать физику на слабых процессорах.
И, наконец, оно просто красиво с алгоритмической точки зрения, объединяя в себе матрицы, дифференциальные уравнения и численные методы.

Теперь, когда вы увидите в научной статье фразу «пусть — решение уравнения Ляпунова», вы будете знать: за этими словами скрывается не просто сухая формула, а мощнейший инструмент укрощения хаоса.