In [1]:
%pylab inline
from nesode import *
rcParams['figure.figsize']=(6,6)
from itertools import product
Populating the interactive namespace from numpy and matplotlib

Дифференциальные уравнения

Совместный бакалавриат ВШЭ-РЭШ, 2013-14 учебный год

И. В. Щуров, П. Ф. Соломатин, И. А. Хованская, А. Петрин, Н. Солодовников

Многомерные линейные уравнения с постоянными коэффициентами

Рассмотрим линейное уравнение $$\tag{1} \dot z=Az,\quad z(t)\in\mathbb R^n$$

Ранее мы разобрали случаи уравнений на плоскости (когда $n=2$) почти полностью (за исключением пары вырожденых случаев). Сегодня мы рассмотрим общий случай.

Когда мы рассматривали линейные системы на плоскости, у нас получалось, что решение задается в виде

$$z(t)=M(t) z^0,$$

где $M(t)$ — некоторая матрица (зависящая от $t$), к тому же $M(0)=E$ (единичная матрица). Это неспроста.

Пусть $\varphi(t;z^0)$ — решение уравнения (1) с начальным условием $z(0)=z^0$. Зафиксируем некоторое $t$ и рассмотрим отображение последования (оно также называется преобразованием фазового потока) за время $t$:

$$g^t(z^0)=\varphi(t;z^0)$$

То есть для каждой точки фазового пространства $z^0$ отобразим её туда, где она окажется через время $t$.

Пример. Рассмотрим систему $\dot x = 2x, \quad \dot y = 3y$, тогда $g^1\begin{pmatrix} x_0 \\ y_0 \end{pmatrix} = \begin{pmatrix} e^2 & 0 \\ 0 & e^3 \end{pmatrix} $, а в общем виде $$g^t\begin{pmatrix} x_0 \\ y_0 \end{pmatrix} = \begin{pmatrix} e^{2t} & 0 \\ 0 & e^{3t} \end{pmatrix} \begin{pmatrix} x_0\\ y_0 \end{pmatrix}$$ или просто $$g^t=\begin{pmatrix} e^{2t} & 0 \\ 0 & e^{3t} \end{pmatrix} $$

Теорема. Для системы (1) и любого $t$ отображение последования является линейным.

Доказательство. Пусть $u$ и $v$ — некоторые векторы из фазового пространства $\mathbb R^n$. Тогда

$$g^t(u+v)=\varphi(t;u+v)$$

Заметим, что $\psi(t)=\varphi(t;u)+\varphi(t;v)$ — решение уравнения (1) (в силу линейности). При этом $\psi(0)=\varphi(0;u)+\varphi(0;v)=u+v$. Значит, это решение с начальным условием $u+v$. Но решение с начальным условием $u+v$ задаётся функцией $\varphi(t;u+v)$. Значит (по теореме о существовании и единственности), $\psi(t)=\varphi(t;u+v)$. Таким образом,

$$g^t(u+v)=\varphi(t;u)+\varphi(t;v)=g^t (u)+g^t (v)$$

Аналогично проверяется и вторая аксиома линейности для $g^t$. (Упраженние: провести его.)

Доказательство завершено.

Следствие. Пространство решений линейного дифференциального уравнения имеет такую же размерность, как фазовое простнаство. (Вообще-то, и не только линейного.)

Матричная экспонента

Таким образом, для решения уравнения (1) нам нужно будет найти матрицу $M(t)$ (она иногда назыавется «матрицей монодромии»), которая задаёт решение. Как её найти? Если бы $A$ была не матрицей, а числом (вещественным и комплексным), мы бы мгновенно записали решение в виде экспоненты. Нельзя ли с матрицей сделать то же самое, то есть записать решение уравнения (1) в виде экспоненты от матрицы?

$$\tag{2} z(t)=e^{At} z^0$$

На первый взгляд, это кажется безумием. (Хотя возможно вас с ним уже познакомили на курсе линейной алгебры.) Что значит «возвести число $e$ в степень матрицы»? Ерунда какая-то. Впрочем, не большая ерунда, чем возведение числа $e$ в степень $\pi$ (изначально возвести число в степень — это умножить его на себя сколько-то раз).

Напомним, что экспоненту от числа $x$ можно определить как сумму ряда

$$e^x=\sum_{n=0}^{\infty} \frac{x^n}{n!}$$

