In [1]:
%pylab inline
from nesode import *
Populating the interactive namespace from numpy and matplotlib

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

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

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

Лекция 11. Системы линейных уравнений с постоянными коэффициентами

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

$$\tag{1} \dot z=Az,\quad z(t)\in \mathbb R^n,$$ где $A\colon \mathbb R^n\to\mathbb R^n$ — фиксированный линейный оператор.

Контрольный вопрос. Найти все особые точки уравнения (1), если $A$ — невырожденная матрица.

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

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

Лемма. Пусть $f(t)$ $f\colon \mathbb R\to \mathbb R^n$ — некоторая дифференцируемая вектор-функция, $A\colon \mathbb R^n\to\mathbb R^n$ — фиксированный линейный оператор. Тогда

$$\frac{d}{dt} (Af)(t)=A\frac{df(t)}{dt}=A\dot f(t)$$

Доказательство. По определению производной,

$$\frac{d}{dt} (Af)(t)=\lim_{h\to 0}\frac{(Af)(t+h)-(Af)(t)}{h}=…$$ в силу линейности оператора $A$ $$…=\lim_{h\to 0}\frac{A(f(t+h)-f(t))}{h}=A\lim_{h\to 0}\frac{f(t+h)-f(t)}{h}=A\frac{df(t)}{dt}$$ Что и требовалось.

Следствие. В дифференциальных линейных уравнениях можно делать замену базиса, как в обычных линейных. Действительно, пусть $w(t)$ — новая неизвестная функция, и пусть $z=Cw$, где $C$ — фиксированная невырожденная матрица. Тогда $w=C^{-1}z$ и

$$\dot w=C^{-1}\dot z=C^{-1} Az=C^{-1}AC w,$$ где первое равенство следует из леммы.

Получили такое уравнение для новой неизвестной функции $w$:

$$\dot w=C^{-1}AC w$$

Это уравнение $(1)$ в новых координатах. Видим, что линейное дифференциальное уравнение преобразуется при замене базиса точно так же, как обычное линейное уравнение. Это хорошая новость.

Линейные уравнения на плоскости: случай вещественных собственных значений

Разберем самый простой случай. Пусть $n=2$ и матрица $A$ имеет два различных вещественных собственных значения $\lambda,\mu\in \mathbb R$, $\lambda\ne\mu$. Великая наука линейная алгебра учит нас: в этом случае у матрицы $A$ есть диагонализирующий базис. То есть существует такая матрица перехода $C$, что $C^{-1}AC=D$, где

$$D=\mathop{\mathrm{diag}} (\lambda,\mu)= \begin{pmatrix}\lambda & 0\\ 0&\mu \end{pmatrix} $$

Сделаем в уравнении (1) замену $z=Cw$, $w=C^{-1} z$. Получим новое уравнение

$$\dot w=D w$$

Пусть $w=(u,v)$. Тогда новое уравнение запишется в виде системы:

$$\begin{cases} \dot u=\lambda u,\\ \dot v=\mu v \end{cases}$$

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

$$u(t)=u_0 e^{\lambda t},\quad v(t)=v_0 e^{\mu t},$$

где $(u_0, v_0)$ — начальное условие (в координатах $(u,v)$). Заметим, что уравнение автономное, и поэтому указывать конкретное $t$ не нужно. Можно ещё записать решение в таком виде:

$$w(t)= \begin{pmatrix} u(t)\\ v(t) \end{pmatrix} = \begin{pmatrix} e^{\lambda t} & 0\\ 0 & e^{\mu t} \end{pmatrix} \begin{pmatrix} u_0\\ v_0 \end{pmatrix} $$

Нам бы надо записать решение в исходных координатах, но это сделать тоже легко.

$$z(t)=Cw(t)=C\begin{pmatrix} e^{\lambda t} & 0\\ 0 & e^{\mu t} \end{pmatrix} \begin{pmatrix} u_0\\ v_0 \end{pmatrix} $$

Нам здесь даже не пришлось вычислять обратную матрицу к $C$! Впрочем, это ненадолго: если нам поставлена задача Коши в исходных координатах, то есть нужно найти решение $z(t)$ с начальным условием $z(0)=z^0=(x_0, y_0)$, то для её решения придётся выразить начальное условие через новые координаты, а тут уже потребуется $C^{-1}$:

