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

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

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

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

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

Отступление: уравнения в Гамильтоновой форме

Пусть $H(x,y)$ — какая-то гладкая функция, $H(x,y):\mathbb{ R}^2 \to \mathbb{ R}$. Тогда для системы

$$\begin{cases} \dot x=\frac{\partial H}{\partial y}\\ \dot y=-\frac{\partial H}{\partial x} \end{cases}$$

функция $H(x,y)$ является первым интегралом.

Доказательство очевидно: посчитаем производную функции $H$ вдоль векторного поля $v=(\frac{\partial H}{\partial y}, -\frac{\partial H}{\partial x})$. Имеем:

$$L_v H=\frac{\partial H}{\partial y} \dot x-\frac{\partial H}{\partial x}\dot y=\frac{\partial H}{\partial y}\frac{\partial H}{\partial x}-\frac{\partial H}{\partial x}\frac{\partial H}{\partial y}=0$$

Говорят, что система записана в «гамильтоновой форме», функция $H$ называется «функцией Гамильтона». Функция Гамильтона является аналогом полной энергии — которая, как известно, сохраняется. (Слово «консервативные» в заголовке лекции означает, что что-то сохраняется.) Конечно, не всякая система, обладающая первым интегралом, является гамильтоновой.

Уравнение Ньютона

Рассмотрим уравнение:

$$\begin{equation}\tag{*} \ddot x=F(x) \end{equation}$$

Оно задает поведение частицы, движущейся в одномерном фазовом пространстве под действием силы $F$, которая зависит только от положения частицы.

In [2]:
rcParams['figure.figsize']=(5,1)
axis([-6,6,-1,2])
axis('off')
arrow( -5, 0, 10, 0, fc="k", ec="k", head_width=.5, head_length=1  )
plot([0],[0],'o',linewidth=1)
arrow( 0, 1, -1, 0, fc="k", ec="k", head_width=.5, head_length=.5)
text(-1,1.6,"$F(x)$",fontsize=16)
plt.show()

Перепишем уравнение (*) в виде системы:

$$\begin{equation}\tag{**} \begin{cases} \dot x=y\\ \dot y=F(x) \end{cases} \end{equation}$$

Введем следующие функции:

$U(x)=-\int F(x)\, dx$ — потенциальная энергия

$V(y)=\frac{y^2}{2}$ — кинетическая энергия

$E(x,y)=U(x)+V(y)$ — полная энергия

Теорема. $E$ — первый интеграл системы (**)

Доказательство. Можно заметить, что система (**) — гамильтонова. Можно посчитать производную явно: $$L_{(y,F(x))} E = \frac{\partial E}{\partial x}y+\frac{\partial E}{\partial y} F(x)=-F(x)y+yF(x)=0$$

Доказано.

Итак, полная энергия системы является первым интегралом. Это позволяет нам рисовать фазовые кривые любых уравнений вида $\ddot x = F(x)$

Примеры.

  1. $F(x)=-x$ — уравнение осциллятора, потенциальная энергия: $U(x)=\frac{x^2}{2}$, полная энергия $E(x,y)=\frac{x^2}{2}+\frac{y^2}{2}$. Построим график полной энергии:
In [3]:
from mpl_toolkits.mplot3d import Axes3D
rcParams['figure.figsize']=(7,7)
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
X, Y = meshgrid(X, Y)
Z=X**2/2+Y**2/2
fig = plt.figure()
ax = fig.gca(projection='3d')
surf1 = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, vmin=0,vmax=25,cmap=cm.coolwarm,
        linewidth=0, antialiased=False,alpha=0.7)
ax.contour(X,Y,Z, colors='black',linewidths=2,linestyles='solid',levels=[0.01,-0.5,0.5,0,-3,3,6,-6,12])
plt.show()

И его линии уровня (они же — фазовые кривые):

In [4]:
axes4x4(labels=('x','y'))
F=lambda x,y:y**2/2+x**2/2
mcontour(linspace(-5,5,300),linspace(-5,5,300),F, 
         levels=[0.01,-0.5,0.5,0,-3,3,6,-6,12],linewidths=2,
         vmin=0,vmax=25,cmap=cm.coolwarm)
vectorfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda x,y:(y,-x))

Заметим, что тут есть единственная стационарная (точнее особая) точка, где правая часть системы (**) равна нулевому вектору — это точка $(0,0)$.

2 . $F(x)=x$ — уравнение перевернутого маятника, потенциальная энергия $U(x)=-\frac{x^2}{2}$, $E(x,y)=\frac{y^2}{2}-\frac{x^2}{2}$.

In [5]:
from mpl_toolkits.mplot3d import Axes3D
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
X, Y = meshgrid(X, Y)
Z=-X**2/2+Y**2/2
fig = plt.figure(1)
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, 
        linewidth=0, antialiased=True,cmap=cm.coolwarm,vmin=-15,vmax=15,alpha=0.8)
ax.contour(X,Y,Z, colors='black',linewidths=2,linestyles='solid',levels=[-0.5,0.5,0,-3,3,6,-6])
plt.show()

И её линии уровня:

In [6]:
axes4x4(labels=('x','y'))
F=lambda x,y:y**2/2-x**2/2
mcontour(linspace(-5,5,300),linspace(-5,5,300),F, 
         levels=[-0.5,0.5,0,-3,3,6,-6],
         linewidths=2,cmap=cm.coolwarm,vmin=-15,vmax=15)
vectorfield(arange(-4,4,0.5),arange(-4,4,0.5),lambda x,y:(y,x))

Интерпретация. Точки минимума функции потенциальной энергии соответствуют устойчивым равновесиям (пример 1), точки максимума — неустойчивого (пример 2).

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

Допустим теперь, что мы начинаем с положения, находящегося левее положения равновесия ($x<0$). Если начальная скорость не слишком велика, маятник приблизится к положению равновесия, но не достигнет его, затормозит (в момент пересечения оси $x$) и станет снова падать влево. Если же начальная скорость большая, маятник перелетит через положение равновесия, и упадёт направо. Наконец, есть промежуточная траектория (она называется сепаратрисой — в данном случае, это прямолинейный луч) — вдоль неё маятник будет обречен вечно приближаться к положению равновесия, но никогда его не достигнет (по теореме существования и единственности).

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

In [10]:
rcParams['font.family']="Arial"
rcParams['figure.figsize']=(12,12)
fig = plt.figure()
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
X, Y = meshgrid(X, Y)
lev=[-0.5,0.5,0,-3,3,6,-6,12]
cmap=cm.coolwarm

data=[(X**2/2+Y**2/2,221,u'Устойчивое равновесие'),(-X**2/2+Y**2/2,222,u'Неустойчивое равновесие')]
for Z,loc,title in data:
    ax=fig.add_subplot(loc, projection='3d')
    ax.set_title(title)
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cmap,
        linewidth=0, antialiased=False,alpha=0.5,vmin=-10,vmax=10)
    ax.contour(X,Y,Z, colors='black',linewidths=2,linestyles='solid',levels=lev)
for Z, loc in zip([d[0] for d in data], [223,224]):
    ax=fig.add_subplot(loc)#FIGURE 3
    axes4x4(labels=('x','y'))
    contour(X,Y,Z, levels=lev,cmap=cmap,vmin=-10,vmax=10,linewidths=2)

plt.show()

Математический маятник

Из школьного курса физики мы знаем, что сила, действующая на математический маятник равна $F(x)=-\sin(x) mg$ (это проекция силы тяжести $mg$ на направление, касательное к окружности, которую может описывать маятник. В подходящих координатах получаем уравнение:

$$\ddot x=-\sin x$$

$$U(x)=-\cos x$$

Построим сначала функцию полной энергии, а затем ее линии уровня.

In [8]:
rcParams['figure.figsize']=(8,8)
X, Y = arange(-2, 6, 0.25), arange(-3, 3, 0.25)
X, Y = meshgrid(X, Y)
Z=Y**2/2-cos(X)
fig = plt.figure(1)
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=2, cstride=2, alpha=0.7,cmap=cm.coolwarm)
cset = ax.contourf(X, Y, Z, zdir='z', offset=-1, cmap=cm.coolwarm)
plt.show()