Оказывается, нет ничего невозможного в том, чтобы подставить в этот ряд матрицу $A$. Действительно, чтобы посчитать значение этого ряда, нам нужно только уметь складывать матрицы и возводить их в натуральные степени — а это мы делать умеем. Итак, по определению,

$$e^A=\sum_{n=0}^{\infty} \frac{A^n}{n!}, \quad A^0=E$$

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

Матричная экспонента даёт решение линейного уравнения

Покажем, что (2) даёт решение уравнения (1). Действительно,

$$\frac{d}{dt}e^{At}z^0=\frac{d}{dt}\sum_{n=0}^{\infty} \frac{t^n A^n z^0}{n!}=\sum_{n=1}^\infty \frac{nt^{n-1} A^n z^0}{n!}=A\sum_{n=0}^\infty \frac{t^n A^n}{n!}=Ae^{At}z^0$$

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

Контрольный вопрос.

Правда ли, что $e^{A+B}=e^A e^B$?

Ответ: нет, если только матрицы $A$ и $B$ не коммутируют.

Нахождение матричной экспоненты

Утверждение: $e^{CAC^{-1}}=Ce^{A}C{-1}$.

Доказательство. $$Ce^AC^{-1} = C\left( \sum_{n=0}^{\infty} \frac{ A^n}{n!}\right) C^{-1} = \sum_{n=0}^{\infty} \frac{ C A^n C^{-1}}{n!}=e^{CAC^{-1}} $$

  1. $A$ диагонализируема. То есть $\exists C: CA C^{-1}$ диагональная. $$\exp {\begin{pmatrix} \lambda_1 &\dots& 0 \\ \vdots& \ddots & \vdots \\ 0& \dots & \lambda_n \end{pmatrix} } = {\begin{pmatrix} \sum_{k=0}^\infty \frac {\lambda_1^k}{k!} &\dots& 0 \\ \vdots& \ddots & \vdots \\ 0& \dots & \sum_{k=0}^\infty \frac {\lambda_n^k}{k!} \end{pmatrix} } = {\begin{pmatrix} e^\lambda_1 &\dots& 0 \\ \vdots& \ddots & \vdots \\ 0& \dots & e^\lambda_n \end{pmatrix} } $$ И теперь всё понятно, мы можем записать ответ: $X = \exp(At)\cdot X_0 = C \cdot \exp(Dt) \cdot C^{-1} \cdot X_0$
  2. Жорданова клетка.

Для этого посчитаем экспоненту от одной жордановой клетки $$\exp \begin{pmatrix} \lambda & 1 &\dots& 0 \\ \vdots& \ddots & \ddots & \vdots \\ \vdots& \dots & \lambda &1 \\ 0& \dots & 0 & \lambda \end{pmatrix} = e^{\lambda E + N} ,$$ где $N$ — это нильпотентный оператор. Но мы знаем, что $\exists k: N^k=0$, то есть $\sum_{i=0}^\infty \frac{N^i}{i!} = \sum_{i=0}^k \frac{N^i}{i!}$, и это можно посчитать! Поэтому мы получили, что $$e^{\lambda E + N} = e^\lambda E e^N,$$ помня, что $e^N$ - это многочлен.

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

Пример: $$\exp \left( \begin{pmatrix} \lambda & 1 \\ 0 & \lambda \end{pmatrix} t \right)= \exp \left( \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix} t \right) \cdot \begin{pmatrix} 1 & t \\ 0 & 1 \end{pmatrix} = \begin{pmatrix} e^{\lambda t} & t e^{\lambda t} \\ 0 & e^{\lambda t} \end{pmatrix} $$

In [3]:
axes4x4(('u','v'))
rcParams['font.family']="Arial"
title(u"Вырожденный узел", fontsize=14)
phaseportrait(lambda X, t=0: array([X[1]-X[0],-X[1]]), [[x,x] for x in linspace(-5,5,18)], [-10,10], n=90)

На рисунке выше изображен случай, когда $\lambda < 0$

А также нельзя пропустить вырожденный случай $\begin{pmatrix} \lambda &0 \\ 0 &\lambda \end{pmatrix}$. Мы его обсуждали, и картинка нам знакома (случай $\lambda>0$).

In [5]:
from itertools import product
axes4x4(('u','v'))
title(u"Дикритический узел", fontsize=14)
phaseportrait(lambda X, t=0: array([X[0],X[1]]), product(linspace(-5,5,13),[-2,0,2]), [-10,10], n=90)

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