$$z(t)=C\begin{pmatrix} e^{\lambda t} & 0\\ 0 & e^{\mu t} \end{pmatrix} C^{-1} \begin{pmatrix} x_0\\ y_0 \end{pmatrix}$$

Пример

В теории всё довольно просто. Время перейти к практике. Построим фазовый портрет, например, системы $$\dot x=2x+3y,\quad \dot y=x+4y$$ Её можно записать в виде $$ \begin{pmatrix} \dot x\\ \dot y \end{pmatrix} = \begin{pmatrix} 2&3\\ 1&4 \end{pmatrix} \begin{pmatrix} x\\ y \end{pmatrix}$$

Диагонализируем матрицу. Характеристический многочлен имеет вид $$(2-\lambda)(4-\lambda)-3=\lambda^2-6\lambda+5=(\lambda-5)(\lambda -1),$$

собственные значени $\lambda=1$ и $\mu=5$, матрица перехода к диагонализирующему базису

$$C=\begin{pmatrix} -3&1\\ 1&1 \end{pmatrix}$$

Проверим на всякий случай, что мы нигде не ошиблись, с помощью SymPy.

In [2]:
from sympy import Matrix, init_printing #нам нужны эти функции
init_printing() #это чтобы красивенько печаталось
In [3]:
A=Matrix([[2,3],[1,4]])
A.diagonalize()
Out[3]:
$$\begin{pmatrix}\left[\begin{matrix}-3 & 1\\1 & 1\end{matrix}\right], & \left[\begin{matrix}1 & 0\\0 & 5\end{matrix}\right]\end{pmatrix}$$

Это значит, что в новом базисе наша система принимает вид

$$\dot u=u,\quad \dot v=5v$$

Её фазовый портрет выглядит примерно так.

In [9]:
from itertools import product
axes4x4(labels=('u','v'))
rcParams['figure.figsize']=(6,6)
phaseportrait(lambda X, t=0: array([X[0],5*X[1]]), product(linspace(-4,4,15),[-2,2]), [-2,2], n=50)

По оси $v$ приближение к нулю происходит гораздо быстрее, чем по оси $u$ (поскольку 5>1), поэтому фазовые кривые стремятся к нулю, касаясь оси $u$. Чтобы нарисовать фазовый портрет в исходном базисе, нужно изобразить новые базисные векторы в старом базисе. Получается примерно так:

In [5]:
axes4x4(labels=('u','v'))
rcParams['figure.figsize']=(6,6)
phaseportrait(lambda X, t=0: array([2*X[0]+3*X[1],X[0]+4*X[1]]), product(linspace(-5,5,15),[-1,1]), [-2,2], n=50)
arrow(0,0,-3,1,head_width=0.2, head_length=0.3,color='red',lw=2)
arrow(0,0,1,1,head_width=0.2, head_length=0.3,color='red',lw=2)
Out[5]:
<matplotlib.patches.FancyArrow at 0x5bd9810>

Классификация линейных особых точек на плоскости

Пока мы рассматриваем случаи, когда собственные значения вещественные и различные. $$A \sim \lambda, \mu, \quad \lambda \neq \mu, \quad \lambda \mu \neq 0$$

  1. $\lambda \mu >0$
    1. $\lambda>0, \mu >0$: Неусточивый узел
    2. $\lambda<0, \mu <0$: Устойчивый узел
  2. $\lambda \mu <0$: Седло
In [11]:
fig = plt.figure()
rcParams['font.family']="Arial"
rcParams['figure.figsize']=(18,6)
fig.suptitle(u"Классификация линейных особых точек на плоскости (начало)", fontsize=20)
data=[(lambda X, t=0: array([X[0],5*X[1]]), 131, u'Неустойчивый узел'),
      (lambda X, t=0: array([-X[0],-5*X[1]]), 132, u'Устойчивый узел'),
      (lambda X, t=0: array([-X[0],X[1]]), 133, u'Седло')]
for f,loc,title in data:
    ax=fig.add_subplot(loc)
    axes4x4(('u','v'))
    ax.set_title(title,fontsize=14)
    phaseportrait(f, product(linspace(-5,5,13),[-2,0,2]), [-4,4], n=50)
plt.show()